1* 2* $Id$ 3* 4 5 6* *********************************** 7* * * 8* * paw_psi_MakeOrtho * 9* * * 10* *********************************** 11 12* This routine orthonormalizes the orbitals using a modified 13* Gram-Schmidt algorithm. 14* 15 subroutine paw_psi_MakeOrtho(npack,ne,psi) 16 implicit none 17 integer npack,ne 18 double complex psi(npack,ne) 19 20* **** local variables **** 21 integer j,k 22 real*8 w 23 24 do k=1,ne 25 call paw_overlap_matrix_gen(1,1,psi(1,k),psi(1,k),w) 26 w = 1.0d0/dsqrt(w) 27c call Pack_c_SMul(1,w,psi(1,k),psi(1,k)) 28 call Pack_c_SMul1(1,w,psi(1,k)) 29 30 do j=k+1,ne 31 call paw_overlap_matrix_gen(1,1,psi(1,k),psi(1,j),w) 32 w = -w 33 call Pack_cc_daxpy(1,w,psi(1,k),psi(1,j)) 34 end do 35 end do 36 37 return 38 end 39 40* *********************************** 41* * * 42* * paw_psi_CheckOrtho * 43* * * 44* *********************************** 45 46* This routine return true if the orbitals are 47* orthonormal. 48 49 real*8 function paw_psi_CheckOrtho(npack,ne,psi) 50 implicit none 51 integer npack,ne 52 double complex psi(npack,ne) 53 54 55* **** local variables **** 56 integer j,k 57 real*8 w,error 58 59 error = 0.0d0 60 do k=1,ne 61 call paw_overlap_matrix_gen(1,1,psi(1,k),psi(1,k),w) 62 error = error + dabs(1.0d0-w) 63 64 do j=k+1,ne 65 call paw_overlap_matrix_gen(1,1,psi(1,j),psi(1,k),w) 66 error = error + dabs(w) 67 end do 68 end do 69 70 paw_psi_CheckOrtho = error 71 return 72 end 73 74 75* *********************************** 76* * * 77* * paw_psi_CheckOrtho2 * 78* * * 79* *********************************** 80 81* This routine return true if the orbitals are 82* orthonormal. 83 84 subroutine paw_psi_CheckOrtho2(npack,ne,psi) 85 implicit none 86 integer npack,ne 87 double complex psi(npack,ne) 88 89 90* **** local variables **** 91 integer j,k 92 real*8 w,error 93 94 95 write(*,*) 96 do k=1,ne 97 call paw_overlap_matrix_gen(1,1,psi(1,k),psi(1,k),w) 98 write(*,*) "CheckOrtho2:",k,k,w 99 100 do j=k+1,ne 101 call paw_overlap_matrix_gen(1,1,psi(1,j),psi(1,k),w) 102 write(*,*) "CheckOrtho2:",j,k,w 103 end do 104 end do 105 write(*,*) 106 107 return 108 end 109 110 111 112* *********************************** 113* * * 114* * paw_psi_lmbda * 115* * * 116* *********************************** 117 118 subroutine paw_psi_lmbda(ispin,ne,nemax,npack1, 119 > psi1,psi2, 120 > dte, 121 > lmbda,tmp,ierr) 122 123 implicit none 124 integer ispin,ne(2),nemax,npack1 125 complex*16 psi1(npack1,nemax) 126 complex*16 psi2(npack1,nemax) 127 real*8 dte 128 real*8 lmbda(*) 129 real*8 tmp(*) 130 integer ierr 131 132 integer MASTER 133 parameter (MASTER=0) 134 135* **** local variables **** 136 logical failed 137 integer taskid 138 integer n1(2),n2(2) 139 integer i,j,ii,jj,ms 140 integer n,nn,index 141 integer st1,st2 142 integer A,B,C,U,D,Ba,Bs,fnm 143 integer sl(2) 144 real*8 alpha 145 146 147 148 call nwpw_timing_start(3) 149 150 call Parallel_taskid(taskid) 151 152 n = ne(1) 153 nn = n**2 154 155 A = 0*nn + 1 156 B = 1*nn + 1 157 C = 2*nn + 1 158 Ba = 3*nn + 1 159 Bs = 4*nn + 1 160 fnm = 5*nn + 1 161 st1 = 6*nn + 1 162 D = 7*nn + 1 163 164 U = Bs 165 st2 = B 166 167 call dcopy(8*nn,0.0d0,0,tmp,1) 168 169 sl(1) = 0*nn + 1 170 sl(2) = 1*nn + 1 171 call dcopy(2*nn,0.0d0,0,lmbda,1) 172 173 n1(1)=1 174 n2(1)=ne(1) 175 n1(2)=ne(1)+1 176 n2(2)=ne(1)+ne(2) 177 178 179 do ms=1,ispin 180 IF(ne(ms).le.0) go to 640 181 182 183* ***** compute the overlap matrices **** 184 call paw_overlap_sym_matrix_gen(n,ne(ms), 185 > psi2(1,n1(ms)), 186 > psi2(1,n1(ms)), 187 > tmp(A)) 188 call paw_overlap_matrix_gen(n,ne(ms), 189 > psi1(1,n1(ms)), 190 > psi2(1,n1(ms)), 191 > tmp(B)) 192 call paw_overlap_sym_matrix_gen(n,ne(ms), 193 > psi1(1,n1(ms)), 194 > psi1(1,n1(ms)), 195 > tmp(C)) 196 197 198 call paw_psi_gen_Ba_Bs(n,ne(ms),tmp(B),tmp(Bs),tmp(Ba)) 199 call paw_psi_gen_UD(n,ne(ms),tmp(Bs),tmp(D),lmbda) 200 201 202 call paw_psi_gen_X(n,ne(ms),tmp(st1),tmp(st2), 203 > tmp(A),tmp(Ba),tmp(C), 204 > tmp(U),tmp(D),tmp(fnm), 205 > failed) 206 207 if (failed) then 208 if (taskid.eq.MASTER) then 209 write(6,*) 210 > 'Warning: Lagrange Multiplier generation failed.' 211 write(6,*) ' +Try using a smaller time step' 212 write(6,*) ' +Gram-Schmidt being performed, spin:',ms 213 end if 214 call paw_psi_MakeOrtho(npack1,ne(ms),psi2(1,n1(ms))) 215 else 216 call dcopy(n*ne(ms),tmp(st1),1,lmbda(sl(ms)),1) 217 call dscal(n*ne(ms),(1.0d0/dte),lmbda(sl(ms)),1) 218 219 220* **** correction due to the constraint **** 221 call dgemm('N','N',2*npack1,ne(ms),ne(ms), 222 > (1.0d0), 223 > psi1(1,n1(ms)),2*npack1, 224 > tmp(st1),n, 225 > (1.0d0), 226 > psi2(1,n1(ms)),2*npack1) 227 228 end if 229 640 continue 230 end do !*ms* 231 call nwpw_timing_end(3) 232 233 return 234 end 235 236 237 238* *********************************** 239* * * 240* * paw_psi_gen_Ba_Bs * 241* * * 242* *********************************** 243 subroutine paw_psi_gen_Ba_Bs(n_max,n,B,Bs,Ba) 244 implicit none 245 integer n_max,n 246 real*8 B(n_max,n) 247 real*8 Bs(n_max,n) 248 real*8 Ba(n_max,n) 249 250 !*** local variables *** 251 integer i,j 252 253 do i=1,n 254 do j=1,n 255 Bs(i,j) = 0.5d0*(B(i,j)+B(j,i)) 256 Ba(i,j) = 0.5d0*(B(i,j)-B(j,i)) 257 end do 258 end do 259 return 260 end 261 262* *********************************** 263* * * 264* * paw_psi_gen_UD * 265* * * 266* *********************************** 267 subroutine paw_psi_gen_UD(n_max,n,Bs,D,work) 268 implicit none 269 integer n_max,n 270 real*8 Bs(n_max,n) 271 real*8 D(n_max,n) 272 real*8 Work(n_max,n) 273 274 !*** local variables *** 275 integer ierr 276 277 !call eigen(n_max,n,Bs,D,D(1,2)) 278 call dsyev('V','U',n,Bs,n_max, D,Work,2*n_max*n_max,ierr) 279 return 280 end 281 282 283 284 285* *********************************** 286* * * 287* * paw_psi_gen_X * 288* * * 289* *********************************** 290 subroutine paw_psi_gen_X(n_max,n, 291 > X1,tmp, 292 > A,Ba,C, 293 > U,D,fnm, 294 > failed) 295 > 296 implicit none 297 integer n_max,n 298 real*8 X1(n_max,n) 299 real*8 tmp(*) 300 real*8 A(n_max,n) 301 real*8 Ba(n_max,n) 302 real*8 C(n_max,n) 303 real*8 U(n_max,n) 304 real*8 D(n_max,n) 305 real*8 fnm(n_max,n) 306 logical failed 307 308 !**** local variables **** 309 integer itrlmd 310 real*8 convg 311 parameter (itrlmd=40, convg=1.0d-15) 312 313 integer i,it 314 real*8 adiff 315 316 !**** external functions **** 317 integer idamax 318 external idamax 319 320 321 !**** A = I-A *** 322 call dscal(n_max*n,(-1.0d0),A,1) 323 do i=1,n 324 A(i,i) = A(i,i) + 1.0d0 325 end do 326 327 !*** fnm = I-A **** 328 call dcopy(n_max*n,A,1,fnm,1) 329 330 !*** solve U*D*Ut*X + X*U*D*Ut = fnm for X *** 331 call paw_psi_fnm_to_X(n_max,n,fnm,U,D,tmp) 332 call dcopy(n_max*n,fnm,1,X1,1) 333 334 335 it = 0 336 failed = .true. 337 do while (failed .and. (it.lt.itrlmd)) 338 it = it + 1 339 340 !*** fnm = X*C*X *** 341 call DMMUL(n_max,n,C,X1,tmp) 342 call DMMUL(n_max,n,X1,tmp,fnm) 343 344 345 !*** fnm = Ba*X - X*C*X *** 346 call DMMUL(n_max,n,Ba,X1,tmp) 347 call DMSUB(n_max,n,tmp,fnm,fnm) 348 349 350 !*** fnm = Ba*X - X*Ba - X*C*X *** 351 call DMMUL(n_max,n,X1,Ba,tmp) 352 call DMSUB(n_max,n,fnm,tmp,fnm) 353 354 !*** fnm = I-A + Ba*X - X*Ba - X*C*X *** 355 call DMADD(n_max,n,fnm,A,fnm) 356 357 358 !*** solve U*D*Ut*X + X*U*D*Ut = fnm for X *** 359 call paw_psi_fnm_to_X(n_max,n,fnm,U,D,tmp) 360 361 call DMSUB(n_max,n,fnm,X1,tmp) 362 adiff = tmp(idamax(n_max*n,tmp,1)) 363 call dcopy(n_max*n,fnm,1,X1,1) 364 365 adiff = dabs(adiff) 366 367 if (adiff.lt.convg) failed = .false. 368 end do 369 370 return 371 end 372 373 374* *********************************** 375* * * 376* * paw_psi_fnm_to_X * 377* * * 378* *********************************** 379 subroutine paw_psi_fnm_to_X(n_max,n, 380 > fnm,U,D,tmp) 381 implicit none 382 integer n_max,n 383 real*8 fnm(n_max,n) 384 real*8 U(n_max,n) 385 real*8 D(n_max,n) 386 real*8 tmp(n_max,n) 387 388 !**** local variables **** 389 integer i,j 390 real*8 d2 391 392 393 !**** fnm = Ut*fnm*U *** 394 call dgemm('N','N',n,n,n,1.0d0, 395 > fnm,n_max, 396 > U,n_max, 397 > 0.0d0, 398 > tmp,n_max) 399 call dgemm('T','N',n,n,n,1.0d0, 400 > U,n_max, 401 > tmp,n_max, 402 > 0.0d0, 403 > fnm,n_max) 404 405 406 !**** fnm = (Ut*fnm*U)_nm/(d_n+d_m) *** 407 do j=1,n 408 do i=1,n 409 d2 = D(i,1)+D(j,1) 410 fnm(i,j) = fnm(i,j)/d2 411 end do 412 end do 413 414 !**** fnm = X = U*{(Ut*fnm*U)_nm/(d_n+d_m)}*Ut *** 415 call dgemm('N','N',n,n,n,1.0d0, 416 > U,n_max, 417 > fnm,n_max, 418 > 0.0d0, 419 > tmp,n_max) 420 call dgemm('N','T',n,n,n,1.0d0, 421 > tmp,n_max, 422 > U,n_max, 423 > 0.0d0, 424 > fnm,n_max) 425 426 return 427 end 428