wfdd.f90 21 KB
Newer Older
1 2 3 4 5 6 7
!
! Copyright (C) 2005  MANU/YUDONG WU/NICOLA MARZARI/ROBERTO CAR
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
8
MODULE wanpar
9
! nw:       the number of the G vectors
10 11
! nit:      the number of total iteration during searching
! nsd:      the number of steepest descent iterations
12
! ibrav:    the structure index, the same as ibrav in CP code.
13 14
  INTEGER :: nw,  nit, nsd, ibrav
  LOGICAL adapt, restart
15 16 17
! wfdt:     time step during searching
! maxwfdt:  the maximum time step during searching
! b1,b2,b3: the reciprocal lattice
18
! alat:     the lattice parameter
19
! a1,a2,a3: the real-space lattice
giannozz's avatar
giannozz committed
20 21
  real(kind=8) :: wfdt, maxwfdt, b1(3), b2(3), b3(3), alat
  real(kind=8) :: a1(3), a2(3), a3(3), tolw
22 23 24

! wfg:      the G vectors involoved in general symmetry calculation
!           the units are b1, b2, b3.
25
!           For example:
26
!            the ith G vector: wfg(i,1)*b1+wfg(i,2)*b2+wfg(i,3)*b3
27
  INTEGER, ALLOCATABLE :: wfg(:,:)
28 29

! weight:   the weight of each G vectors
30
  real(kind=8), ALLOCATABLE :: weight(:)
31
!
32
!       These are the Input variables for Damped Dynamics
33
!
34 35
! q:            imaginary mass of the Unitary Matrix
! dt:           Time Step for damped dynamics
36
! cgordd:       1=conjugate gradient/SD
37 38 39
!               any other number = damped dynamics
! fric:         damping coefficient, b/w 0 and 1
! nsteps:       Max No. of MD Steps
giannozz's avatar
giannozz committed
40
  real(kind=8) :: q, dt, fric
41
  INTEGER :: cgordd, nsteps
42 43


44
END MODULE wanpar
45 46

!----------------------------------------------------------------------
47
PROGRAM wfdd
48 49
!----------------------------------------------------------------------
!
50
!    This program works on the overlap matrix calculated
51 52 53 54 55 56 57 58
!        from parallel machine and search the unitary transformation
!        Uall corresponding to the Maximally localized Wannier functions.
!
!    The overlap matrix and lattice information are read from fort.38.
!
!
!    Searching parameters are in the input file:
!
59
!       cgordd  wfdt   maxwfdt   nit   nsd  q dt fric nsteps
60 61 62 63 64
!
!
!    The final unitary matrix Uall is output to fort.39.
!    Some running information is output to fort.24.
!
65
!                                            Yudong Wu
giannozz's avatar
giannozz committed
66 67
!                                            June 28,2001
!
68 69 70 71 72
!       This code has been modified to include Damped dynamics to
!       find the maximally localized wannier functions.
!
!                                                Manu
!                                                September 16,2001
giannozz's avatar
giannozz committed
73 74 75 76 77
!
!
!     copyright MANU/YUDONG WU/NICOLA MARZARI/ROBERTO CAR
!
!----------------------------------------------------------------------
78
  USE wanpar
giannozz's avatar
giannozz committed
79
!
80
  IMPLICIT NONE
giannozz's avatar
giannozz committed
81

82 83 84 85 86 87 88
  INTEGER :: i, j, inw, n, nspin, nupdwn(2)
  COMPLEX(kind=8), ALLOCATABLE :: O(:, :, :), Ospin(:, :, :)
  real(kind=8), ALLOCATABLE :: Uall(:,:), Uspin(:,:), u1(:,:)

        READ (5,*) cgordd, wfdt, maxwfdt, nit, nsd
        READ (5,*)  q, dt, fric, adapt, nsteps, tolw
        READ (5,*) restart
89

giannozz's avatar
giannozz committed
90 91 92 93

!
!    input the overlap matrix from fort.38
!
94 95 96 97 98 99
  REWIND 38
  READ(38, '(i5, 2i2, i3, f9.5)') n, nw, nspin, ibrav, alat
  ALLOCATE(wfg(nw, 3), weight(nw), O(nw,n,n), Uall(n,n), u1(n,n))
  IF (nspin==2) THEN
     READ(38, '(i5)') nupdwn(1)
  ENDIF
giannozz's avatar
giannozz committed
100
  nupdwn(2)=n-nupdwn(1)
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
  READ(38, *) a1
  READ(38, *) a2
  READ(38, *) a3
  READ(38, *) b1
  READ(38, *) b2
  READ(38, *) b3
  DO inw=1, nw
     READ(38, *) wfg(inw, :), weight(inw)
  ENDDO
  DO inw=1, nw
     DO i=1, n
        DO j=1, n
           READ(38, *) O(inw, i, j)
        ENDDO
     ENDDO
  ENDDO
  IF(restart) THEN
  DO i=1, n
     DO j=1, n
        READ(39, *) Uall(i, j)
     ENDDO
  ENDDO
  ELSE
giannozz's avatar
giannozz committed
124
  Uall=0.0d0
125
     DO i=1,n
giannozz's avatar
giannozz committed
126
        Uall(i,i)=1.d0
127 128
     ENDDO
  ENDIF
giannozz's avatar
giannozz committed
129

130
  REWIND 24
giannozz's avatar
giannozz committed
131

132 133 134 135
  IF(cgordd==1) THEN
    IF (nspin==1) THEN
      CALL searchwf(n, O, Uall)
    ELSE
giannozz's avatar
giannozz committed
136 137 138 139 140
!
!   For those spin-polarized calculation,
!    spin up and spin down parts are dealt with seperately
!    and the total unitary matrices are put together.
!
141 142 143 144
     WRITE(24, *) "spin up:"
     ALLOCATE(Uspin(nupdwn(1), nupdwn(1)), Ospin(nw, nupdwn(1), nupdwn(1)))
     DO i=1, nupdwn(1)
        DO j=1, nupdwn(1)
giannozz's avatar
giannozz committed
145 146
           Uspin(i, j)=Uall(i, j)
           Ospin(:, i, j)=O(:, i, j)
147 148 149 150 151
        ENDDO
     ENDDO
     CALL searchwf(nupdwn(1), Ospin, Uspin)
     DO i=1, nupdwn(1)
        DO j=1, nupdwn(1)
giannozz's avatar
giannozz committed
152
           Uall(i, j)=Uspin(i, j)
153 154 155 156 157 158 159
        ENDDO
     ENDDO
     DEALLOCATE(Uspin, Ospin)
     WRITE(24, *) "spin down:"
     ALLOCATE(Uspin(nupdwn(2), nupdwn(2)), Ospin(nw, nupdwn(2), nupdwn(2)))
     DO i=1, nupdwn(2)
        DO j=1, nupdwn(2)
giannozz's avatar
giannozz committed
160 161
           Uspin(i, j)=Uall(i+nupdwn(1), j+nupdwn(1))
           Ospin(:, i, j)=O(:, i+nupdwn(1), j+nupdwn(1))
162 163 164 165 166
        ENDDO
     ENDDO
     CALL searchwf(nupdwn(2), Ospin, Uspin)
     DO i=1, nupdwn(2)
        DO j=1, nupdwn(2)
giannozz's avatar
giannozz committed
167
           Uall(i+nupdwn(1), j+nupdwn(1))=Uspin(i, j)
168 169 170 171
        ENDDO
     ENDDO
     DEALLOCATE(Uspin, Ospin)
    ENDIF
giannozz's avatar
giannozz committed
172 173


174
  ELSE
giannozz's avatar
giannozz committed
175

176 177 178
    IF (nspin==1) THEN
      CALL ddyn(n,O,Uall)
    ELSE
giannozz's avatar
giannozz committed
179 180 181 182 183
!
!   For those spin-polarized calculation,
!    spin up and spin down parts are dealt with seperately
!    and the total unitary matrices are put together.
!
184 185 186 187
     WRITE(24, *) "spin up:"
     ALLOCATE(Uspin(nupdwn(1), nupdwn(1)), Ospin(nw, nupdwn(1), nupdwn(1)))
     DO i=1, nupdwn(1)
        DO j=1, nupdwn(1)
giannozz's avatar
giannozz committed
188 189
           Uspin(i, j)=Uall(i, j)
           Ospin(:, i, j)=O(:, i, j)
190 191 192 193 194
        ENDDO
     ENDDO
     CALL ddyn(nupdwn(1), Ospin, Uspin)
     DO i=1, nupdwn(1)
        DO j=1, nupdwn(1)
giannozz's avatar
giannozz committed
195
           Uall(i, j)=Uspin(i, j)
196 197 198 199 200 201 202
        ENDDO
     ENDDO
     DEALLOCATE(Uspin, Ospin)
     WRITE(24, *) "spin down:"
     ALLOCATE(Uspin(nupdwn(2), nupdwn(2)), Ospin(nw, nupdwn(2), nupdwn(2)))
     DO i=1, nupdwn(2)
        DO j=1, nupdwn(2)
giannozz's avatar
giannozz committed
203 204
           Uspin(i, j)=Uall(i+nupdwn(1), j+nupdwn(1))
           Ospin(:, i, j)=O(:, i+nupdwn(1), j+nupdwn(1))
205 206 207 208 209
        ENDDO
     ENDDO
     CALL ddyn(nupdwn(2), Ospin, Uspin)
     DO i=1, nupdwn(2)
        DO j=1, nupdwn(2)
giannozz's avatar
giannozz committed
210
           Uall(i+nupdwn(1), j+nupdwn(1))=Uspin(i, j)
211 212 213 214 215
        ENDDO
     ENDDO
     DEALLOCATE(Uspin, Ospin)
    ENDIF
  ENDIF
216

giannozz's avatar
giannozz committed
217 218


219 220 221 222 223 224
  REWIND 39
  DO i=1, n
     DO j=1, n
        WRITE(39, *) Uall(i, j)
     ENDDO
  ENDDO
giannozz's avatar
giannozz committed
225

226
!u1=matmul(Uall,transpose(Uall))
giannozz's avatar
giannozz committed
227 228 229 230 231 232 233

! do i=1, n
!     do j=1, n
!        write(6, *) u1(i, j)
!     end do
!  end do

234
  DEALLOCATE(wfg, weight, O, Uall,u1)
giannozz's avatar
giannozz committed
235

236
CONTAINS
giannozz's avatar
giannozz committed
237
!-------------------------------------------------------------------------
238
SUBROUTINE ddyn(m,Omat,Umat)
giannozz's avatar
giannozz committed
239 240 241 242 243 244
!    (m,m) is the size of the matrix Ospin.
!    Ospin is input overlap matrix.
!    Uspin is the output unitary transformation.
!             Rough guess for Uspin can be carried in.
!
!
245 246
!                                        MANU
!                                        SEPTEMBER 17, 2001
247 248
!-------------------------------------------------------------------------

249
  USE wanpar
250

251
  USE constants, ONLY : tpi, autoaf => BOHR_RADIUS_ANGS
252
!  implicit none
253 254 255 256 257 258 259 260
  INTEGER :: f3(nw), f4(nw), i,j,inw
  INTEGER ,INTENT(in) :: m
  real(kind=8), INTENT(inout) :: Umat(m,m)
  COMPLEX(kind=8), INTENT(inout) :: Omat(nw,m,m)
  COMPLEX(kind=8) :: U2(m,m),U3(m,m)
  INTEGER :: adjust,ini, ierr1
  real(kind=8), ALLOCATABLE, DIMENSION(:) :: wr
  real(kind=8), ALLOCATABLE, DIMENSION(:,:) :: W
giannozz's avatar
giannozz committed
261 262 263 264 265
  real(kind=8) :: t0, U(m,m), t2
  real(kind=8) :: A(m,m),oldt0,Wm(m,m),U1(m,m)
  real(kind=8) :: Aminus(m,m), Aplus(m,m),f2(3*m-2)
!  real(kind=8) :: Aminus(m,m), Aplus(m,m),f2(4*m)
  real(kind=8) :: temp(m,m)
266 267 268 269 270
  COMPLEX(kind=8) :: d(m,m), alpha, beta1, ci
  COMPLEX(kind=8) :: f1(2*m-1), wp(m*(m+1)/2),z(m,m)
  COMPLEX(kind=8), ALLOCATABLE, DIMENSION(:, :) :: X1
  COMPLEX(kind=8), ALLOCATABLE, DIMENSION(:, :, :) :: Oc
  real(kind=8) , ALLOCATABLE , DIMENSION(:) :: mt
giannozz's avatar
giannozz committed
271
  real(kind=8) :: spread, sp, oldspread
giannozz's avatar
giannozz committed
272
  real(kind=8) :: wfc(3,n), gr(nw,3)
273 274 275 276 277

  alpha=(1.d0,0.d0)
  beta1=(0.d0,0.d0)
  ci   =(0.d0,1.d0)

278 279 280
  ALLOCATE(mt(nw))
  ALLOCATE(X1(m,m))
  ALLOCATE(Oc(nw,m,m))
281 282

!   fric=friction
283
  ALLOCATE (W(m,m),wr(m))
284 285 286 287 288

!   Umat=0.d0
!   do i=1,m
!       Umat(i,i)=1.d0
!   end do
289 290

        U2=Umat*alpha
291 292 293 294

!
! update Oc using the initial guess of Uspin
!
295
  DO inw=1, nw
296 297
    X1(:, :)=Omat(inw, :, :)
     U3=beta1
298 299
!    call ZGEMUL(U2, m, 'T', X1, m, 'N', U3, m, m,m,m)
     CALL zgemm ('T', 'N', m,m,m,alpha,U2,m,X1,m,beta1,U3,m)
300
    X1=beta1
301 302
!    call ZGEMUL(U3, m, 'N', U2, m, 'N', X1, m, m,m,m)
     CALL zgemm ('N','N', m,m,m,alpha,U3,m,U2,m,beta1,X1,m)
303
    Oc(inw, :, :)=X1(:, :)
304
  ENDDO
305

306 307
        U2=beta1
        U3=beta1
308

giannozz's avatar
giannozz committed
309
 oldspread=0.0d0
310 311 312 313 314
  WRITE(24, *) "spread: (unit \AA^2)"
  DO i=1, m
     mt=1.d0-dble(Oc(:,i,i)*conjg(Oc(:,i,i)))
     sp= (alat*autoaf/tpi)**2*sum(mt*weight)
     WRITE(24, '(f10.7)') (alat*autoaf/tpi)**2*sum(mt*weight)
315
     oldspread=oldspread+sp
316
  ENDDO
317 318

  oldspread=oldspread/m
319
     WRITE(51, '(f10.7)') oldspread
320 321 322 323 324 325 326

    oldt0=0.d0
    A=0.d0
    Aminus=A
    temp=Aminus


327
!        START ITERATIONS HERE
328

329
  DO ini=1, nsteps
330 331

    t0=0.d0     !use t0 to store the value of omega
332 333 334 335 336
    DO inw=1, nw
       DO i=1, m
          t0=t0+dble(conjg(Oc(inw, i, i))*Oc(inw, i, i))
       ENDDO
    ENDDO
337

338
    WRITE(6, *) t0
339 340


341 342 343 344
        IF(abs(t0-oldt0)<tolw) THEN
           WRITE(6,*) "MLWF Generated at Step",ini
           GOTO 241
        ENDIF
345

346 347
        IF(adapt) THEN
          IF(oldt0<t0) THEN
giannozz's avatar
giannozz committed
348
            fric=fric/2.d0
349 350
            A=Aminus
            Aminus=temp
351 352
          ENDIF
        ENDIF
353 354 355 356 357 358

!   calculate d(omega)/dA and store result in W
!   this is the force for the damped dynamics
!

    W=0.d0
359
    DO inw=1, nw
360
       t2=weight(inw)
361 362 363
       DO i=1,m
          DO j=1,m
             W(i,j)=W(i,j)+t2*dble(Oc(inw,i,j)*conjg(Oc(inw,i,i)        &
giannozz's avatar
giannozz committed
364
                  -Oc(inw,j,j))+conjg(Oc(inw,j,i))*(Oc(inw,i,i)-Oc(inw,j,j)))
365 366 367
          ENDDO
       ENDDO
    ENDDO
368

369 370

!   the verlet scheme to calculate A(t+dt)
371 372

        Aplus=0.d0
373

374 375
   DO i=1,m
     DO j=i+1,m
376
         Aplus(i,j)=Aplus(i,j)+(2*dt/(2*dt+fric))*(2*A(i,j)               &
377
         -Aminus(i,j)+(dt*dt/q)*W(i,j)) + (fric/(2*dt+fric))*Aminus(i,j)
378 379
     ENDDO
   ENDDO
380

381 382
        Aplus=Aplus-transpose(Aplus)
        Aplus=(Aplus-A)
383

384 385 386 387 388
    DO i=1, m
       DO j=i,m
        wp(i + (j-1)*j/2) = cmplx(0.0d0, Aplus(i,j), kind=8)
       ENDDO
    ENDDO
389

390
    CALL zhpev('V','U',m,wp,wr,z,m,f1,f2,ierr1)
391 392
!    call zhpev(21, wp, wr, z, m, m, f2, 4*m)

393 394 395 396
    IF (ierr1/=0) THEN
   WRITE(6,*) "failed to diagonalize W!"
    STOP
    ENDIF
397 398

    d=0.d0
399
    DO i=1, m
400
       d(i, i)=exp(ci*wr(i)*dt)
401
    ENDDO      !d=exp(d)
402 403

!   U=z*exp(d)*z+
404
!
405
     U3=beta1
406
     CALL zgemm ('N', 'N', m,m,m,alpha,z,m,d,m,beta1,U3,m)
407
     U2=beta1
408 409
     CALL zgemm ('N','c', m,m,m,alpha,U3,m,z,m,beta1,U2,m)
    U=dble(U2)
410 411 412 413 414 415 416 417 418 419 420
    U2=beta1
    U3=beta1

   temp=Aminus
   Aminus=A
   A=Aplus


!   update Umat
!
    U1=beta1
421 422
    CALL dgemm ('N', 'N', m,m,m,alpha,Umat,m,U,m,beta1,U1,m)
    Umat=U1
423 424 425 426 427

!   update Oc
!
    U2=Umat*alpha
    U3=beta1
428
  DO inw=1, nw
429
    X1(:, :)=Omat(inw, :, :)
430
    CALL zgemm ('T', 'N', m,m,m,alpha,U2,m,X1,m,beta1,U3,m)
431
    X1=beta1
432
    CALL zgemm ('N','N',m,m,m,alpha,U3,m,U2,m,beta1,X1,m)
433
    Oc(inw, :, :)=X1(:, :)
434
  ENDDO
435 436 437
    U2=beta1
    U3=beta1

438 439 440
        IF(abs(t0-oldt0)>=tolw.and.ini>=nsteps) THEN
        GOTO 241
        ENDIF
441 442 443

    oldt0=t0

444
   ENDDO
giannozz's avatar
giannozz committed
445
241  spread=0.0d0
446 447 448 449 450
  WRITE(24, *) "spread: (unit \AA^2)"
  DO i=1, m
     mt=1.d0-dble(Oc(:,i,i)*conjg(Oc(:,i,i)))
     sp= (alat*autoaf/tpi)**2*sum(mt*weight)
     WRITE(24, '(f10.7)') (alat*autoaf/tpi)**2*sum(mt*weight)
451
     spread=spread+sp
452
  ENDDO
453 454

  spread=spread/m
455
     WRITE(51, '(f10.7)') spread
456 457


458
  DEALLOCATE(wr, W)
459 460 461 462


!   output wfc's and spreads of the max. loc. wf's
!
463 464
  ALLOCATE(wr(nw), W(nw, nw))
  DO inw=1, nw
465
     gr(inw, :)=wfg(inw,1)*b1(:)+wfg(inw,2)*b2(:)+wfg(inw,3)*b3(:)
466
  ENDDO
467
!
degironc's avatar
degironc committed
468
! set up a matrix with the element (i,j) is G_i * G_j * weight(j)
469 470
! to check the correctness of choices on G vectors
!
471 472 473
  DO i=1, nw
     DO j=1, nw
        W(i,j)=sum(gr(i,:)*gr(j,:))*weight(j)
474
!        write(6, *) i,j,W(i,j)
475 476
     ENDDO
  ENDDO
477
!  write(24, *) "wannier function centers: (unit:\AA)"
478
  DO i=1, m
giannozz's avatar
giannozz committed
479
     mt=-aimag(log(Oc(:,i,i)))/tpi
480 481 482 483 484 485
     wfc(1, i)=sum(mt*weight*gr(:,1))
     wfc(2, i)=sum(mt*weight*gr(:,2))
     wfc(3, i)=sum(mt*weight*gr(:,3))
     DO inw=1, nw
        wr(inw)=sum(wfc(:,i)*gr(inw,:))-mt(inw)
     ENDDO
486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502
     mt=wr
     f3=0
     adjust=0
!
!   balance the phase factor if necessary
!
!     do while(SUM((mt-f3)**2).gt.0.01d0)
!        f4=f3
!        f3=nint(mt-mt(1))
!        if (adjust.gt.200) f3=f3-1
!        if (adjust.gt.100.and.adjust.le.200) f3=f3+1
!        mt=wr+matmul(W, f3)
!        write(6,*) "mt:", mt
!        write(6,*) "f3:", f3
!        adjust=adjust+1
!        if (adjust.gt.300) stop "unable to balance the phase!"
!     end do
503 504 505 506 507
     wfc(1,i)=(wfc(1,i)+sum(mt*weight*gr(:,1)))*alat
     wfc(2,i)=(wfc(2,i)+sum(mt*weight*gr(:,2)))*alat
     wfc(3,i)=(wfc(3,i)+sum(mt*weight*gr(:,3)))*alat
  ENDDO

508 509 510 511 512 513 514
!  if (ibrav.eq.1.or.ibrav.eq.6.or.ibrav.eq.8) then
!     do i=1, m
!        if (wfc(1, i).lt.0) wfc(1, i)=wfc(1, i)+a1(1)
!        if (wfc(2, i).lt.0) wfc(2, i)=wfc(2, i)+a2(2)
!        if (wfc(3, i).lt.0) wfc(3, i)=wfc(3, i)+a3(3)
!     end do
!  end if
515 516 517 518 519 520 521
  DO i=1, m
     WRITE(26, '(3f11.6)') wfc(:,i)*autoaf
  ENDDO

     WRITE(6,*) "Friction =", fric
     WRITE(6,*) "Mass =", q

522

523
  DEALLOCATE(wr, W)
524

525
  RETURN
526

527
  END SUBROUTINE ddyn
528 529

!-----------------------------------------------------------------------
530
  SUBROUTINE searchwf(m, Omat, Umat)
531 532 533 534 535 536
!-----------------------------------------------------------------------
!    (m,m) is the size of the matrix Ospin.
!    Ospin is input overlap matrix.
!    Uspin is the output unitary transformation.
!             Rough guess for Uspin can be carried in.
!
537 538
  USE wanpar
  USE constants, ONLY : tpi, autoaf => BOHR_RADIUS_ANGS
539 540 541 542
!
!
!     conjugated gradient to search maximization
!
543
  IMPLICIT NONE
544
!
545 546 547
  INTEGER, INTENT(in) :: m
  COMPLEX(kind=8), INTENT(in) :: Omat(nw, m, m)
  real(kind=8), INTENT(inout) :: Umat(m,m)
548
!
549 550
  INTEGER :: i, j, k, l, ig, ierr, ti, tj, tk, inw, ir
  INTEGER ::  istep
giannozz's avatar
giannozz committed
551
  real(kind=8) :: slope, slope2, t1, t2, t3, mt(nw),t21,temp1,maxdt
giannozz's avatar
giannozz committed
552 553
  real(kind=8) :: U(m,m), wfc(3, m), Wm(m,m), schd(m,m), f2(3*m-2), gr(nw, 3)
  real(kind=8) :: Uspin2(m,m),temp2,wfdtold,oldt1,t01, d3(m,m), d4(m,m), U1(m,m)
554 555 556 557 558
  real(kind=8), ALLOCATABLE, DIMENSION(:) :: wr
  real(kind=8), ALLOCATABLE, DIMENSION(:,:) :: W
  COMPLEX(kind=8) :: ci, ct1, ct2, ct3, z(m, m), X(m, m), d(m,m), d2(m,m)
  COMPLEX(kind=8) :: f1(2*m-1), wp(m*(m+1)/2), Oc(nw, m, m), alpha, beta1
  COMPLEX(kind=8) ::  Oc2(nw, m, m),wp1(m*(m+1)/2), X1(m,m), U2(m,m), U3(m,m)
559 560 561 562 563 564

!
  ci=(0.d0,1.d0)
  alpha=(1.0d0, 0.0d0)
  beta1=(0.0d0, 0.0d0)
!
565
  ALLOCATE(W(m,m), wr(m))
566 567 568 569 570 571 572 573 574 575 576 577 578 579


!  Umat=0.d0
!  do i=1,m
!    Umat(i,i)=1.d0
!  end do
  Oc=beta1
  Oc2=beta1
  X1=beta1
  U2=Umat*alpha

!
! update Oc using the initial guess of Uspin
!
580
  DO inw=1, nw
581 582
     X1(:, :)=Omat(inw, :, :)
     U3=beta1
583
     CALL zgemm ('T', 'N', m,m,m,alpha,U2,m,X1,m,beta1,U3,m)
584
     X1=beta1
585
     CALL zgemm ('N','N', m,m,m,alpha,U3,m,U2,m,beta1,X1,m)
586
     Oc(inw, :, :)=X1(:, :)
587
  ENDDO
588 589 590 591 592 593 594 595 596

     U2=beta1
     U3=beta1

  W=0.d0
  schd=0.d0
  oldt1=0.d0
  wfdtold=0.d0

597
  DO k=1, nit
598
    t01=0.d0     !use t1 to store the value of omiga
599 600 601 602 603
    DO inw=1, nw
       DO i=1, m
          t01=t01+dble(conjg(Oc(inw, i, i))*Oc(inw, i, i))
       ENDDO
    ENDDO
604

605 606 607
    WRITE(6,*) t01

    IF(abs(oldt1-t01)<tolw) GOTO 40
608 609

    oldt1=t01
610

611 612 613 614 615
!   calculate d(omiga)/dW and store result in W
!   W should be a real symmetric matrix for gamma-point calculation
!
    Wm=W
    W=0.d0
616
    DO inw=1, nw
617
       t2=weight(inw)
618 619 620
       DO i=1,m
          DO j=i+1,m
             W(i,j)=W(i,j)+t2*dble(Oc(inw,i,j)*conjg(Oc(inw,i,i)        &
giannozz's avatar
giannozz committed
621
              -Oc(inw,j,j))+conjg(Oc(inw,j,i))*(Oc(inw,i,i)-Oc(inw,j,j)))
622 623 624
          ENDDO
       ENDDO
    ENDDO
625
    W=W-transpose(W)
626

627
!   calculate slope=d(omiga)/d(lamda)
628 629
    slope=sum(W**2)

630 631
!   calculate slope2=d2(omiga)/d(lamda)2
    slope2=0.d0
632 633 634
    DO ti=1, m
       DO tj=1, m
          DO tk=1, m
635
             t2=0.d0
636 637
             DO inw=1, nw
                t2=t2+dble(Oc(inw,tj,tk)*conjg(Oc(inw,tj,tj)+Oc(inw,tk,tk) &
638
                          -2.d0*Oc(inw,ti,ti))-4.d0*Oc(inw,ti,tk)          &
giannozz's avatar
giannozz committed
639
                          *conjg(Oc(inw,ti,tj)))*weight(inw)
640
             ENDDO
641
             slope2=slope2+W(tk,ti)*W(ti,tj)*2.d0*t2
642 643 644
          ENDDO
       ENDDO
     ENDDO
645
    slope2=2.d0*slope2
646

647
!   use parabola approximation. Defined by 1 point and 2 slopes
648 649 650 651
    IF (slope2<0) wfdt=-slope/2.d0/slope2
    IF (maxwfdt>0.and.wfdt>maxwfdt) wfdt=maxwfdt

    IF (k<nsd) THEN
652 653 654
       schd=W    !use steepest-descent technique

!   calculate slope=d(omiga)/d(lamda)
655
    slope=sum(schd**2)
656 657

!       schd=schd*maxwfdt
658 659 660 661 662
    DO i=1, m
       DO j=i, m
        wp1(i + (j-1)*j/2) = cmplx(0.0d0, schd(i,j), kind=8)
       ENDDO
    ENDDO
663

664 665
    CALL zhpev('V','U',m,wp1,wr,z,m,f1,f2,ierr)
    IF (ierr/=0) STOP 'failed to diagonalize W!'
666

667
    ELSE
668 669 670 671 672 673
!
!     construct conjugated gradient
!        d(i)=-g(i)+belta(i)*d(i-1)
!        belta^FR(i)=g(i)t*g(i)/(g(i-1)t*g(i-1))
!        belta^PR(i)=g(i)t*(g(i)-g(i-1))/(g(i-1)t*g(i-1))
!
674 675
        CALL dgemm ('T','N', m,m,m,alpha,Wm,m,Wm,m,beta1,d3,m)

676 677

       t1=0.d0
678
       DO i=1, m
679
          t1=t1+d3(i, i)
680 681
       ENDDO
       IF (t1/=0) THEN
682
          d4=(W-Wm)
683
          CALL dgemm ('T','N', m,m,m,alpha,W,m,d4,m,beta1,d3,m)
684
          t2=0.d0
685
          DO i=1, m
686
             t2=t2+d3(i, i)
687
          ENDDO
688 689
          t3=t2/t1
          schd=W+schd*t3
690
       ELSE
691
          schd=W
692
       ENDIF
693
!
694 695
!        calculate the new d(Lambda) for the new Search Direction
!        added by Manu. September 19, 2001
696 697
!
!   calculate slope=d(omiga)/d(lamda)
698
    slope=sum(schd**2)
699
!------------------------------------------------------------------------
700
!   schd=schd*maxwfdt
701 702 703 704 705
    DO i=1, m
       DO j=i, m
        wp1(i + (j-1)*j/2) = cmplx(0.0d0, schd(i,j), kind=8)
       ENDDO
    ENDDO
706

707 708
    CALL zhpev('V','U',m,wp1,wr,z,m,f1,f2,ierr)
    IF (ierr/=0) STOP 'failed to diagonalize W!'
709 710 711 712

      maxdt=maxwfdt

11    d=0.d0
713
    DO i=1, m
714
       d(i, i)=exp(ci*(maxwfdt)*wr(i))
715 716
    ENDDO      !d=exp(d)

717 718
!   U=z*exp(d)*z+
     U3=beta1
719
     CALL zgemm ('N', 'N', m,m,m,alpha,z,m,d,m,beta1,U3,m)
720
     U2=beta1
721 722
     CALL zgemm ('N','c', m,m,m,alpha,U3,m,z,m,beta1,U2,m)
     U=dble(U2)
723 724 725 726 727
     U2=beta1
     U3=beta1
!
!   update Uspin
    U1=beta1
728
    CALL dgemm ('N', 'N', m,m,m,alpha,Umat,m,U,m,beta1,U1,m)
729
    Umat=U1
730

731 732 733 734 735 736
!    Uspin2=matmul(Uspin, U2)
!
!   update Oc
!
     U2=Umat*alpha
     U3=beta1
737
     DO inw=1, nw
738
      X1(:,:)=Omat(inw,:,:)
739
      CALL zgemm ('T', 'N', m,m,m,alpha,U2,m,X1,m,beta1,U3,m)
740
      X1=beta1
741
      CALL zgemm ('N','N',m,m,m,alpha,U3,m,U2,m,beta1,X1,m)
742
      Oc2(inw, :,:)=X(:,:)
743 744
     ENDDO
     U2=beta1
745 746 747
     U3=beta1
!
    t21=0.d0     !use t21 to store the value of omiga
748 749 750 751 752 753
    DO inw=1, nw
       DO i=1, m
          t21=t21+dble(conjg(Oc2(inw, i, i))*Oc2(inw, i, i))
       ENDDO
    ENDDO

754 755 756 757
      temp1=-((t01-t21)+slope*maxwfdt)/(maxwfdt**2)
      temp2=slope
      wfdt=-temp2/(2*temp1)

758
        IF (wfdt>maxwfdt.or.wfdt<0.d0) THEN
759
        maxwfdt=2*maxwfdt
760 761
        GOTO 11
        ENDIF
762 763 764 765 766 767 768 769 770 771 772 773 774

        maxwfdt=maxdt
!
!
!   use parabola approximation. Defined by 2 point and 1 slopes
!    if (slope2.lt.0) wfdt=-slope/2.d0/slope2
!    if (maxwfdt.gt.0.and.wfdt.gt.maxwfdt) wfdt=maxwfdt
!
!    write(6, '(e12.5E2,1x,e11.5E2,1x,f6.2)') slope2, slope, wfdt
!-------------------------------------------------------------------------
!
!      schd is the new searching direction
!
775 776
    ENDIF

777
    d=0.d0
778
    DO i=1, m
779
       d(i, i)=exp(ci*wfdt*wr(i))
780 781
    ENDDO          !d=exp(d)

782 783 784 785

!   U=z*exp(d)*z+
!
     U3=beta1
786
     CALL zgemm ('N', 'N', m,m,m,alpha,z,m,d,m,beta1,U3,m)
787
     U2=beta1
788 789
     CALL zgemm ('N','c', m,m,m,alpha,U3,m,z,m,beta1,U2,m)
     U=dble(U2)
790 791 792 793 794 795
     U2=beta1
     U3=beta1

!   update Uspin
!
    U1=beta1
796
    CALL dgemm ('N', 'N', m,m,m,alpha,Umat,m,U,m,beta1,U1,m)
797 798 799 800 801 802
    Umat=U1

!   update Oc
!
       U2=Umat*alpha
       U3=beta1
803
     DO inw=1, nw
804
       X1(:, :)=Omat(inw, :, :)
805
       CALL zgemm ('T', 'N', m,m,m,alpha,U2,m,X1,m,beta1,U3,m)
806
       X1=beta1
807
       CALL zgemm ('N','N',m,m,m,alpha,U3,m,U2,m,beta1,X1,m)
808
       Oc(inw, :, :)=X1(:, :)
809
     ENDDO
810 811
    U2=beta1
    U3=beta1
812
  ENDDO
813

814
40  DEALLOCATE(W, wr)
815 816 817 818

!
! calculate the spread
!
819 820 821 822 823
  WRITE(24, *) "spread: (unit \AA^2)"
  DO i=1, m
     mt=1.d0-dble(Oc(:,i,i)*conjg(Oc(:,i,i)))
     WRITE(24, '(f10.7)') (alat*autoaf/tpi)**2*sum(mt*weight)
  ENDDO
824 825 826 827

!
! calculate wannier-function centers
!
828 829
  ALLOCATE(wr(nw), W(nw, nw))
  DO inw=1, nw
830
     gr(inw, :)=wfg(inw,1)*b1(:)+wfg(inw,2)*b2(:)+wfg(inw,3)*b3(:)
831
  ENDDO
832
!
degironc's avatar
degironc committed
833
! set up a matrix with the element (i,j) is G_i * G_j * weight(j)
834 835
! to check the correctness of choices on G vectors
!
836 837 838
  DO i=1, nw
     DO j=1, nw
        W(i,j)=sum(gr(i,:)*gr(j,:))*weight(j)
839
!        write(6, *) i,j,W(i,j)
840 841
     ENDDO
  ENDDO
842
! write(24, *) "wannier function centers: (unit:\AA)"
843
  DO i=1, m
giannozz's avatar
giannozz committed
844
     mt=-aimag(log(Oc(:,i,i)))/tpi
845 846 847 848 849 850
     wfc(1, i)=sum(mt*weight*gr(:,1))
     wfc(2, i)=sum(mt*weight*gr(:,2))
     wfc(3, i)=sum(mt*weight*gr(:,3))
     DO inw=1, nw
        wr(inw)=sum(wfc(:,i)*gr(inw,:))-mt(inw)
     ENDDO
851
     mt=wr
852 853 854 855
     wfc(1, i)=(wfc(1,i)+sum(mt*weight*gr(:,1)))*alat
     wfc(2, i)=(wfc(2,i)+sum(mt*weight*gr(:,2)))*alat
     wfc(3, i)=(wfc(3,i)+sum(mt*weight*gr(:,3)))*alat
  ENDDO
856 857


858 859 860 861 862 863
  DO i=1, m
     WRITE(26, '(3f11.6)') wfc(:,i)*autoaf
  ENDDO
  DEALLOCATE(wr, W)
  RETURN
  END SUBROUTINE searchwf
864 865


866
 END PROGRAM wfdd