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