1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2! This module contains the following active subroutines: 3! besselps, EvaluateP, EvaluateTp, genOrthog, hatpotL, hatL, 4! dqij, dtij, altdtij, dvij, avij, calcwij, FillHat, 5! Get_Energy_EXX_pseudo, Calc_Moment, Get_Energy_EXX_pseudo_one, 6! SPMatrixElements, COREVAL_EXX, CORECORE_EXX 7! This module contains the following active functions: 8! sepnorm, genoverlap 9!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 10 11#if defined HAVE_CONFIG_H 12#include "config.h" 13#endif 14 15MODULE pseudo_sub 16 17 USE io_tools 18 USE globalmath 19 USE gridmod 20 USE pseudodata 21 USE excor 22 USE fock 23 24 IMPLICIT NONE 25 26CONTAINS 27 28 !*************************************************************** 29 ! SUBROUTINE besselps(lmax,Grid,Pot) 30 ! Creates screened pseudopotential by simply pseudizing the 31 ! AE potential with a l=0 spherical Bessel function: 32 ! Vps(r) = a.sin(qr)/r 33 !*************************************************************** 34 SUBROUTINE besselps(Grid,Pot,PAW) 35 TYPE(Gridinfo), INTENT(IN) :: Grid 36 TYPE(Potentialinfo), INTENT(IN) :: Pot 37 TYPE(Pseudoinfo), INTENT(INOUT) :: PAW 38 39 INTEGER :: i,irc,l,n 40 REAL(8) :: e,rc,alpha,beta,vv,vvp,AA,QQ,xx(1) 41 REAL(8),ALLOCATABLE :: VNC(:),wfn(:) 42 REAL(8),POINTER :: r(:),rv(:) 43 CHARACTER(132) :: line 44 45 n=Grid%n 46 r=>Grid%r 47 rv=>Pot%rv 48 irc=PAW%irc_vloc 49 rc=PAW%rc_vloc 50 51 vv=rv(irc);vvp=Gfirstderiv(Grid,irc,rv) 52 53 alpha=1.D0-rc*vvp/vv;beta=1.D0 54 call solvbes(xx,alpha,beta,0,1);QQ=xx(1) 55 AA=vv/sin(QQ);QQ=QQ/rc 56 57 PAW%rveff(1)=0.d0 58 PAW%rveff(irc+1:n)=rv(irc+1:n) 59 do i=2,irc 60 PAW%rveff(i)=AA*sin(QQ*r(i)) 61 enddo 62 63 END SUBROUTINE besselps 64 65 66 !*************************************************************** 67 ! SUBROUTINE EvaluateP 68 ! Inverts 4x4 matrix used by kerker subroutine 69 !*************************************************************** 70 SUBROUTINE EvaluateP(m,A,B,C,D,coef) 71 INTEGER, INTENT(IN) :: m(4) 72 REAL(8), INTENT(IN) :: A,B,C,D 73 REAL(8), INTENT(OUT) :: coef(4) 74 75 REAL(8) :: t(4,4) 76 INTEGER :: i,n 77 78 t=0 79 Coef(1)=A; Coef(2)=B; Coef(3)=C; Coef(4)=D 80 t(1,1:4)=1 81 t(2,1:4)=m(1:4) 82 DO i=1,4 83 t(3,i)=m(i)*(m(i)-1) 84 ENDDO 85 DO i=1,4 86 t(4,i)=m(i)*(m(i)-1)*(m(i)-2) 87 ENDDO 88 n=4 89 CALL linsol(t,Coef,n,4,4,4) 90 END SUBROUTINE EvaluateP 91 92 !*************************************************************** 93 ! SUBROUTINE EvaluateTp 94 ! Inverts 5x5 matrix used by troullier subroutine 95 !*************************************************************** 96 SUBROUTINE EvaluateTp(l,A,B,C,D,F,coef) 97 INTEGER, INTENT(IN) :: l 98 REAL(8), INTENT(IN) :: A,B,C,D,F 99 REAL(8), INTENT(OUT) :: coef(6) 100 101 REAL(8) :: t(6,6),coef10,old 102 REAL(8), PARAMETER :: small=1.e-10 103 INTEGER :: i,n,iter 104 INTEGER, PARAMETER :: niter=1000 105 106 old=-1.e30; Coef10=-1; iter=-1 107 DO WHILE (iter < niter .AND. ABS(old-coef10)> small) 108 iter=iter+1 109 t=0 110 Coef(1)=A-Coef10; Coef(2)=B-2*Coef10; Coef(3)=C-2*Coef10; 111 Coef(4)=D; Coef(5)=F 112 Coef(6)=-Coef10**2 113 DO i=1,6 114 t(1,i)=1 115 t(2,i)=2*i 116 t(3,i)=2*i*(2*i-1) 117 t(4,i)=2*i*(2*i-1)*(2*i-2) 118 t(5,i)=2*i*(2*i-1)*(2*i-2)*(2*i-3) 119 ENDDO 120 t(6,1)=2*Coef10; t(6,2)=2*l+5 121 122 n=6 123 CALL linsol(t,Coef,n,6,6,6) 124 125 old=Coef10; Coef10=Coef10+Coef(1) 126 !WRITE(STD_OUT,'("EvaluateTp: iter",i5,1p,2e15.7)') iter,Coef(1),Coef10 127 !WRITE(STD_OUT,'("Coef: ",1p,6e15.7)')Coef10,(Coef(i),i=2,6) 128 Coef(1)=Coef10 129 ENDDO 130 131 IF (iter >= niter) THEN 132 WRITE(STD_OUT,*) 'Error in EvaluateTP -- no convergence' 133 STOP 134 ENDIF 135 END SUBROUTINE EvaluateTp 136 137 !******************************************************************* 138 ! function to calculated <wfn|O|wfn> for smooth paw wavefunction 139 !******************************************************************* 140 FUNCTION sepnorm(Grid,PAW,nr,l,wfn) 141 REAL(8) :: sepnorm 142 TYPE(GridInfo), INTENT(IN) :: Grid 143 TYPE(PseudoInfo), INTENT(IN) :: PAW 144 INTEGER, INTENT(IN) :: nr,l 145 REAL(8), INTENT(IN) :: wfn(:) 146 147 INTEGER :: n,ib,ic,nbase,irc 148 REAL(8) :: h 149 REAL(8), ALLOCATABLE :: b(:) 150 151 n=Grid%n; h=Grid%h; nbase=PAW%nbase; irc=PAW%irc 152 ALLOCATE(b(nbase),stat=ib) 153 IF (ib/=0) THEN 154 WRITE(STD_OUT,*) 'Error in sepnorm allocation', nbase,ib 155 STOP 156 ENDIF 157 158 IF (nr<irc) THEN 159 WRITE(STD_OUT,*) 'Error in sepnorm -- nr < irc' 160 STOP 161 ENDIF 162 sepnorm=overlap(Grid,wfn,wfn,1,nr) 163 b=0 164 DO ib=1,nbase 165 IF (l==PAW%l(ib)) b(ib)=overlap(Grid,wfn,PAW%otp(:,ib),1,nr) 166 ENDDO 167 DO ib=1,nbase 168 DO ic=1,nbase 169 sepnorm=sepnorm+b(ib)*PAW%oij(ib,ic)*b(ic) 170 ENDDO 171 ENDDO 172 173 DEALLOCATE(b) 174 END FUNCTION sepnorm 175 176 !******************************************************************* 177 ! function to calculated <wfn1|O|wfn2> for smooth paw functions 178 !******************************************************************* 179 FUNCTION genoverlap(Grid,PAW,nr,l,wfn1,wfn2) 180 REAL(8) :: genoverlap 181 TYPE(GridInfo), INTENT(IN) :: Grid 182 TYPE(PseudoInfo), INTENT(IN) :: PAW 183 INTEGER, INTENT(IN) :: nr,l 184 REAL(8), INTENT(IN) :: wfn1(:),wfn2(:) 185 186 INTEGER :: n,ib,ic,nbase,irc,p 187 REAL(8) :: h 188 REAL(8), ALLOCATABLE :: a(:),b(:) 189 190 n=Grid%n; h=Grid%h; nbase=PAW%nbase; irc=PAW%irc;p=irc+1 191 ALLOCATE(a(nbase),b(nbase),stat=ib) 192 IF (ib/=0) THEN 193 WRITE(STD_OUT,*) 'Error in genoverlap allocation', nbase,ib 194 STOP 195 ENDIF 196 197 IF (nr<irc) THEN 198 WRITE(STD_OUT,*) 'Error in genoverlap -- nr < irc' 199 STOP 200 ENDIF 201 genoverlap=overlap(Grid,wfn1,wfn2,1,nr) 202 a=0; b=0 203 DO ib=1,nbase 204 IF (l==PAW%l(ib)) THEN 205 a(ib)=overlap(Grid,wfn1,PAW%otp(:,ib)) 206 b(ib)=overlap(Grid,wfn2,PAW%otp(:,ib)) 207 ENDIF 208 ENDDO 209 DO ib=1,nbase 210 DO ic=1,nbase 211 genoverlap=genoverlap+a(ib)*PAW%oij(ib,ic)*b(ic) 212 ENDDO 213 ENDDO 214 215 DEALLOCATE(a,b) 216 END FUNCTION genoverlap 217 218 219 ! generalized Graham-Schmidt orthogonalization of wfn1 to wfn2 220 ! on output <wfn1|O|wfn2>=0 221 SUBROUTINE genOrthog(Grid,PAW,nr,l,wfn1,wfn2) 222 !REAL(8) :: genoverlap 223 TYPE(GridInfo), INTENT(IN) :: Grid 224 TYPE(PseudoInfo), INTENT(IN) :: PAW 225 INTEGER, INTENT(IN) :: nr,l 226 REAL(8), INTENT(INOUT) :: wfn1(:) 227 REAL(8), INTENT(IN) :: wfn2(:) 228 229 REAL(8) :: x,y 230 231 ! Orthogonalize 232 x=genoverlap(Grid,PAW,nr,l,wfn1,wfn2) 233 y=genoverlap(Grid,PAW,nr,l,wfn2,wfn2) 234 235 WRITE(STD_OUT,*) 'overlap ', x,y 236 wfn1(:)=wfn1(:)-(x/y)*wfn2(:) 237 238 END SUBROUTINE genOrthog 239 240 !************************************************************** 241 ! subroutine hatpotL 242 ! Calculates potential associated with L component 243 ! of unit hat density 244 !************************************************************** 245 SUBROUTINE hatpotL(Grid,PAW,l,vhat) 246 TYPE(GridInfo), INTENT(IN) :: Grid 247 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 248 INTEGER, INTENT(IN) :: l 249 REAL(8), INTENT(OUT) :: vhat(:) 250 251 INTEGER :: n,irc,i 252 REAL(8), POINTER :: r(:) 253 REAL(8), ALLOCATABLE :: den(:),a(:) 254 REAL(8) :: h,con 255 REAL(8) :: qr,jbes1,jbes2,dum1,dum2,al(2),ql(2) 256 257 n=Grid%n 258 h=Grid%h 259 r=>Grid%r 260 261 irc=PAW%irc 262 263 ALLOCATE(den(n),a(n),stat=i) 264 IF (i/=0) THEN 265 WRITE(STD_OUT,*) 'Error in hatpotL allocation',n,i 266 STOP 267 ENDIF 268 269 270 IF (besselshapefunction) THEN 271 CALL shapebes(al,ql,l,PAW%rc_shap) 272 DO i=1,PAW%irc_shap 273 qr=ql(1)*r(i);CALL jbessel(jbes1,dum1,dum2,l,0,qr) 274 qr=ql(2)*r(i);CALL jbessel(jbes2,dum1,dum2,l,0,qr) 275 den(i)=(al(1)*jbes1+al(2)*jbes2)*r(i)**2 276 ENDDO 277 IF (n>PAW%irc_shap) den(PAW%irc_shap+1:n)=0.d0 278 ELSE 279 DO i=1,n 280 den(i)=(r(i)**l)*PAW%hatden(i) 281 ENDDO 282 a(1:n)=den(1:n)*(r(1:n)**l) 283 con=integrator(Grid,a,1,PAW%irc_shap) 284 den=den/con 285 ENDIF 286 287 a(1:n)=den(1:n)*(r(1:n)**l) 288 write(std_out,*) 'hatpotL check l ', l,integrator(Grid,a) 289 290 vhat=0 291 CALL apoisson(Grid,l,n,den,vhat(1:n)) 292 293 ! apoisson returns vhat*r 294 !DO i=1,n 295 ! WRITE (78+l,'(i5,1p,5e15.7)') i,Grid%r(i),den(i),vhat(i) 296 !ENDDO 297 vhat(2:n)=vhat(2:n)/r(2:n) 298 CALL extrapolate(Grid,vhat) 299 300 DEALLOCATE(den,a) 301 END SUBROUTINE hatpotL 302 303 !************************************************************** 304 ! subroutine hatL 305 ! Calculates density associated with L component 306 ! normalized to unity 307 !************************************************************** 308 SUBROUTINE hatL(Grid,PAW,l,dhat) 309 TYPE(GridInfo), INTENT(IN) :: Grid 310 TYPE(PseudoInfo), INTENT(IN) :: PAW 311 INTEGER, INTENT(IN) :: l 312 REAL(8), INTENT(OUT) :: dhat(:) 313 314 INTEGER :: n,irc,i 315 REAL(8), POINTER :: r(:) 316 REAL(8), ALLOCATABLE :: den(:),a(:) 317 REAL(8) :: h,con 318 REAL(8) :: qr,jbes1,jbes2,dum1,dum2,al(2),ql(2) 319 320 n=Grid%n 321 h=Grid%h 322 r=>Grid%r 323 324 irc=PAW%irc 325 326 ALLOCATE(den(n),a(n),stat=i) 327 IF (i/=0) THEN 328 WRITE(STD_OUT,*) 'Error in hatL allocation',irc,i 329 STOP 330 ENDIF 331 332 IF (besselshapefunction) THEN 333 CALL shapebes(al,ql,l,PAW%rc_shap) 334 DO i=1,PAW%irc_shap 335 qr=ql(1)*r(i);CALL jbessel(jbes1,dum1,dum2,l,0,qr) 336 qr=ql(2)*r(i);CALL jbessel(jbes2,dum1,dum2,l,0,qr) 337 den(i)=(al(1)*jbes1+al(2)*jbes2)*r(i)**2 338 ENDDO 339 IF (n>PAW%irc_shap) den(PAW%irc_shap+1:n)=0.d0 340 ELSE 341 DO i=1,n 342 den(i)=(r(i)**l)*PAW%hatden(i) 343 ENDDO 344 a(1:n)=den(1:n)*(r(1:n)**l) 345 con=integrator(Grid,a,1,PAW%irc_shap) 346 den=den/con 347 ENDIF 348 349 dhat=0 350 351 dhat(1:n)=den(1:n) 352 353 DEALLOCATE(den,a) 354 END SUBROUTINE hatL 355 356 !***********************************************************************88 357 ! on input: f1(i) and f2(i) are radial wfn * r for angular momentum l 358 ! on input: t1(i) and t2(i) are smooth radial wfn * r for angular momentum l 359 ! for r > rc, f1=t1, f2=t2 360 ! on output: qqqq is difference overlap matrix element 361 ! qqqq=<f1|f2>-<t1|t2> 362 !***********************************************************************88 363 SUBROUTINE dqij(Grid,PAW,ib,ic,qqqq) 364 TYPE(GridInfo), INTENT(IN) :: Grid 365 TYPE(PseudoInfo), INTENT(IN) :: PAW 366 INTEGER, INTENT(IN) :: ib,ic 367 REAL(8), INTENT(OUT) :: qqqq 368 369 INTEGER :: n,i,ok,irc 370 REAL(8) :: h 371 REAL(8), ALLOCATABLE :: dum(:) 372 373 qqqq=0 374 IF (PAW%l(ib)/=PAW%l(ic)) RETURN 375 n=Grid%n; h=Grid%h; irc=PAW%irc 376 ALLOCATE(dum(n),stat=ok) 377 IF (ok /=0) THEN 378 WRITE(STD_OUT,*) 'Error in dqij allocation', n,ok 379 STOP 380 ENDIF 381 DO i=1,n 382 dum(i)=PAW%ophi(i,ib)*PAW%ophi(i,ic)-PAW%otphi(i,ib)*PAW%otphi(i,ic) 383 ENDDO 384 qqqq=integrator(Grid,dum,1,irc) 385 386 DEALLOCATE(dum) 387 END SUBROUTINE dqij 388 389 !*********************************************************************** 390 ! SUBROUTINE dtij 391 ! on input: f1(i) and f2(i) are radial wfn * r for angular momentum l 392 ! on input: t1(i) and t2(i) are smooth radial wfn * r for angular momentum l 393 ! for r > rc, f1=t1, f2=t2 394 ! on output: tij is difference kinetic energy matrix element in Rydberg units 395 ! tij =<f1|T|f2>-<t1|T|t2> 396 !************************************************************************ 397 SUBROUTINE dtij(Grid,PAW,ib,ic,tij) 398 TYPE(GridInfo), INTENT(IN) :: Grid 399 TYPE(PseudoInfo), INTENT(IN) :: PAW 400 INTEGER, INTENT(IN) :: ib,ic 401 REAL(8), INTENT(OUT) :: tij 402 403 INTEGER :: n,i,ok,l,irc 404 REAL(8) :: angm 405 REAL(8), POINTER :: r(:) 406 REAL(8), ALLOCATABLE :: dum(:),del1(:),del2(:),tdel1(:),tdel2(:) 407 408 tij=0 409 IF (PAW%l(ib)/=PAW%l(ic)) RETURN 410 n=Grid%n; r=>Grid%r; l=PAW%l(ib); irc=PAW%irc 411 ALLOCATE(dum(n),del1(n),tdel1(n),del2(n),tdel2(n),stat=ok) 412 IF (ok /=0) THEN 413 WRITE(STD_OUT,*) 'Error in dtij allocation', n,ok 414 STOP 415 ENDIF 416 CALL derivative(Grid,PAW%ophi(:,ib),del1) 417 CALL derivative(Grid,PAW%ophi(:,ic),del2) 418 CALL derivative(Grid,PAW%otphi(:,ib),tdel1) 419 CALL derivative(Grid,PAW%otphi(:,ic),tdel2) 420 dum=0 ; angm=l*(l+1) 421 DO i=1,irc 422 dum(i)=del1(i)*del2(i)-tdel1(i)*tdel2(i) 423 ENDDO 424 del1=0;del2=0;tdel1=0;tdel2=0 425 del1(2:irc)=PAW%ophi(2:irc,ib)/Grid%r(2:irc) 426 del2(2:irc)=PAW%ophi(2:irc,ic)/Grid%r(2:irc) 427 tdel1(2:irc)=PAW%otphi(2:irc,ib)/Grid%r(2:irc) 428 tdel2(2:irc)=PAW%otphi(2:irc,ic)/Grid%r(2:irc) 429 DO i=1,irc 430 dum(i)=dum(i)+angm*(del1(i)*del2(i)-tdel1(i)*tdel2(i)) 431 ENDDO 432 tij=integrator(Grid,dum,1,irc) 433 434 DEALLOCATE(dum,del1,del2,tdel1,tdel2) 435 END SUBROUTINE dtij 436 437 !*********************************************************************** 438 ! SUBROUTINE altdtij 439 ! on input: f1(i) and f2(i) are radial wfn * r for angular momentum l 440 ! on input: t1(i) and t2(i) are smooth radial wfn * r for angular momentum l 441 ! for r > rc, f1=t1, f2=t2 442 ! on output: tij is difference kinetic energy matrix element in Rydberg units 443 ! tij =<f1|T|f2>-<t1|T|t2> 444 !************************************************************************ 445 SUBROUTINE altdtij(Grid,PAW,ib,ic,tij) 446 TYPE(GridInfo), INTENT(IN) :: Grid 447 TYPE(PseudoInfo), INTENT(IN) :: PAW 448 INTEGER, INTENT(IN) :: ib,ic 449 REAL(8), INTENT(OUT) :: tij 450 451 INTEGER :: n,i,ok,l,irc 452 REAL(8) :: angm 453 REAL(8), POINTER :: r(:) 454 REAL(8), ALLOCATABLE :: dum(:),tdel1(:),tdel2(:) 455 456 tij=0 457 IF (PAW%l(ib)/=PAW%l(ic)) RETURN 458 n=Grid%n; r=>Grid%r; l=PAW%l(ib); irc=PAW%irc 459 ALLOCATE(dum(n),tdel1(n),tdel2(n),stat=ok) 460 IF (ok /=0) THEN 461 WRITE(STD_OUT,*) 'Error in dtij allocation', n,ok 462 STOP 463 ENDIF 464 dum=0 465 DO i=2,irc 466 !dum(i)=(PAW%eig(ic)-AEPot%rv(i)/Grid%r(i))*PAW%ophi(i,ib)*PAW%ophi(i,ic) 467 dum(i)=PAW%ophi(i,ib)*PAW%Kop(i,ic) 468 ENDDO 469 CALL derivative(Grid,PAW%otphi(:,ic),tdel1) 470 CALL derivative(Grid,tdel1,tdel2) 471 angm=l*(l+1) 472 DO i=2,irc 473 dum(i)=dum(i)+PAW%otphi(i,ib)*(tdel2(i)-& 474& angm*PAW%otphi(i,ic)/(Grid%r(i)**2)) 475 ENDDO 476 tij=integrator(Grid,dum,1,irc) 477 478 DEALLOCATE(dum,tdel1,tdel2) 479 END SUBROUTINE altdtij 480 481 SUBROUTINE dvij(Grid,PAW,FC,nz,ib,ic,vij) 482 TYPE(GridInfo), INTENT(IN) :: Grid 483 TYPE(PseudoInfo), INTENT(IN) :: PAW 484 TYPE(FCInfo), INTENT(IN) :: FC 485 INTEGER, INTENT(IN) :: ib,ic 486 REAL(8), INTENT(IN) :: nz 487 REAL(8), INTENT(OUT) :: vij 488 489 INTEGER :: n,i,ok,irc 490 REAL(8) :: h,en,q,qt 491 REAL(8), ALLOCATABLE :: dum(:),d1(:) 492 REAL(8), POINTER :: r(:) 493 494 vij=0 495 IF (PAW%l(ib)/=PAW%l(ic)) RETURN 496 n=Grid%n; h=Grid%h; r=>Grid%r; irc=PAW%irc 497 ALLOCATE(dum(n),d1(n),stat=ok) 498 IF (ok /=0) THEN 499 WRITE(STD_OUT,*) 'Error in dvij allocation', n,ok 500 STOP 501 ENDIF 502 503 q=integrator(Grid,FC%coreden) 504 WRITE(STD_OUT,*) 'core electrons ',q,FC%zcore 505 CALL poisson(Grid,q,FC%coreden,dum,en) 506 dum=dum-2*nz 507 qt=integrator(Grid,PAW%tcore) 508 WRITE(STD_OUT,*) 'coretail electrons ',qt 509 CALL poisson(Grid,qt,PAW%tcore,d1,en) 510 dum(1)=0 511 DO i=2,irc 512 dum(i)=PAW%ophi(i,ib)*PAW%ophi(i,ic)*dum(i)/r(i)-& 513& PAW%otphi(i,ib)*PAW%otphi(i,ic)*(PAW%vloc(i)+d1(i)/r(i)) 514 ENDDO 515 vij=integrator(Grid,dum,1,irc) 516 517 DEALLOCATE(dum) 518 END SUBROUTINE dvij 519 520 SUBROUTINE avij(Grid,PAW,ib,ic,vij) 521 TYPE(GridInfo), INTENT(IN) :: Grid 522 TYPE(PseudoInfo), INTENT(IN) :: PAW 523 INTEGER, INTENT(IN) :: ib,ic 524 REAL(8), INTENT(OUT) :: vij 525 526 INTEGER :: n,i,ok,irc 527 REAL(8) :: h,en,q,qt 528 REAL(8), ALLOCATABLE :: dum(:),d1(:) 529 REAL(8), POINTER :: r(:) 530 531 vij=0 532 IF (PAW%l(ib)/=PAW%l(ic)) RETURN 533 n=Grid%n; h=Grid%h; r=>Grid%r; irc=PAW%irc 534 ALLOCATE(dum(n),d1(n),stat=ok) 535 IF (ok /=0) THEN 536 WRITE(STD_OUT,*) 'Error in dvij allocation', n,ok 537 STOP 538 ENDIF 539 540 dum(1)=0 541 DO i=2,irc 542 dum(i)=(PAW%ophi(i,ib)*PAW%ophi(i,ic)*PAW%AErefrv(i)-& 543& PAW%otphi(i,ib)*PAW%otphi(i,ic)*PAW%rveff(i))/r(i) 544 ENDDO 545 vij=integrator(Grid,dum,1,irc) 546 547 DEALLOCATE(dum) 548 END SUBROUTINE avij 549 550 !******************************************************** 551 ! SUBROUTINE calcwij 552 ! subroutine to accumulate the wij coefficience for an input 553 ! smooth wavefunction twfn and occupancy and l 554 !******************************************************** 555 SUBROUTINE calcwij(Grid,PAW,l,occ,twfn,wij) 556 TYPE(GridInfo), INTENT(IN) :: Grid 557 TYPE(PseudoInfo), INTENT(IN) :: PAW 558 INTEGER, INTENT(IN) :: l 559 REAL(8), INTENT(IN) :: twfn(:),occ 560 REAL(8), INTENT(INOUT) :: wij(:,:) 561 562 INTEGER :: n,i,ok,irc,ib,ic,nbase 563 REAL(8) :: h 564 REAL(8), ALLOCATABLE :: bm(:) 565 566 n=Grid%n; h=Grid%h ; irc=PAW%irc 567 568 nbase=PAW%nbase 569 570 ALLOCATE(bm(nbase)) 571 572 bm=0 573 DO ib=1,nbase 574 IF (l==PAW%l(ib)) bm(ib)=overlap(Grid,PAW%otp(:,ib),twfn,1,irc) 575 !IF (l==PAW%l(ib)) write(std_out,*) 'accum wij',l,ib,bm(ib) 576 ENDDO 577 DO ib=1,nbase 578 DO ic=1,nbase 579 wij(ib,ic)=wij(ib,ic) + occ*bm(ib)*bm(ic) 580 ENDDO 581 ENDDO 582 DEALLOCATE(bm) 583 END SUBROUTINE calcwij 584 585 SUBROUTINE FillHat(Grid,PAW) 586 TYPE(GridInfo) , INTENT(IN):: Grid 587 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 588 589 INTEGER :: ll,n,l 590 591 ll=MAXVAL(PAW%TOCCWFN%l(:)); ll=MAX(ll,PAW%lmax); ll=2*ll 592 n=Grid%n 593 594 ALLOCATE(PAW%g(n,ll+1)) 595 596 DO l=0,ll 597 CALL hatL(Grid,PAW,l,PAW%g(:,l+1)) 598 ENDDO 599 600 END SUBROUTINE FillHat 601 602!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 603 !! Get_Energy_EXX_smoothpseudo !!!! 604!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 605 SUBROUTINE Get_Energy_EXX_pseudo(Grid,PAW,eex) 606 TYPE(gridinfo), INTENT(IN) :: Grid 607 TYPE(pseudoinfo), INTENT(IN) :: PAW 608 REAL(8), INTENT(OUT) :: eex 609 610 REAL(8), POINTER :: phi(:,:),r(:) 611 REAL(8), ALLOCATABLE :: wfp(:),dum(:),ddum(:),ch(:),ar(:) 612 REAL(8), ALLOCATABLE :: Fnu(:,:),Lnu(:,:),Snu(:),Mnunup(:,:),Cnu(:) 613 INTEGER :: i,j,k,n,m,l,ll,li,lj,lmin,lmax,io,jo,ok,norbit,last,nodes 614 INTEGER :: nu,nup,ni,nj 615 REAL(8) :: occ,term,term1,wgt,occi,occj,rc,kappa,test 616 REAL(8), PARAMETER :: threshold=1.d-8 617 618 n=Grid%n 619 norbit=PAW%TOCCwfn%norbit 620 phi=>PAW%TOCCwfn%wfn 621 r=>grid%r 622 623 ALLOCATE(wfp(n),dum(n),ddum(n),ch(n),ar(n),stat=i) 624 IF (i/=0) THEN 625 WRITE(STD_OUT,*) 'allocation error in Get_Energy_EXX', i,n 626 STOP 627 ENDIF 628 629 eex=0;test=0 630 DO io=1,norbit 631 occ=PAW%TOCCwfn%occ(io) 632 IF (occ>threshold.AND.(.NOT.PAW%TOCCwfn%iscore(io))) THEN 633 li=PAW%TOCCwfn%l(io);ni=PAW%TOCCwfn%np(io) 634 DO jo=1,norbit 635 occj=PAW%TOCCwfn%occ(jo) 636 IF (occj>threshold.AND.(.NOT.PAW%TOCCwfn%iscore(jo))) THEN 637 lj=PAW%TOCCwfn%l(jo);nj=PAW%TOCCwfn%np(jo) 638 lmax=li+lj 639 lmin=ABS(li-lj) 640 wfp(1:n)=phi(1:n,io)*phi(1:n,jo) 641 ar(1:n)=PAW%OCCwfn%wfn(1:n,io)*PAW%OCCwfn%wfn(1:n,jo) 642 DO ll=lmin,lmax,2 643 CALL EXXwgt(occ,occj,io,li,jo,lj,ll,wgt) 644 IF (wgt>threshold) THEN 645 wgt=wgt*(2*ll+1) ! because of apoisson convention 646 call Calc_Moment(Grid,PAW,phi(1:n,io),phi(1:n,jo),& 647& li,lj,ll,dum) 648 ddum=wfp+dum 649 CALL apoisson(Grid,ll,n,ddum,dum) 650 CALL apoisson(Grid,ll,n,ar,ch) 651 ddum(2:n)=(dum(2:n)*ddum(2:n))/r(2:n) 652 ddum(1)=0 653 eex=eex-wgt*integrator(Grid,ddum)/2 654 ch(2:n)=(ch(2:n)*ar(2:n))/r(2:n); ch(1)=0 655 test=test-wgt*integrator(Grid,ch)/2 656 ENDIF 657 ENDDO 658 ENDIF 659 ENDDO !jo 660 661 ENDIF 662 ENDDO 663 664 write(std_out,*) 'Get_Energy_EXX_pseudo', eex,test 665 DEALLOCATE(wfp,dum,ddum,ch,ar) 666 667 END SUBROUTINE Get_Energy_EXX_pseudo 668 669 670!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 671 !! Calc_Moment 672!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 673 SUBROUTINE Calc_Moment(Grid,PAW,wfn1,wfn2,l1,l2,ll,m) 674 TYPE(gridinfo), INTENT(IN) :: Grid 675 TYPE(pseudoinfo), INTENT(IN) :: PAW 676 REAL(8), INTENT(IN) :: wfn1(:),wfn2(:) 677 INTEGER, INTENT(IN) :: l1,l2,ll 678 REAL(8), INTENT(INOUT) :: m(:) 679 680 REAL(8) :: dot1,dot2 681 INTEGER :: i,j,k,n,nbasis,ib,jb 682 REAL(8), PARAMETER :: threshold=1.d-8 683 REAL(8), ALLOCATABLE :: dum(:) 684 685 nbasis=PAW%nbase 686 n=Grid%n 687 ALLOCATE(dum(n)) 688 689 m=0 690 DO ib=1,nbasis 691 if (PAW%l(ib)==l1) then 692 dot1=overlap(Grid,wfn1,PAW%otp(:,ib)); write(std_out,*)'dot1',dot1 693 do jb=1,nbasis 694 if (PAW%l(jb)==l2) then 695 dot2=overlap(Grid,wfn2,PAW%otp(:,jb)); write(std_out,*)'dot2',dot2 696 m=m+dot1*dot2*PAW%mLij(ib,jb,ll+1)*PAW%g(:,ll+1) 697 endif 698 enddo 699 endif 700 Enddo 701 702 dum=m*(Grid%r**ll) 703 write(std_out,*)'Check moment',l1,l2,ll,integrator(Grid,dum) 704 deallocate(dum) 705 END SUBROUTINE Calc_Moment 706 707 708!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 709 !! Get_Energy_EXX_onecenter_pseudo !!!! 710!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 711 SUBROUTINE Get_Energy_EXX_pseudo_one(Grid,PAW,eex) 712 TYPE(gridinfo), INTENT(IN) :: Grid 713 TYPE(pseudoinfo), INTENT(IN) :: PAW 714 REAL(8), INTENT(OUT) :: eex 715 716 REAL(8), POINTER :: phi(:,:),r(:) 717 INTEGER :: i,j,k,n,m,l,ll,li,lj,lmin,lmax,io,jo,ok,norbit,last,nodes 718 INTEGER :: nu,nup,ni,nj,ib,jb,kb,lb,nbasis 719 REAL(8) :: occ,term,term1,wgt,occi,occj,rc,kappa,doti,dotj,dotk,dotl 720 REAL(8), PARAMETER :: threshold=1.d-8 721 722 n=Grid%n 723 norbit=PAW%TOCCwfn%norbit 724 phi=>PAW%TOCCwfn%wfn 725 r=>grid%r 726 nbasis=PAW%nbase 727 728 729 eex=0 730 DO io=1,norbit 731 occ=PAW%TOCCwfn%occ(io) 732 IF (occ>threshold.AND.(.NOT.PAW%TOCCwfn%iscore(io))) THEN 733 li=PAW%TOCCwfn%l(io);ni=PAW%TOCCwfn%np(io) 734 DO jo=1,norbit 735 occj=PAW%TOCCwfn%occ(jo) 736 IF (occj>threshold.AND.(.NOT.PAW%TOCCwfn%iscore(jo))) THEN 737 lj=PAW%TOCCwfn%l(jo);nj=PAW%TOCCwfn%np(jo) 738 lmax=li+lj 739 lmin=ABS(li-lj) 740 DO ib=1,nbasis 741 if (PAW%l(ib)==li) then 742 doti=overlap(Grid,PAW%otp(:,ib),phi(:,io)) 743 Do jb=1,nbasis 744 If (PAW%l(jb)==lj) then 745 dotj=overlap(Grid,PAW%otp(:,jb),phi(:,jo)) 746 Do kb=1,nbasis 747 If (PAW%l(kb)==lj) then 748 dotk=overlap(Grid,PAW%otp(:,kb),phi(:,jo)) 749 Do lb=1,nbasis 750 If (PAW%l(lb)==li) then 751 dotl=overlap(Grid,PAW%otp(:,lb),phi(:,io)) 752 term=doti*dotj*dotk*dotl 753 DO ll=lmin,lmax,2 754 CALL EXXwgt(occ,occj,io,li,jo,lj,ll,wgt) 755 IF (wgt>threshold) THEN 756 wgt=wgt*(2*ll+1) !Not in DR 757 eex=eex-wgt*term*PAW%DR(ib,jb,kb,lb,ll+1)/2 758 endif 759 enddo !ll 760 endif 761 enddo !ib 762 endif 763 enddo !kb 764 endif 765 enddo !jb 766 endif 767 enddo !ib 768 endif 769 enddo !jo 770 endif 771 enddo !io 772 773 write(std_out,*) 'Get_Energy_EXX_pseudo_one val-val: ', eex 774 775 If (PAW%ncoreshell>0) then 776 do k=1,PAW%ncoreshell 777 do ib=1,nbasis 778 do jb=1,nbasis 779 if (PAW%l(ib)==PAW%l(jb)) then 780 eex=eex-PAW%wij(ib,jb)*PAW%DRVC(k,ib,jb) 781 endif 782 enddo 783 enddo 784 enddo 785 write(std_out,*) 'Get_Energy_EXX_pseudo_one corrected: ', eex 786 endif 787 END SUBROUTINE Get_Energy_EXX_pseudo_one 788 789 790 791!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 792! Fill stored matrix elements and calculate atomic energy 793!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 794 SUBROUTINE SPMatrixElements(Grid,Pot,FC,PAW) 795 TYPE(GridInfo) , INTENT(IN):: Grid 796 TYPE(PotentialInfo), INTENT(IN) :: Pot 797 TYPE(FCInfo), INTENT(IN) :: FC 798 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 799 800 INTEGER :: io,jo,ko,lo,ll,lmax,i,j,k,l,n,norbit,nbase,ib,jb,kb,lb,irc 801 INTEGER :: llmin, llmax,lp,lr 802 REAL(8), ALLOCATABLE :: arg(:),dum(:),rh(:),trh(:),d(:),td(:),wij(:,:) 803 REAL(8) :: x,y,z,rr,occ 804 REAL(8) :: Qcore,tQcore 805 LOGICAL :: even 806 807 n=Grid%n ; irc=PAW%irc; lr=min(irc+20,n) ; nbase=PAW%nbase 808 ll=MAXVAL(PAW%TOCCWFN%l(:)); ll=MAX(ll,PAW%lmax); ll=2*ll 809 ALLOCATE(arg(n),dum(n),rh(n),trh(n),d(n),td(n)) 810 arg=0; dum=0 811 ALLOCATE(PAW%mLij(nbase,nbase,ll+1),PAW%DR(nbase,nbase,nbase,nbase,ll+1)) 812 ALLOCATE(wij(nbase,nbase)) 813 PAW%mLij=0.d0;PAW%DR=0.d0 814 815 CALL poisson(Grid,Qcore,PAW%core,dum,y) 816 PAW%rVf=-2*Pot%nz+dum 817 818 CALL poisson(Grid,tQcore,PAW%tcore,dum,y) 819 PAW%rtVf=Grid%r*PAW%vloc + dum + & 820& (-Pot%nz + Qcore - tQcore)*PAW%hatpot 821 !DO l=0,ll 822 ! CALL hatL(Grid,PAW,l,PAW%g(:,l+1)) 823 !ENDDO 824 825 OPEN(1001,file='rvf',form='formatted') 826 DO i=1,n 827 WRITE(1001,'(1p,30e15.7)') Grid%r(i),PAW%rVf(i),PAW%rtVf(i),& 828& PAW%hatpot(i), PAW%hatden(i),(PAW%g(i,l+1),l=0,ll) 829 ENDDO 830 CLOSE(1001) 831 832 arg=PAW%tcore+(-Pot%nz + Qcore - tQcore)*PAW%hatden 833 CALL poisson(Grid,y,arg,dum,PAW%Eaion) 834 write(std_out,*) 'Eaion ', y, PAW%Eaion 835 arg=dum*PAW%hatden; arg(1)=0; arg(2:n)=arg(2:n)/Grid%r(2:n) 836 PAW%Eaionhat=integrator(Grid,arg) 837 write(std_out,*) 'Eaionhat ', PAW%Eaionhat 838 839 DO ib=1,nbase 840 l=PAW%l(ib) 841 DO jb=1,nbase 842 IF (PAW%l(jb)==l) THEN 843 !CALL kinetic_ij(Grid,PAW%ophi(:,ib),PAW%ophi(:,jb),l,x,lr) 844 !CALL kinetic_ij(Grid,PAW%otphi(:,ib),PAW%otphi(:,jb),l,y,lr) 845 If (scalarrelativistic) then 846 call altdtij(Grid,PAW,ib,jb,x) 847 Else 848 CALL deltakinetic_ij(Grid,PAW%ophi(:,ib),PAW%ophi(:,jb), & 849& PAW%otphi(:,ib),PAW%otphi(:,jb),l,x,PAW%irc) 850 Endif 851 WRITE(STD_OUT,'(" Kinetic ", 3i5, 1p,3e15.7)') ib,jb,l,x 852 PAW%Kij(ib,jb)=x 853 dum=PAW%ophi(:,ib)*PAW%ophi(:,jb)*PAW%rVf - & 854& PAW%otphi(:,ib)*PAW%otphi(:,jb)*PAW%rtVf 855 dum(2:n)=dum(2:n)/Grid%r(2:n) 856 dum(1)=0 857 PAW%VFij(ib,jb)=integrator(Grid,dum,1,lr) 858 WRITE(STD_OUT,'(" Fix pot ", 3i5, 1p,3e15.7)') ib,jb,l,PAW%VFij(ib,jb) 859 ENDIF 860 ENDDO 861 ENDDO 862 863 ! Average equivalent terms 864 DO ib=1,nbase 865 DO jb=ib,nbase 866 IF(jb>ib) THEN 867 x=PAW%Kij(ib,jb); y=PAW%Kij(jb,ib) 868 x=0.5d0*(x+y) 869 PAW%Kij(ib,jb)=x; PAW%Kij(jb,ib)=x 870 x=PAW%VFij(ib,jb); y=PAW%VFij(jb,ib) 871 x=0.5d0*(x+y) 872 PAW%VFij(ib,jb)=x; PAW%VFij(jb,ib)=x 873 ENDIF 874 ENDDO 875 ENDDO 876 877 DO ib=1,nbase 878 DO jb=1,nbase 879 llmin=ABS(PAW%l(ib)-PAW%l(jb)) 880 llmax=PAW%l(ib)+PAW%l(jb) 881 DO l=llmin,llmax,2 882 DO i=1,irc 883 rr=Grid%r(i) 884 dum(i)=(rr**l)*(PAW%ophi(i,ib)*PAW%ophi(i,jb) & 885& -PAW%otphi(i,ib)*PAW%otphi(i,jb)) 886 ENDDO 887 PAW%mLij(ib,jb,l+1)=integrator(Grid,dum,1,irc) 888 WRITE(STD_OUT,'("mLij ",3i5,1p,1e15.7)') ib,jb,l,PAW%mLij(ib,jb,l+1) 889 ENDDO 890 ENDDO 891 ENDDO 892 893 ! Average equivalent terms 894 DO ib=1,nbase 895 DO jb=ib,nbase 896 IF (jb>ib) THEN 897 llmin=ABS(PAW%l(ib)-PAW%l(jb)) 898 llmax=PAW%l(ib)+PAW%l(jb) 899 DO l=llmin,llmax,2 900 x=PAW%mLij(ib,jb,l+1); y=PAW%mLij(jb,ib,l+1) 901 x=0.5d0*(x+y) 902 PAW%mLij(ib,jb,l+1)=x; PAW%mLij(jb,ib,l+1)=x 903 ENDDO 904 ENDIF 905 ENDDO 906 ENDDO 907 908 WRITE(STD_OUT,*) 'DR matrix elements ' 909 DO ib=1,nbase 910 DO jb=1,nbase 911 d=PAW%ophi(:,ib)*PAW%ophi(:,jb) 912 td=PAW%otphi(:,ib)*PAW%otphi(:,jb) 913 llmin=ABS(PAW%l(ib)-PAW%l(jb)) 914 llmax=PAW%l(ib)+PAW%l(jb) 915 DO l=llmin,llmax,2 916 CALL apoisson(Grid,l,lr,d,rh) 917 arg=td+PAW%mLij(ib,jb,l+1)*PAW%g(:,l+1) 918 CALL apoisson(Grid,l,lr,arg,trh) 919 DO kb=1,nbase 920 DO lb=1,nbase 921 lp=llmax+PAW%l(kb)+PAW%l(lb) 922 even=.FALSE. 923 IF (2*(lp/2)==lp) even=.TRUE. 924 IF (even.AND.l.GE.ABS(PAW%l(kb)-PAW%l(lb)).AND. & 925& l.LE.(PAW%l(kb)+PAW%l(lb))) THEN 926 arg=PAW%otphi(:,kb)*PAW%otphi(:,lb)+& 927& PAW%mLij(kb,lb,l+1)*PAW%g(:,l+1) 928 dum(1:lr)=PAW%ophi(1:lr,kb)*PAW%ophi(1:lr,lb)*rh(1:lr)& 929& - arg(1:lr)*trh(1:lr) 930 dum(1)=0; dum(2:lr)=dum(2:lr)/Grid%r(2:lr) 931 PAW%DR(ib,jb,kb,lb,l+1)=integrator(Grid,dum,1,lr) 932 !if (abs(PAW%DR(ib,jb,kb,lb,l+1))>100.d0) then 933 ! write(std_out,*) 'problem', ib,jb,kb,lb,l+1,PAW%DR(ib,jb,kb,lb,l+1),& 934 !& PAW%mLij(kb,lb,l+1) 935 ! open(unit=1001,file='stuff',form='formatted') 936 ! do i=1,lr 937 ! write(1001,'(1p,20e15.7)') & 938 !& Grid%r(i),PAW%ophi(i,kb),PAW%ophi(i,lb), & 939 !& PAW%otphi(i,kb),PAW%otphi(i,lb),PAW%g(i,l+1),rh(i),trh(i) 940 ! enddo 941 ! close(1001) 942 ! stop 943 ! endif 944 WRITE(STD_OUT,'(5i5,1p,1e20.10)') ib,jb,kb,lb,l, & 945& PAW%DR(ib,jb,kb,lb,l+1) 946 ENDIF 947 ENDDO !lb 948 ENDDO !kb 949 ENDDO !l 950 ENDDO !jb 951 ENDDO !ib 952 953 ! Average equivalent terms 954 DO ib=1,nbase 955 DO jb=ib,nbase 956 llmin=ABS(PAW%l(ib)-PAW%l(jb)) 957 llmax=PAW%l(ib)+PAW%l(jb) 958 DO l=llmin,llmax,2 959 DO kb=1,nbase 960 DO lb=kb,nbase 961 lp=llmax+PAW%l(kb)+PAW%l(lb) 962 even=.FALSE. 963 IF (2*(lp/2)==lp) even=.TRUE. 964 IF (even.AND.l.GE.ABS(PAW%l(kb)-PAW%l(lb)).AND. & 965& l.LE.(PAW%l(kb)+PAW%l(lb))) THEN 966 arg(1)=PAW%DR(ib,jb,kb,lb,l+1) 967 arg(2)=PAW%DR(ib,jb,lb,kb,l+1) 968 arg(3)=PAW%DR(jb,ib,kb,lb,l+1) 969 arg(4)=PAW%DR(jb,ib,lb,kb,l+1) 970 x=0.25d0*(arg(1)+arg(2)+arg(3)+arg(4)) 971 PAW%DR(ib,jb,kb,lb,l+1)=x 972 PAW%DR(ib,jb,lb,kb,l+1)=x 973 PAW%DR(jb,ib,kb,lb,l+1)=x 974 PAW%DR(jb,ib,lb,kb,l+1)=x 975 ENDIF 976 ENDDO !lb 977 ENDDO !kb 978 ENDDO !l 979 ENDDO !jb 980 ENDDO !ib 981 982! If (TRIM(PAW%exctype)=="EXXKLI".OR.TRIM(PAW%exctype)=="HF") & 983!& Call COREVAL_EXX(Grid,PAW) 984 985! COREVAL_EXX now called in order to output core-valence exchange terms 986! for possible use in exact exchange and hybrid treatments 987 988 Call COREVAL_EXX(Grid,PAW) 989 Call CORECORE_EXX(Grid,PAW) 990 991 IF (TRIM(PAW%exctype)=="HF") THEN 992 norbit=PAW%TOCCWFN%norbit 993 ALLOCATE(PAW%mLic(nbase,norbit,ll+1),PAW%DRC(nbase,nbase,norbit,ll+1),& 994& PAW%mLcc(norbit,norbit,ll+1),PAW%DRCC(nbase,norbit,norbit,ll+1),& 995& PAW%DRCjkl(norbit,nbase,nbase,nbase,ll+1),& 996& PAW%Dcj(norbit,nbase),PAW%TOCCWFN%lqp(norbit,norbit)) 997 PAW%mLic=0.d0; PAW%DRC=0.d0; PAW%mLcc=0.d0; PAW%DRCC=0.d0; PAW%Dcj=0.d0 998 PAW%DRCjkl=0.d0; PAW%TOCCWFN%lqp=0.d0 999 lr=MAX(lr,PAW%coretailpoints) 1000 write(std_out,*) 'lr for core treatment ', lr 1001 DO io=1,norbit 1002 IF (PAW%TOCCWFN%iscore(io)) THEN 1003 DO jb=1,nbase 1004 llmin=ABS(PAW%TOCCWFN%l(io)-PAW%l(jb)) 1005 llmax=PAW%TOCCWFN%l(io)+PAW%l(jb) 1006 DO l=llmin,llmax,2 1007 DO i=1,n 1008 rr=Grid%r(i) 1009 dum(i)=(rr**l)*(PAW%OCCWFN%wfn(i,io)*PAW%ophi(i,jb) & 1010& -PAW%TOCCWFN%wfn(i,io)*PAW%otphi(i,jb)) 1011 ENDDO 1012 PAW%mLic(jb,io,l+1)=integrator(Grid,dum,1,lr) 1013 WRITE(STD_OUT,'("mLic ",3i5,1p,1e15.7)') jb,io,l,PAW%mLic(jb,io,l+1) 1014 ENDDO 1015 ENDDO 1016 ENDIF 1017 ENDDO 1018 1019 DO io=1,norbit 1020 IF (PAW%TOCCWFN%iscore(io)) THEN 1021 DO jo=1,norbit 1022 IF (PAW%TOCCWFN%iscore(jo)) THEN 1023 llmin=ABS(PAW%TOCCWFN%l(io)-PAW%TOCCWFN%l(jo)) 1024 llmax=PAW%TOCCWFN%l(io)+PAW%TOCCWFN%l(jo) 1025 DO l=llmin,llmax,2 1026 DO i=1,n 1027 rr=Grid%r(i) 1028 dum(i)=(rr**l)*& 1029& (PAW%OCCWFN%wfn(i,io)*PAW%OCCWFN%wfn(i,jo) & 1030& -PAW%TOCCWFN%wfn(i,io)*PAW%TOCCWFN%wfn(i,jo)) 1031 ENDDO 1032 PAW%mLcc(jo,io,l+1)=integrator(Grid,dum,1,lr) 1033 WRITE(STD_OUT,'("mLcc ",3i5,1p,1e15.7)') jo,io,l,& 1034& PAW%mLcc(jo,io,l+1) 1035 ENDDO 1036 ENDIF 1037 ENDDO 1038 ENDIF 1039 ENDDO 1040 1041 WRITE(STD_OUT,*) 'DRC matrix elements ' 1042 DO io=1,norbit 1043 IF (PAW%TOCCWFN%iscore(io)) THEN 1044 DO ib=1,nbase 1045 d=PAW%ophi(:,ib)*PAW%OCCWFN%wfn(:,io) 1046 td=PAW%otphi(:,ib)*PAW%TOCCWFN%wfn(:,io) 1047 llmin=ABS(PAW%l(ib)-PAW%TOCCWFN%l(io)) 1048 llmax=PAW%l(ib)+PAW%TOCCWFN%l(io) 1049 DO l=llmin,llmax,2 1050 CALL apoisson(Grid,l,lr,d,rh) 1051 arg=td+PAW%mLic(ib,io,l+1)*PAW%g(:,l+1) 1052 CALL apoisson(Grid,l,n,arg,trh) 1053 DO jb=1,nbase 1054 lp=llmax+PAW%l(jb)+PAW%TOCCWFN%l(io) 1055 even=.FALSE. 1056 IF (2*(lp/2)==lp) even=.TRUE. 1057 IF (even.AND.l.GE.ABS(PAW%l(jb)-PAW%TOCCWFN%l(io)).AND. & 1058& l.LE.(PAW%l(jb)+PAW%TOCCWFN%l(io))) THEN 1059 arg=PAW%otphi(:,jb)*PAW%TOCCWFN%wfn(:,io) +& 1060& PAW%mLic(jb,io,l+1)*PAW%g(:,l+1) 1061 dum=PAW%ophi(:,jb)*PAW%OCCWFN%wfn(:,io)*rh - arg*trh 1062 dum(1)=0; dum(2:n)=dum(2:n)/Grid%r(2:n) 1063 PAW%DRC(ib,jb,io,l+1)=integrator(Grid,dum,1,lr) 1064 WRITE(STD_OUT,'(4i5,1p,1e20.10)') ib,jb,io,l, & 1065& PAW%DRC(ib,jb,io,l+1) 1066 ENDIF 1067 ENDDO !jb 1068 ENDDO !l 1069 ENDDO !ib 1070 ENDIF 1071 ENDDO !io 1072 1073 ! Average equivalent terms 1074 DO io=1,norbit 1075 IF (PAW%TOCCWFN%iscore(io)) THEN 1076 DO ib=1,nbase 1077 llmin=ABS(PAW%l(ib)-PAW%TOCCWFN%l(io)) 1078 llmax=PAW%l(ib)+PAW%TOCCWFN%l(io) 1079 DO l=llmin,llmax,2 1080 DO jb=ib,nbase 1081 lp=llmax+PAW%l(jb)+PAW%TOCCWFN%l(io) 1082 even=.FALSE. 1083 IF (2*(lp/2)==lp) even=.TRUE. 1084 IF (even.AND.l.GE.ABS(PAW%l(jb)-PAW%TOCCWFN%l(io)).AND. & 1085& l.LE.(PAW%l(jb)+PAW%TOCCWFN%l(io))) THEN 1086 arg(1)=PAW%DRC(ib,jb,io,l+1) 1087 arg(2)=PAW%DRC(jb,ib,io,l+1) 1088 x=0.5d0*(arg(1)+arg(2)) 1089 PAW%DRC(ib,jb,io,l+1)=x 1090 PAW%DRC(jb,ib,io,l+1)=x 1091 ENDIF 1092 ENDDO 1093 ENDDO 1094 ENDDO 1095 ENDIF 1096 ENDDO 1097 1098 WRITE(STD_OUT,*) 'DRCC matrix elements ' 1099 DO io=1,norbit 1100 IF (PAW%TOCCWFN%iscore(io)) THEN 1101 DO jo=1,norbit 1102 IF (PAW%TOCCWFN%iscore(jo)) THEN 1103 d=PAW%OCCWFN%wfn(:,jo)*PAW%OCCWFN%wfn(:,io) 1104 td=PAW%TOCCWFN%wfn(:,jo)*PAW%TOCCWFN%wfn(:,io) 1105 llmin=ABS(PAW%TOCCWFN%l(jo)-PAW%TOCCWFN%l(io)) 1106 llmax=PAW%TOCCWFN%l(jo)+PAW%TOCCWFN%l(io) 1107 DO l=llmin,llmax,2 1108 CALL apoisson(Grid,l,lr,d,rh) 1109 arg=td+PAW%mLcc(jo,io,l+1)*PAW%g(:,l+1) 1110 CALL apoisson(Grid,l,lr,arg,trh) 1111 DO jb=1,nbase 1112 lp=llmax+PAW%l(jb)+PAW%TOCCWFN%l(io) 1113 even=.FALSE. 1114 IF (2*(lp/2)==lp) even=.TRUE. 1115 IF (even.AND.l.GE.ABS(PAW%l(jb)-PAW%TOCCWFN%l(io)) & 1116& .AND.l.LE.(PAW%l(jb)+PAW%TOCCWFN%l(io))) THEN 1117 arg=PAW%otphi(:,jb)*PAW%TOCCWFN%wfn(:,io) +& 1118& PAW%mLic(jb,io,l+1)*PAW%g(:,l+1) 1119 dum=PAW%ophi(:,jb)*PAW%OCCWFN%wfn(:,io)*rh - arg*trh 1120 dum(1)=0; dum(2:n)=dum(2:n)/Grid%r(2:n) 1121 PAW%DRCC(jb,jo,io,l+1)=integrator(Grid,dum,1,lr) 1122 WRITE(STD_OUT,'(4i5,1p,1e20.10)') jb,jo,io,l,PAW%DRCC(jb,jo,io,l+1) 1123 ENDIF 1124 ENDDO !jb 1125 ENDDO !l 1126 ENDIF 1127 ENDDO !jo 1128 ENDIF 1129 ENDDO !io 1130 1131 WRITE(STD_OUT,*) 'Dcj matrix elements ' 1132 DO io=1,norbit 1133 IF (PAW%TOCCWFN%iscore(io)) THEN 1134 l=PAW%TOCCWFN%l(io) 1135 DO jb=1,nbase 1136 IF (PAW%l(jb)==l) THEN 1137 !CALL kinetic_ij(Grid,PAW%OCCWFN%wfn(:,io),& 1138 !& PAW%ophi(:,jb),l,x,lr) 1139 !CALL kinetic_ij(Grid,PAW%TOCCWFN%wfn(:,io),& 1140 !& PAW%otphi(:,jb),l,y,lr) 1141 CALL deltakinetic_ij(Grid,PAW%OCCWFN%wfn(:,io),PAW%ophi(:,jb), & 1142& PAW%TOCCWFN%wfn(:,io),PAW%otphi(:,jb),l,x,PAW%irc) 1143 !WRITE(STD_OUT,'(" Kinetic ", 3i5, 1p,3e15.7)') io,jb,l,x,y,x-y 1144 PAW%Dcj(io,jb)=x 1145 dum=PAW%OCCWFN%wfn(:,io)*PAW%ophi(:,jb)*PAW%rVf - & 1146& PAW%TOCCWFN%wfn(:,io)*PAW%otphi(:,jb)*PAW%rtVf 1147 dum(2:n)=dum(2:n)/Grid%r(2:n) 1148 dum(1)=0 1149 x=integrator(Grid,dum,1,lr) 1150 !WRITE(STD_OUT,'(" Fix pot ", 3i5, 1p,3e15.7)') io,jb,l,x 1151 PAW%Dcj(io,jb)=PAW%Dcj(io,jb)+x 1152 WRITE(STD_OUT,'(3i5,1p,e15.7)') io,jb,l,PAW%Dcj(io,jb) 1153 ENDIF 1154 ENDDO 1155 ENDIF 1156 ENDDO 1157 1158 WRITE(STD_OUT,*) 'DRCjkl matrix elements ' 1159 DO io=1,norbit 1160 IF (PAW%TOCCWFN%iscore(io)) THEN 1161 DO jb=1,nbase 1162 d=PAW%OCCWFN%wfn(:,io)*PAW%ophi(:,jb) 1163 td=PAW%TOCCWFN%wfn(:,io)*PAW%otphi(:,jb) 1164 llmin=ABS(PAW%TOCCWFN%l(io)-PAW%l(jb)) 1165 llmax=PAW%TOCCWFN%l(io)+PAW%l(jb) 1166 DO l=llmin,llmax,2 1167 CALL apoisson(Grid,l,lr,d,rh) 1168 arg=td+PAW%mLic(jb,io,l+1)*PAW%g(:,l+1) 1169 CALL apoisson(Grid,l,lr,arg,trh) 1170 DO kb=1,nbase 1171 DO lb=1,nbase 1172 lp=llmax+PAW%l(kb)+PAW%l(lb) 1173 even=.FALSE. 1174 IF (2*(lp/2)==lp) even=.TRUE. 1175 IF (even.AND.l.GE.ABS(PAW%l(kb)-PAW%l(lb)).AND. & 1176& l.LE.(PAW%l(kb)+PAW%l(lb))) THEN 1177 arg=PAW%otphi(:,kb)*PAW%otphi(:,lb)+& 1178& PAW%mLij(kb,lb,l+1)*PAW%g(:,l+1) 1179 dum=PAW%ophi(:,kb)*PAW%ophi(:,lb)*rh - arg*trh 1180 dum(1)=0; dum(2:n)=dum(2:n)/Grid%r(2:n) 1181 PAW%DRCjkl(io,jb,kb,lb,l+1)=integrator(Grid,dum,1,lr) 1182 WRITE(STD_OUT,'(5i5,1p,1e20.10)') io,jb,kb,lb,l, & 1183& PAW%DRCjkl(io,jb,kb,lb,l+1) 1184 ENDIF 1185 ENDDO !lb 1186 ENDDO !kb 1187 ENDDO !l 1188 ENDDO !jb 1189 ENDIF 1190 ENDDO !io 1191 ENDIF 1192 1193 WRITE(STD_OUT,*) ' Completed SPMatrixElements'; CALL flush_unit(std_out) 1194 1195! Calculate atomic energy from PAW matrix elements 1196 PAW%tkin=0; PAW%tion=0; PAW%tvale=0;PAW%txc=0;PAW%Ea=0 1197 PAW%Etotal=0;PAW%Eaxc=0;PAW%den=0; PAW%tden=0 1198 norbit=PAW%TOCCWFN%norbit 1199 DO io=1,norbit 1200 if (.NOT.PAW%TOCCWFN%iscore(io)) then 1201 l=PAW%TOCCWFN%l(io) ; occ=PAW%TOCCWFN%occ(io) 1202 CALL kinetic(Grid,PAW%TOCCWFN%wfn(:,io),l,x) 1203 PAW%tkin=PAW%tkin+occ*x 1204 PAW%den=PAW%den+occ*(PAW%OCCWFN%wfn(:,io))**2 1205 PAW%tden=PAW%tden+occ*(PAW%TOCCWFN%wfn(:,io))**2 1206 endif 1207 ENDDO 1208 write(std_out,*) 'smooth kinetic ', PAW%tkin 1209 1210 arg=PAW%den+FC%coreden-PAW%tden-PAW%tcore 1211 x=-Pot%nz+integrator(Grid,arg,1,irc) 1212 write(std_out,*) 'q00 for atom ', x 1213 arg=PAW%tden+PAW%tcore+x*PAW%hatden 1214 write(std_out,*) ' Total charge check ', integrator(Grid,arg) 1215 call poisson(Grid,x,arg,dum,y) 1216 write(std_out,*) ' Smooth coulomb ', y 1217 PAW%tvale=PAW%tkin+y 1218 arg=PAW%vloc*PAW%tden 1219 PAW%tion=integrator(Grid,arg,1,irc) 1220 write(std_out,*) ' Vloc energy ', PAW%tion 1221 PAW%tvale=PAW%tvale+PAW%tion 1222 1223 IF (TRIM(PAW%exctype)=="HF".or.TRIM(PAW%exctype)=="EXX".or.& 1224& TRIM(PAW%exctype)=="EXXKLI") THEN 1225 !CALL Get_Energy_EXX(Grid,PAW%TOCCWFN,x) ! not correct need moment 1226 CALL Get_Energy_EXX_pseudo(Grid,PAW,x) 1227 write(std_out,*) 'Warning: does not include core contributions' 1228 ELSE 1229 arg=PAW%tden+PAW%tcore 1230 CALL exch(Grid,arg,dum,y,x) 1231 ENDIF 1232 write(std_out,*) 'Smooth exchange-correlation contribution ', x 1233 PAW%txc=x ; PAW%tvale=PAW%tvale+PAW%txc 1234 1235! one center terms 1236 arg=PAW%den-PAW%tden 1237 x=integrator(Grid,arg,1,irc) 1238 write(std_out,*) 'valence Q', x 1239 PAW%Ea=-PAW%Eaion-x*PAW%Eaionhat 1240 1241 wij=0 1242 do io=1,norbit 1243 l=PAW%TOCCWFN%l(io); occ=PAW%TOCCWFN%occ(io) 1244 if (occ>1.d-8.and..NOT.PAW%TOCCWFN%iscore(io)) & 1245& call calcwij(Grid,PAW,l,occ,PAW%TOCCWFN%wfn(:,io),wij) 1246 enddo 1247 1248 1249 do ib=1,nbase 1250 do jb=1,nbase 1251 PAW%wij(lb,jb)=wij(ib,jb) 1252 enddo 1253 enddo 1254 DO ib=1,nbase 1255 do jb=1,nbase 1256 if (PAW%l(ib)==PAW%l(jb)) then 1257 PAW%Ea=PAW%Ea+PAW%wij(ib,jb)*(PAW%Kij(ib,jb)+PAW%VFij(ib,jb)) 1258 x=0 1259 do kb=1,nbase 1260 do lb=1,nbase 1261 if (PAW%l(kb)==PAW%l(lb)) then 1262 x=x+PAW%wij(kb,lb)*PAW%DR(ib,jb,kb,lb,1) 1263 endif 1264 enddo 1265 enddo 1266 PAW%Ea=PAW%Ea+0.5d0*x*PAW%wij(ib,jb) 1267 endif 1268 enddo 1269 enddo 1270 1271 write(std_out,*) 'Ea up to exchange-correlation terms ', PAW%Ea 1272 1273 IF (TRIM(PAW%exctype)=="HF".or.TRIM(PAW%exctype)=="EXX".or.& 1274& TRIM(PAW%exctype)=="EXXKLI") THEN 1275 CALL Get_Energy_EXX_pseudo_one(Grid,PAW,x) 1276 write(std_out,*) 'one-center HF exchange', x 1277 PAW%Ea=PAW%Ea+x; PAW%Eaxc=x 1278 ELSE 1279 arg=PAW%tden+PAW%tcore 1280 !CALL exch(Grid,arg(1:irc),dum(1:irc),y,x,fin=irc) 1281 CALL exch(Grid,arg,dum,y,x) 1282 arg=PAW%den+FC%coreden 1283 !CALL exch(Grid,arg(1:irc),dum(1:irc),y,z,fin=irc) 1284 CALL exch(Grid,arg,dum,y,z) 1285 write(std_out,*) ' one center xc ', z,x,z-x 1286 PAW%Ea=PAW%Ea+z-x; PAW%Eaxc=z-x 1287 ENDIF 1288 1289 PAW%Etotal=PAW%tvale+PAW%Ea 1290 write(std_out,*) 'Energy terms ', PAW%tvale, PAW%Ea, PAW%Etotal 1291 1292 PAW%mesh_size=PAW%irc+Grid%ishift 1293 PAW%coretailpoints=MAX(PAW%coretailpoints,PAW%mesh_size) 1294 1295 1296 DEALLOCATE(arg,dum,rh,trh,d,td,wij) 1297 1298 END SUBROUTINE SPMatrixElements 1299 1300!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1301! Calculate the core-valence matrix elements with correct (-) sign 1302!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1303 SUBROUTINE COREVAL_EXX(Grid,PAW) 1304 TYPE(GridInfo), INTENT(IN) :: Grid 1305 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 1306 INTEGER :: i,k, io,ib,ic,nbase,li,lj,lc,l 1307 REAL(8) :: accum,term,occ 1308 REAL(8), ALLOCATABLE :: f(:) 1309 1310 write(std_out,*) 'Entering Core-valence exchange subroutine ' 1311 write(std_out,*) '!!!WARNING -- further testing needed!!!' 1312 nbase=PAW%nbase 1313 k=0 ! core-valence terms 1314 do io=1,PAW%OCCWFN%norbit 1315 if (PAW%OCCWFN%iscore(io)) then 1316 li=PAW%OCCWFN%l(io) 1317 do ib=1,PAW%nbase 1318 do ic=1,PAW%nbase 1319 if (PAW%l(ib)==PAW%l(ic)) k=k+1 1320 enddo 1321 enddo 1322 endif 1323 enddo 1324 PAW%ncoreshell=k 1325 ALLOCATE (PAW%TXVC(nbase,nbase)); PAW%TXVC=0; 1326 1327 if(k>0) then 1328 ALLOCATE (PAW%DRVC(k,nbase,nbase),f(k)); PAW%DRVC=0 1329 !!!!WRITE(ifatompaw,'(" COREVAL_LIST ",i10)') k 1330 !!!!WRITE(ifatompaw,'(" COREVAL_R ")') 1331 i=0;f=0 1332 do io=1,PAW%OCCWFN%norbit 1333 if (PAW%OCCWFN%iscore(io)) then 1334 i=i+1 1335 lc=PAW%OCCWFN%l(io) 1336 occ=PAW%OCCWFN%occ(io) 1337 do ib=1,PAW%nbase 1338 lj=PAW%l(ib) 1339 do ic=1,PAW%nbase 1340 if (PAW%l(ib)==PAW%l(ic)) then 1341 f(i)=0 1342 do l=abs(lc-lj),(lc+lj),2 1343 call EXXwgt(1.d0,1.d0,1,lc,2,lj,l,accum) 1344 call CondonShortley(Grid,l,PAW%OCCWFN%wfn(:,io), & 1345& PAW%ophi(:,ib),PAW%OCCWFN%wfn(:,io), & 1346& PAW%ophi(:,ic),term) 1347 f(i)=f(i)+2*accum*term !EXXwgt returns 1/2*3J 1348 write(std_out,*) 'core-val CondonShortley',& 1349& i,li,lj,l,2*accum,term 1350 enddo 1351 !!!!!WRITE(ifatompaw,'(3i10,1p,1e25.17)') i, ib,ic,f(i) 1352 PAW%DRVC(i,ib,ic)=f(i) 1353 PAW%TXVC(ib,ic)=PAW%TXVC(ib,ic)-0.5d0*occ*f(i) 1354 endif 1355 enddo !ic 1356 enddo !ib 1357 endif 1358 enddo !norbit 1359 1360 Write(STD_OUT,*) 'Core - valence accumlated terms; output will be symmetrized' 1361 do ib=1,PAW%nbase 1362 do ic=ib,PAW%nbase 1363 write(std_out,'(2i5,2x,1p,2e15.7)') ib,ic,PAW%TXVC(ib,ic), & 1364 PAW%TXVC(ib,ic)-PAW%TXVC(ic,ib) 1365 PAW%TXVC(ib,ic)=0.5d0*(PAW%TXVC(ib,ic)+PAW%TXVC(ic,ib)) 1366 PAW%TXVC(ic,ib)=PAW%TXVC(ib,ic) 1367 enddo 1368 enddo 1369 1370 DEALLOCATE(f) 1371 endif 1372 1373 END SUBROUTINE COREVAL_EXX 1374 1375 SUBROUTINE CORECORE_EXX(Grid,PAW) 1376 1377!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1378! calculate core-core exchange interaction from all-electron core 1379! (assumed to be contained in augmentation radius) 1380! Added 8/9/2014 NAWH 1381!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1382 TYPE(GridInfo), INTENT(IN) :: Grid 1383 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 1384 INTEGER :: i,k, io,jo,l,lci,lcj 1385 REAL(8) :: accum,term,occi,occj,x 1386 1387 write(std_out,*) 'Entering Core-Core exchange subroutine ' 1388 write(std_out,*) '!!!WARNING -- further testing needed!!!' 1389 PAW%XCORECORE=0.d0;x=0.d0 1390 do io=1,PAW%OCCWFN%norbit 1391 if (PAW%OCCWFN%iscore(io)) then 1392 lci=PAW%OCCWFN%l(io) 1393 occi=PAW%OCCWFN%occ(io) 1394 do jo=1,PAW%OCCWFN%norbit 1395 if (PAW%OCCWFN%iscore(jo)) then 1396 lcj=PAW%OCCWFN%l(jo) 1397 occj=PAW%OCCWFN%occ(jo) 1398 do l=abs(lci-lcj),(lci+lcj),2 1399 call EXXwgt(occi,occj,io,lci,jo,lcj,l,accum) 1400 call CondonShortley(Grid,l,PAW%OCCWFN%wfn(:,io), & 1401& PAW%OCCWFN%wfn(:,jo),PAW%OCCWFN%wfn(:,io), & 1402& PAW%OCCWFN%wfn(:,jo),term) 1403 x=x-0.5d0*accum*term 1404 write(std_out,*) io,jo,l,accum,term,x 1405 enddo !l 1406 endif !jo core 1407 enddo !jo 1408 endif !io core 1409 enddo !io 1410 1411 write(std_out,*) 'Core-core exchange calculated: ', x 1412 PAW%XCORECORE=x 1413 END SUBROUTINE CORECORE_EXX 1414 1415END MODULE pseudo_sub 1416