1!
2! Copyright (C) 2005  MANU/YUDONG WU/NICOLA MARZARI/ROBERTO CAR
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8MODULE wanpar
9! nw:       the number of the G vectors
10! nit:      the number of total iteration during searching
11! nsd:      the number of steepest descent iterations
12! ibrav:    the structure index, the same as ibrav in CP code.
13  INTEGER :: nw,  nit, nsd, ibrav
14  LOGICAL adapt, restart
15! wfdt:     time step during searching
16! maxwfdt:  the maximum time step during searching
17! b1,b2,b3: the reciprocal lattice
18! alat:     the lattice parameter
19! a1,a2,a3: the real-space lattice
20  real(kind=8) :: wfdt, maxwfdt, b1(3), b2(3), b3(3), alat
21  real(kind=8) :: a1(3), a2(3), a3(3), tolw
22
23! wfg:      the G vectors involoved in general symmetry calculation
24!           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! q:            imaginary mass of the Unitary Matrix
35! dt:           Time Step for damped dynamics
36! cgordd:       1=conjugate gradient/SD
37!               any other number = damped dynamics
38! fric:         damping coefficient, b/w 0 and 1
39! nsteps:       Max No. of MD Steps
40  real(kind=8) :: q, dt, fric
41  INTEGER :: cgordd, nsteps
42
43
44END MODULE wanpar
45
46!----------------------------------------------------------------------
47PROGRAM wfdd
48!----------------------------------------------------------------------
49!
50!    This program works on the overlap matrix calculated
51!        from parallel machine and search the unitary transformation
52!        Uall corresponding to the Maximally localized Wannier functions.
53!
54!    The overlap matrix and lattice information are read from fort.38.
55!
56!
57!    Searching parameters are in the input file:
58!
59!       cgordd  wfdt   maxwfdt   nit   nsd  q dt fric nsteps
60!
61!
62!    The final unitary matrix Uall is output to fort.39.
63!    Some running information is output to fort.24.
64!
65!                                            Yudong Wu
66!                                            June 28,2001
67!
68!       This code has been modified to include Damped dynamics to
69!       find the maximally localized wannier functions.
70!
71!                                                Manu
72!                                                September 16,2001
73!
74!
75!     copyright MANU/YUDONG WU/NICOLA MARZARI/ROBERTO CAR
76!
77!----------------------------------------------------------------------
78  USE wanpar
79!
80  IMPLICIT NONE
81
82  INTEGER :: i, j, inw, n, nspin, nupdwn(2)
83  COMPLEX(kind=8), ALLOCATABLE :: O(:, :, :), Ospin(:, :, :)
84  real(kind=8), ALLOCATABLE :: Uall(:,:), Uspin(:,:), u1(:,:)
85
86        READ (5,*) cgordd, wfdt, maxwfdt, nit, nsd
87        READ (5,*)  q, dt, fric, adapt, nsteps, tolw
88        READ (5,*) restart
89
90
91!
92!    input the overlap matrix from fort.38
93!
94  REWIND 38
95  READ(38, '(i5, 2i2, i3, f9.5)') n, nw, nspin, ibrav, alat
96  ALLOCATE(wfg(nw, 3), weight(nw), O(nw,n,n), Uall(n,n), u1(n,n))
97  IF (nspin==2) THEN
98     READ(38, '(i5)') nupdwn(1)
99  ENDIF
100  nupdwn(2)=n-nupdwn(1)
101  READ(38, *) a1
102  READ(38, *) a2
103  READ(38, *) a3
104  READ(38, *) b1
105  READ(38, *) b2
106  READ(38, *) b3
107  DO inw=1, nw
108     READ(38, *) wfg(inw, :), weight(inw)
109  ENDDO
110  DO inw=1, nw
111     DO i=1, n
112        DO j=1, n
113           READ(38, *) O(inw, i, j)
114        ENDDO
115     ENDDO
116  ENDDO
117  IF(restart) THEN
118  DO i=1, n
119     DO j=1, n
120        READ(39, *) Uall(i, j)
121     ENDDO
122  ENDDO
123  ELSE
124  Uall=0.0d0
125     DO i=1,n
126        Uall(i,i)=1.d0
127     ENDDO
128  ENDIF
129
130  REWIND 24
131
132  IF(cgordd==1) THEN
133    IF (nspin==1) THEN
134      CALL searchwf(n, O, Uall)
135    ELSE
136!
137!   For those spin-polarized calculation,
138!    spin up and spin down parts are dealt with seperately
139!    and the total unitary matrices are put together.
140!
141     WRITE(24, *) "spin up:"
142     ALLOCATE(Uspin(nupdwn(1), nupdwn(1)), Ospin(nw, nupdwn(1), nupdwn(1)))
143     DO i=1, nupdwn(1)
144        DO j=1, nupdwn(1)
145           Uspin(i, j)=Uall(i, j)
146           Ospin(:, i, j)=O(:, i, j)
147        ENDDO
148     ENDDO
149     CALL searchwf(nupdwn(1), Ospin, Uspin)
150     DO i=1, nupdwn(1)
151        DO j=1, nupdwn(1)
152           Uall(i, j)=Uspin(i, j)
153        ENDDO
154     ENDDO
155     DEALLOCATE(Uspin, Ospin)
156     WRITE(24, *) "spin down:"
157     ALLOCATE(Uspin(nupdwn(2), nupdwn(2)), Ospin(nw, nupdwn(2), nupdwn(2)))
158     DO i=1, nupdwn(2)
159        DO j=1, nupdwn(2)
160           Uspin(i, j)=Uall(i+nupdwn(1), j+nupdwn(1))
161           Ospin(:, i, j)=O(:, i+nupdwn(1), j+nupdwn(1))
162        ENDDO
163     ENDDO
164     CALL searchwf(nupdwn(2), Ospin, Uspin)
165     DO i=1, nupdwn(2)
166        DO j=1, nupdwn(2)
167           Uall(i+nupdwn(1), j+nupdwn(1))=Uspin(i, j)
168        ENDDO
169     ENDDO
170     DEALLOCATE(Uspin, Ospin)
171    ENDIF
172
173
174  ELSE
175
176    IF (nspin==1) THEN
177      CALL ddyn(n,O,Uall)
178    ELSE
179!
180!   For those spin-polarized calculation,
181!    spin up and spin down parts are dealt with seperately
182!    and the total unitary matrices are put together.
183!
184     WRITE(24, *) "spin up:"
185     ALLOCATE(Uspin(nupdwn(1), nupdwn(1)), Ospin(nw, nupdwn(1), nupdwn(1)))
186     DO i=1, nupdwn(1)
187        DO j=1, nupdwn(1)
188           Uspin(i, j)=Uall(i, j)
189           Ospin(:, i, j)=O(:, i, j)
190        ENDDO
191     ENDDO
192     CALL ddyn(nupdwn(1), Ospin, Uspin)
193     DO i=1, nupdwn(1)
194        DO j=1, nupdwn(1)
195           Uall(i, j)=Uspin(i, j)
196        ENDDO
197     ENDDO
198     DEALLOCATE(Uspin, Ospin)
199     WRITE(24, *) "spin down:"
200     ALLOCATE(Uspin(nupdwn(2), nupdwn(2)), Ospin(nw, nupdwn(2), nupdwn(2)))
201     DO i=1, nupdwn(2)
202        DO j=1, nupdwn(2)
203           Uspin(i, j)=Uall(i+nupdwn(1), j+nupdwn(1))
204           Ospin(:, i, j)=O(:, i+nupdwn(1), j+nupdwn(1))
205        ENDDO
206     ENDDO
207     CALL ddyn(nupdwn(2), Ospin, Uspin)
208     DO i=1, nupdwn(2)
209        DO j=1, nupdwn(2)
210           Uall(i+nupdwn(1), j+nupdwn(1))=Uspin(i, j)
211        ENDDO
212     ENDDO
213     DEALLOCATE(Uspin, Ospin)
214    ENDIF
215  ENDIF
216
217
218
219  REWIND 39
220  DO i=1, n
221     DO j=1, n
222        WRITE(39, *) Uall(i, j)
223     ENDDO
224  ENDDO
225
226!u1=matmul(Uall,transpose(Uall))
227
228! do i=1, n
229!     do j=1, n
230!        write(6, *) u1(i, j)
231!     end do
232!  end do
233
234  DEALLOCATE(wfg, weight, O, Uall,u1)
235
236CONTAINS
237!-------------------------------------------------------------------------
238SUBROUTINE ddyn(m,Omat,Umat)
239!    (m,m) is the size of the matrix Ospin.
240!    Ospin is input overlap matrix.
241!    Uspin is the output unitary transformation.
242!             Rough guess for Uspin can be carried in.
243!
244!
245!                                        MANU
246!                                        SEPTEMBER 17, 2001
247!-------------------------------------------------------------------------
248
249  USE wanpar
250
251  USE constants, ONLY : tpi, autoaf => BOHR_RADIUS_ANGS
252!  implicit none
253  INTEGER :: f3(nw), f4(nw), i,j,inw
254  INTEGER ,INTENT(in) :: m
255  real(kind=8), INTENT(inout) :: Umat(m,m)
256  COMPLEX(kind=8), INTENT(inout) :: Omat(nw,m,m)
257  COMPLEX(kind=8) :: U2(m,m),U3(m,m)
258  INTEGER :: adjust,ini, ierr1
259  real(kind=8), ALLOCATABLE, DIMENSION(:) :: wr
260  real(kind=8), ALLOCATABLE, DIMENSION(:,:) :: W
261  real(kind=8) :: t0, U(m,m), t2
262  real(kind=8) :: A(m,m),oldt0,Wm(m,m),U1(m,m)
263  real(kind=8) :: Aminus(m,m), Aplus(m,m),f2(3*m-2)
264!  real(kind=8) :: Aminus(m,m), Aplus(m,m),f2(4*m)
265  real(kind=8) :: temp(m,m)
266  COMPLEX(kind=8) :: d(m,m), alpha, beta1, ci
267  COMPLEX(kind=8) :: f1(2*m-1), wp(m*(m+1)/2),z(m,m)
268  COMPLEX(kind=8), ALLOCATABLE, DIMENSION(:, :) :: X1
269  COMPLEX(kind=8), ALLOCATABLE, DIMENSION(:, :, :) :: Oc
270  real(kind=8) , ALLOCATABLE , DIMENSION(:) :: mt
271  real(kind=8) :: spread, sp, oldspread
272  real(kind=8) :: wfc(3,n), gr(nw,3)
273
274  alpha=(1.d0,0.d0)
275  beta1=(0.d0,0.d0)
276  ci   =(0.d0,1.d0)
277
278  ALLOCATE(mt(nw))
279  ALLOCATE(X1(m,m))
280  ALLOCATE(Oc(nw,m,m))
281
282!   fric=friction
283  ALLOCATE (W(m,m),wr(m))
284
285!   Umat=0.d0
286!   do i=1,m
287!       Umat(i,i)=1.d0
288!   end do
289
290        U2=Umat*alpha
291
292!
293! update Oc using the initial guess of Uspin
294!
295  DO inw=1, nw
296    X1(:, :)=Omat(inw, :, :)
297     U3=beta1
298!    call ZGEMUL(U2, m, 'T', X1, m, 'N', U3, m, m,m,m)
299     CALL zgemm ('T', 'N', m,m,m,alpha,U2,m,X1,m,beta1,U3,m)
300    X1=beta1
301!    call ZGEMUL(U3, m, 'N', U2, m, 'N', X1, m, m,m,m)
302     CALL zgemm ('N','N', m,m,m,alpha,U3,m,U2,m,beta1,X1,m)
303    Oc(inw, :, :)=X1(:, :)
304  ENDDO
305
306        U2=beta1
307        U3=beta1
308
309 oldspread=0.0d0
310  WRITE(24, *) "spread: (unit \AA^2)"
311  DO i=1, m
312     mt=1.d0-dble(Oc(:,i,i)*conjg(Oc(:,i,i)))
313     sp= (alat*autoaf/tpi)**2*sum(mt*weight)
314     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    oldt0=0.d0
322    A=0.d0
323    Aminus=A
324    temp=Aminus
325
326
327!        START ITERATIONS HERE
328
329  DO ini=1, nsteps
330
331    t0=0.d0     !use t0 to store the value of omega
332    DO inw=1, nw
333       DO i=1, m
334          t0=t0+dble(conjg(Oc(inw, i, i))*Oc(inw, i, i))
335       ENDDO
336    ENDDO
337
338    WRITE(6, *) t0
339
340
341        IF(abs(t0-oldt0)<tolw) THEN
342           WRITE(6,*) "MLWF Generated at Step",ini
343           GOTO 241
344        ENDIF
345
346        IF(adapt) THEN
347          IF(oldt0<t0) THEN
348            fric=fric/2.d0
349            A=Aminus
350            Aminus=temp
351          ENDIF
352        ENDIF
353
354!   calculate d(omega)/dA and store result in W
355!   this is the force for the damped dynamics
356!
357
358    W=0.d0
359    DO inw=1, nw
360       t2=weight(inw)
361       DO i=1,m
362          DO j=1,m
363             W(i,j)=W(i,j)+t2*dble(Oc(inw,i,j)*conjg(Oc(inw,i,i)        &
364                  -Oc(inw,j,j))+conjg(Oc(inw,j,i))*(Oc(inw,i,i)-Oc(inw,j,j)))
365          ENDDO
366       ENDDO
367    ENDDO
368
369
370!   the verlet scheme to calculate A(t+dt)
371
372        Aplus=0.d0
373
374   DO i=1,m
375     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     ENDDO
379   ENDDO
380
381        Aplus=Aplus-transpose(Aplus)
382        Aplus=(Aplus-A)
383
384    DO i=1, m
385       DO j=i,m
386        wp(i + (j-1)*j/2) = cmplx(0.0d0, Aplus(i,j), kind=8)
387       ENDDO
388    ENDDO
389
390    CALL zhpev('V','U',m,wp,wr,z,m,f1,f2,ierr1)
391!    call zhpev(21, wp, wr, z, m, m, f2, 4*m)
392
393    IF (ierr1/=0) THEN
394   WRITE(6,*) "failed to diagonalize W!"
395    STOP
396    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     CALL zgemm ('N','c', m,m,m,alpha,U3,m,z,m,beta1,U2,m)
409    U=dble(U2)
410    U2=beta1
411    U3=beta1
412
413   temp=Aminus
414   Aminus=A
415   A=Aplus
416
417
418!   update Umat
419!
420    U1=beta1
421    CALL dgemm ('N', 'N', m,m,m,alpha,Umat,m,U,m,beta1,U1,m)
422    Umat=U1
423
424!   update Oc
425!
426    U2=Umat*alpha
427    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    U2=beta1
436    U3=beta1
437
438        IF(abs(t0-oldt0)>=tolw.and.ini>=nsteps) THEN
439        GOTO 241
440        ENDIF
441
442    oldt0=t0
443
444   ENDDO
445241  spread=0.0d0
446  WRITE(24, *) "spread: (unit \AA^2)"
447  DO i=1, m
448     mt=1.d0-dble(Oc(:,i,i)*conjg(Oc(:,i,i)))
449     sp= (alat*autoaf/tpi)**2*sum(mt*weight)
450     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!   output wfc's and spreads of the max. loc. wf's
462!
463  ALLOCATE(wr(nw), W(nw, nw))
464  DO inw=1, nw
465     gr(inw, :)=wfg(inw,1)*b1(:)+wfg(inw,2)*b2(:)+wfg(inw,3)*b3(:)
466  ENDDO
467!
468! set up a matrix with the element (i,j) is G_i * G_j * weight(j)
469! to check the correctness of choices on G vectors
470!
471  DO i=1, nw
472     DO j=1, nw
473        W(i,j)=sum(gr(i,:)*gr(j,:))*weight(j)
474!        write(6, *) i,j,W(i,j)
475     ENDDO
476  ENDDO
477!  write(24, *) "wannier function centers: (unit:\AA)"
478  DO i=1, m
479     mt=-aimag(log(Oc(:,i,i)))/tpi
480     wfc(1, i)=sum(mt*weight*gr(:,1))
481     wfc(2, i)=sum(mt*weight*gr(:,2))
482     wfc(3, i)=sum(mt*weight*gr(:,3))
483     DO inw=1, nw
484        wr(inw)=sum(wfc(:,i)*gr(inw,:))-mt(inw)
485     ENDDO
486     mt=wr
487     f3=0
488     adjust=0
489!
490!   balance the phase factor if necessary
491!
492!     do while(SUM((mt-f3)**2).gt.0.01d0)
493!        f4=f3
494!        f3=nint(mt-mt(1))
495!        if (adjust.gt.200) f3=f3-1
496!        if (adjust.gt.100.and.adjust.le.200) f3=f3+1
497!        mt=wr+matmul(W, f3)
498!        write(6,*) "mt:", mt
499!        write(6,*) "f3:", f3
500!        adjust=adjust+1
501!        if (adjust.gt.300) stop "unable to balance the phase!"
502!     end do
503     wfc(1,i)=(wfc(1,i)+sum(mt*weight*gr(:,1)))*alat
504     wfc(2,i)=(wfc(2,i)+sum(mt*weight*gr(:,2)))*alat
505     wfc(3,i)=(wfc(3,i)+sum(mt*weight*gr(:,3)))*alat
506  ENDDO
507
508!  if (ibrav.eq.1.or.ibrav.eq.6.or.ibrav.eq.8) then
509!     do i=1, m
510!        if (wfc(1, i).lt.0) wfc(1, i)=wfc(1, i)+a1(1)
511!        if (wfc(2, i).lt.0) wfc(2, i)=wfc(2, i)+a2(2)
512!        if (wfc(3, i).lt.0) wfc(3, i)=wfc(3, i)+a3(3)
513!     end do
514!  end if
515  DO i=1, m
516     WRITE(26, '(3f11.6)') wfc(:,i)*autoaf
517  ENDDO
518
519     WRITE(6,*) "Friction =", fric
520     WRITE(6,*) "Mass =", q
521
522
523  DEALLOCATE(wr, W)
524
525  RETURN
526
527  END SUBROUTINE ddyn
528
529!-----------------------------------------------------------------------
530  SUBROUTINE searchwf(m, Omat, Umat)
531!-----------------------------------------------------------------------
532!    (m,m) is the size of the matrix Ospin.
533!    Ospin is input overlap matrix.
534!    Uspin is the output unitary transformation.
535!             Rough guess for Uspin can be carried in.
536!
537  USE wanpar
538  USE constants, ONLY : tpi, autoaf => BOHR_RADIUS_ANGS
539!
540!
541!     conjugated gradient to search maximization
542!
543  IMPLICIT NONE
544!
545  INTEGER, INTENT(in) :: m
546  COMPLEX(kind=8), INTENT(in) :: Omat(nw, m, m)
547  real(kind=8), INTENT(inout) :: Umat(m,m)
548!
549  INTEGER :: i, j, k, l, ig, ierr, ti, tj, tk, inw, ir
550  real(kind=8) :: slope, slope2, t1, t2, t3, mt(nw),t21,temp1,maxdt
551  real(kind=8) :: U(m,m), wfc(3, m), Wm(m,m), schd(m,m), f2(3*m-2), gr(nw, 3)
552  real(kind=8) :: Uspin2(m,m),temp2,wfdtold,oldt1,t01, d3(m,m), d4(m,m), U1(m,m)
553  real(kind=8), ALLOCATABLE, DIMENSION(:) :: wr
554  real(kind=8), ALLOCATABLE, DIMENSION(:,:) :: W
555  COMPLEX(kind=8) :: ci, ct1, ct2, ct3, z(m, m), X(m, m), d(m,m), d2(m,m)
556  COMPLEX(kind=8) :: f1(2*m-1), wp(m*(m+1)/2), Oc(nw, m, m), alpha, beta1
557  COMPLEX(kind=8) ::  Oc2(nw, m, m),wp1(m*(m+1)/2), X1(m,m), U2(m,m), U3(m,m)
558
559!
560  ci=(0.d0,1.d0)
561  alpha=(1.0d0, 0.0d0)
562  beta1=(0.0d0, 0.0d0)
563!
564  ALLOCATE(W(m,m), wr(m))
565
566
567!  Umat=0.d0
568!  do i=1,m
569!    Umat(i,i)=1.d0
570!  end do
571  Oc=beta1
572  Oc2=beta1
573  X1=beta1
574  U2=Umat*alpha
575
576!
577! update Oc using the initial guess of Uspin
578!
579  DO inw=1, nw
580     X1(:, :)=Omat(inw, :, :)
581     U3=beta1
582     CALL zgemm ('T', 'N', m,m,m,alpha,U2,m,X1,m,beta1,U3,m)
583     X1=beta1
584     CALL zgemm ('N','N', m,m,m,alpha,U3,m,U2,m,beta1,X1,m)
585     Oc(inw, :, :)=X1(:, :)
586  ENDDO
587
588     U2=beta1
589     U3=beta1
590
591  W=0.d0
592  schd=0.d0
593  oldt1=0.d0
594  wfdtold=0.d0
595
596  DO k=1, nit
597    t01=0.d0     !use t1 to store the value of omiga
598    DO inw=1, nw
599       DO i=1, m
600          t01=t01+dble(conjg(Oc(inw, i, i))*Oc(inw, i, i))
601       ENDDO
602    ENDDO
603
604    WRITE(6,*) t01
605
606    IF(abs(oldt1-t01)<tolw) GOTO 40
607
608    oldt1=t01
609
610!   calculate d(omiga)/dW and store result in W
611!   W should be a real symmetric matrix for gamma-point calculation
612!
613    Wm=W
614    W=0.d0
615    DO inw=1, nw
616       t2=weight(inw)
617       DO i=1,m
618          DO j=i+1,m
619             W(i,j)=W(i,j)+t2*dble(Oc(inw,i,j)*conjg(Oc(inw,i,i)        &
620              -Oc(inw,j,j))+conjg(Oc(inw,j,i))*(Oc(inw,i,i)-Oc(inw,j,j)))
621          ENDDO
622       ENDDO
623    ENDDO
624    W=W-transpose(W)
625
626!   calculate slope=d(omiga)/d(lamda)
627    slope=sum(W**2)
628
629!   calculate slope2=d2(omiga)/d(lamda)2
630    slope2=0.d0
631    DO ti=1, m
632       DO tj=1, m
633          DO tk=1, m
634             t2=0.d0
635             DO inw=1, nw
636                t2=t2+dble(Oc(inw,tj,tk)*conjg(Oc(inw,tj,tj)+Oc(inw,tk,tk) &
637                          -2.d0*Oc(inw,ti,ti))-4.d0*Oc(inw,ti,tk)          &
638                          *conjg(Oc(inw,ti,tj)))*weight(inw)
639             ENDDO
640             slope2=slope2+W(tk,ti)*W(ti,tj)*2.d0*t2
641          ENDDO
642       ENDDO
643     ENDDO
644    slope2=2.d0*slope2
645
646!   use parabola approximation. Defined by 1 point and 2 slopes
647    IF (slope2<0) wfdt=-slope/2.d0/slope2
648    IF (maxwfdt>0.and.wfdt>maxwfdt) wfdt=maxwfdt
649
650    IF (k<nsd) THEN
651       schd=W    !use steepest-descent technique
652
653!   calculate slope=d(omiga)/d(lamda)
654    slope=sum(schd**2)
655
656!       schd=schd*maxwfdt
657    DO i=1, m
658       DO j=i, m
659        wp1(i + (j-1)*j/2) = cmplx(0.0d0, schd(i,j), kind=8)
660       ENDDO
661    ENDDO
662
663    CALL zhpev('V','U',m,wp1,wr,z,m,f1,f2,ierr)
664    IF (ierr/=0) STOP 'failed to diagonalize W!'
665
666    ELSE
667!
668!     construct conjugated gradient
669!        d(i)=-g(i)+belta(i)*d(i-1)
670!        belta^FR(i)=g(i)t*g(i)/(g(i-1)t*g(i-1))
671!        belta^PR(i)=g(i)t*(g(i)-g(i-1))/(g(i-1)t*g(i-1))
672!
673        CALL dgemm ('T','N', m,m,m,alpha,Wm,m,Wm,m,beta1,d3,m)
674
675
676       t1=0.d0
677       DO i=1, m
678          t1=t1+d3(i, i)
679       ENDDO
680       IF (t1/=0) THEN
681          d4=(W-Wm)
682          CALL dgemm ('T','N', m,m,m,alpha,W,m,d4,m,beta1,d3,m)
683          t2=0.d0
684          DO i=1, m
685             t2=t2+d3(i, i)
686          ENDDO
687          t3=t2/t1
688          schd=W+schd*t3
689       ELSE
690          schd=W
691       ENDIF
692!
693!        calculate the new d(Lambda) for the new Search Direction
694!        added by Manu. September 19, 2001
695!
696!   calculate slope=d(omiga)/d(lamda)
697    slope=sum(schd**2)
698!------------------------------------------------------------------------
699!   schd=schd*maxwfdt
700    DO i=1, m
701       DO j=i, m
702        wp1(i + (j-1)*j/2) = cmplx(0.0d0, schd(i,j), kind=8)
703       ENDDO
704    ENDDO
705
706    CALL zhpev('V','U',m,wp1,wr,z,m,f1,f2,ierr)
707    IF (ierr/=0) STOP 'failed to diagonalize W!'
708
709      maxdt=maxwfdt
710
71111    d=0.d0
712    DO i=1, m
713       d(i, i)=exp(ci*(maxwfdt)*wr(i))
714    ENDDO      !d=exp(d)
715
716!   U=z*exp(d)*z+
717     U3=beta1
718     CALL zgemm ('N', 'N', m,m,m,alpha,z,m,d,m,beta1,U3,m)
719     U2=beta1
720     CALL zgemm ('N','c', m,m,m,alpha,U3,m,z,m,beta1,U2,m)
721     U=dble(U2)
722     U2=beta1
723     U3=beta1
724!
725!   update Uspin
726    U1=beta1
727    CALL dgemm ('N', 'N', m,m,m,alpha,Umat,m,U,m,beta1,U1,m)
728    Umat=U1
729
730!    Uspin2=matmul(Uspin, U2)
731!
732!   update Oc
733!
734     U2=Umat*alpha
735     U3=beta1
736     DO inw=1, nw
737      X1(:,:)=Omat(inw,:,:)
738      CALL zgemm ('T', 'N', m,m,m,alpha,U2,m,X1,m,beta1,U3,m)
739      X1=beta1
740      CALL zgemm ('N','N',m,m,m,alpha,U3,m,U2,m,beta1,X1,m)
741      Oc2(inw, :,:)=X(:,:)
742     ENDDO
743     U2=beta1
744     U3=beta1
745!
746    t21=0.d0     !use t21 to store the value of omiga
747    DO inw=1, nw
748       DO i=1, m
749          t21=t21+dble(conjg(Oc2(inw, i, i))*Oc2(inw, i, i))
750       ENDDO
751    ENDDO
752
753      temp1=-((t01-t21)+slope*maxwfdt)/(maxwfdt**2)
754      temp2=slope
755      wfdt=-temp2/(2*temp1)
756
757        IF (wfdt>maxwfdt.or.wfdt<0.d0) THEN
758        maxwfdt=2*maxwfdt
759        GOTO 11
760        ENDIF
761
762        maxwfdt=maxdt
763!
764!
765!   use parabola approximation. Defined by 2 point and 1 slopes
766!    if (slope2.lt.0) wfdt=-slope/2.d0/slope2
767!    if (maxwfdt.gt.0.and.wfdt.gt.maxwfdt) wfdt=maxwfdt
768!
769!    write(6, '(e12.5E2,1x,e11.5E2,1x,f6.2)') slope2, slope, wfdt
770!-------------------------------------------------------------------------
771!
772!      schd is the new searching direction
773!
774    ENDIF
775
776    d=0.d0
777    DO i=1, m
778       d(i, i)=exp(ci*wfdt*wr(i))
779    ENDDO          !d=exp(d)
780
781
782!   U=z*exp(d)*z+
783!
784     U3=beta1
785     CALL zgemm ('N', 'N', m,m,m,alpha,z,m,d,m,beta1,U3,m)
786     U2=beta1
787     CALL zgemm ('N','c', m,m,m,alpha,U3,m,z,m,beta1,U2,m)
788     U=dble(U2)
789     U2=beta1
790     U3=beta1
791
792!   update Uspin
793!
794    U1=beta1
795    CALL dgemm ('N', 'N', m,m,m,alpha,Umat,m,U,m,beta1,U1,m)
796    Umat=U1
797
798!   update Oc
799!
800       U2=Umat*alpha
801       U3=beta1
802     DO inw=1, nw
803       X1(:, :)=Omat(inw, :, :)
804       CALL zgemm ('T', 'N', m,m,m,alpha,U2,m,X1,m,beta1,U3,m)
805       X1=beta1
806       CALL zgemm ('N','N',m,m,m,alpha,U3,m,U2,m,beta1,X1,m)
807       Oc(inw, :, :)=X1(:, :)
808     ENDDO
809    U2=beta1
810    U3=beta1
811  ENDDO
812
81340  DEALLOCATE(W, wr)
814
815!
816! calculate the spread
817!
818  WRITE(24, *) "spread: (unit \AA^2)"
819  DO i=1, m
820     mt=1.d0-dble(Oc(:,i,i)*conjg(Oc(:,i,i)))
821     WRITE(24, '(f10.7)') (alat*autoaf/tpi)**2*sum(mt*weight)
822  ENDDO
823
824!
825! calculate wannier-function centers
826!
827  ALLOCATE(wr(nw), W(nw, nw))
828  DO inw=1, nw
829     gr(inw, :)=wfg(inw,1)*b1(:)+wfg(inw,2)*b2(:)+wfg(inw,3)*b3(:)
830  ENDDO
831!
832! set up a matrix with the element (i,j) is G_i * G_j * weight(j)
833! to check the correctness of choices on G vectors
834!
835  DO i=1, nw
836     DO j=1, nw
837        W(i,j)=sum(gr(i,:)*gr(j,:))*weight(j)
838!        write(6, *) i,j,W(i,j)
839     ENDDO
840  ENDDO
841! write(24, *) "wannier function centers: (unit:\AA)"
842  DO i=1, m
843     mt=-aimag(log(Oc(:,i,i)))/tpi
844     wfc(1, i)=sum(mt*weight*gr(:,1))
845     wfc(2, i)=sum(mt*weight*gr(:,2))
846     wfc(3, i)=sum(mt*weight*gr(:,3))
847     DO inw=1, nw
848        wr(inw)=sum(wfc(:,i)*gr(inw,:))-mt(inw)
849     ENDDO
850     mt=wr
851     wfc(1, i)=(wfc(1,i)+sum(mt*weight*gr(:,1)))*alat
852     wfc(2, i)=(wfc(2,i)+sum(mt*weight*gr(:,2)))*alat
853     wfc(3, i)=(wfc(3,i)+sum(mt*weight*gr(:,3)))*alat
854  ENDDO
855
856
857  DO i=1, m
858     WRITE(26, '(3f11.6)') wfc(:,i)*autoaf
859  ENDDO
860  DEALLOCATE(wr, W)
861  RETURN
862  END SUBROUTINE searchwf
863
864
865 END PROGRAM wfdd
866