1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2! The subroutines in this module were mainly written by Xiao Xu as 3! part of his Ph. D. work completed in 2011. They are not currently 4! used for typical runs of atompaw. The subroutine and function names 5! included are as follows: 6! EXXKLI_pseudoVx, Calc_tdexdphi_io, EXX_pseudoVx, InitialVx_pseudo, 7! Init_trvx, CompensationRHO, tVXsub0, tVXsub1, FindVX, multisolv, 8! Calc_w, Calc_u 9!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 10 11#if defined HAVE_CONFIG_H 12#include "config.h" 13#endif 14 15Module exx_pseudo 16 17 USE io_tools 18 USE anderson_driver 19 USE atomdata 20 USE exx_mod 21 USE fock 22 USE globalmath 23 USE gridmod 24 USE pseudodata 25 USE pseudo_sub 26 27 IMPLICIT NONE 28 29 INTEGER, PRIVATE :: S_n,S_nocc,S_nbase,S_irc,S_constr 30 INTEGER, ALLOCATABLE, TARGET, PRIVATE :: S_l(:),S_bmap(:),S_omap(:) 31 REAL(8), ALLOCATABLE, TARGET, PRIVATE :: S_gvv(:,:),S_tgvv(:,:) 32 REAL(8), ALLOCATABLE, TARGET, PRIVATE :: S_W(:,:),S_X(:,:),S_U(:),S_E(:) 33 REAL(8), ALLOCATABLE, TARGET, PRIVATE :: S_Vx(:,:),S_shift(:) 34 REAL(8), ALLOCATABLE, TARGET, PRIVATE :: S_Vxvale(:),S_tVxvale(:) 35 REAL(8), ALLOCATABLE, TARGET, PRIVATE :: S_dExdtphi(:,:),S_UOphi(:,:) 36 37 TYPE(Anderson_context), PRIVATE :: AC 38 TYPE(GridInfo), POINTER, PRIVATE :: Gridwk 39 TYPE(PseudoInfo), POINTER, PRIVATE :: PAW 40 !!!! Note: S_omap(:) maps occupied valence states to atomdata indices 41 !!!! S_bmap(:) maps occupied valence states to PAW basis indices 42 43 44 CONTAINS 45 46 SUBROUTINE EXXKLI_pseudoVx(Grid,PAW,rtvx) 47 TYPE(GridInfo), INTENT(IN) :: Grid 48 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 49 REAL(8), INTENT(INOUT) :: rtvx(:) 50 51 INTEGER :: i,j,k,l,n,io,jo 52 REAL(8), ALLOCATABLE :: dum1(:),dum2(:),dum3(:),den(:),tden(:) 53 REAL(8) :: occ 54 55 n=Grid%n 56 Allocate(dum1(n),dum2(n),dum3(n),den(n),tden(n)) 57 58 rtvx=0;dum1=0;dum2=0;den=0;tden=0 59 do io=1,PAW%TOCCWFN%norbit 60 occ=PAW%TOCCWFN%occ(io) 61 den=den+occ*PAW%OCCwfn%wfn(:,io)**2 62 tden=tden+occ*PAW%tOCCwfn%wfn(:,io)**2 63 call Calc_dexdphi_io(Grid,PAW%OCCwfn,io,dum3);PAW%OCCwfn%X(:,io)=dum3 64 ! not really needed -- just for checking 65 call Calc_tdexdphi_io(Grid,PAW,io,dum3) 66 PAW%tOCCwfn%X(:,io)=dum3 67 dum1=dum1+occ*(PAW%OCCwfn%X(:,io)*PAW%OCCwfn%wfn(:,io) & 68& + (PAW%OCCwfn%wfn(:,io)**2)*HSZ%U(io)*Grid%r) 69 dum2=dum2+occ*(PAW%tOCCwfn%X(:,io)*PAW%tOCCwfn%wfn(:,io) & 70& + (PAW%tOCCwfn%wfn(:,io)**2)*HSZ%U(io)*Grid%r) 71 enddo 72 73 rtvx(1)=0; dum3=0 74 do i=2,n 75 if (den(i)>machine_zero) then 76 dum3(i)=dum1(i)/den(i) 77 else 78 dum3(i)=0 79 endif 80 if (tden(i)>machine_zero) then 81 rtvx(i)=dum2(i)/tden(i) 82 else 83 rtvx(i)=0 84 endif 85 enddo 86 87 OPEN(1001,file='checkpseudovx',form='formatted') 88 do i=1,n 89 write(1001,'(1p,16e15.7)') Grid%r(i),dum1(i),dum2(i),dum3(i),& 90& rtvx(i),den(i),tden(i) 91 enddo 92 CLOSE(1001) 93 94 OPEN(1001,file='checkagain',form='formatted') 95 do i=1,n 96 write(1001,'(1p,50e15.7)') Grid%r(i),(PAW%OCCWFN%wfn(i,io),& 97& PAW%tOCCWFN%wfn(i,io),io=1,PAW%TOCCWFN%norbit) 98 enddo 99 CLOSE(1001) 100 101 102 Deallocate(dum1,dum2,dum3,den,tden) 103 END SUBROUTINE EXXKLI_pseudoVx 104 105!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 106 ! Calc_texpot_io(Grid,Orbit,io,res) 107! Determine exchange potential contribution from single pseudo orbital 108! Appropriate for reference configuration 109!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 110 SUBROUTINE Calc_tdexdphi_io(Grid,PAW,io,res) 111 TYPE(GridInfo), INTENT(IN):: Grid 112 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 113 INTEGER, INTENT(IN) :: io 114 REAL(8), INTENT(OUT) :: res(:) 115 116 REAL(8), POINTER :: r(:) 117 TYPE(OrbitInfo), POINTER :: Orbit, tOrbit 118 REAL(8), ALLOCATABLE :: wfp(:),vl(:),arg(:),dum(:),dum1(:) 119 INTEGER :: i,j,k,n,m,l,ll,li,lj,lmin,lmax,jo,ok 120 INTEGER :: nu,nup,ni,nj 121 REAL(8) :: occ,term,term1,wgt,occi,occj,rc,uvjo 122 REAL(8), PARAMETER :: threshold=1.d-8 123 124 write(std_out,*) 'Calc_tecpot_io ',io; call flush_unit(std_out) 125 n=Grid%n 126 Orbit=>PAW%OCCwfn 127 tOrbit=>PAW%tOCCwfn 128 ALLOCATE(wfp(n),vl(n),arg(n),dum(n),dum1(n),stat=ok) 129 IF (ok/=0) THEN 130 WRITE(STD_OUT,*) 'calc_expot_io allocation error:', n,Orbit%norbit,ok 131 STOP 132 ENDIF 133 134 r=>Grid%r 135 136 res=0 137 occ=orbit%occ(io) 138 if(occ<threshold) return 139 li=Orbit%l(io);ni=Orbit%np(io) 140 DO jo=1,Orbit%norbit 141 occj=Orbit%occ(jo);vl=0 142 lj=Orbit%l(jo);nj=Orbit%np(jo) 143 IF(occj>threshold) THEN 144 lmax=li+lj 145 lmin=ABS(li-lj) 146 wfp(1:n)=tOrbit%wfn(1:n,io)*tOrbit%wfn(1:n,jo) 147 arg(1:n)=Orbit%wfn(1:n,io)*Orbit%wfn(1:n,jo)-wfp(1:n) 148 DO ll=lmin,lmax,2 149 CALL EXXwgt(occ,occj,io,li,jo,lj,ll,wgt) 150 IF (wgt>threshold) THEN 151 wgt=wgt*(2*ll+1) ! because of apoisson convention 152 dum=(r**ll)*arg;term1=integrator(Grid,dum) 153 dum1=wfp+term1*PAW%g(:,ll+1) 154 CALL apoisson(Grid,ll,n,dum1,dum) 155 vl(1:n)=vl(1:n)-wgt*dum(1:n)*tOrbit%wfn(1:n,jo)/occ 156 ENDIF 157 ENDDO 158 ENDIF 159 res=res+vl 160 ENDDO 161 162 163 DEALLOCATE(dum,wfp,arg,vl) 164 END SUBROUTINE Calc_tdexdphi_io 165 166!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 167! Assumes tildeg(r) = \sum_n=1^N C_n r^{lead+n+l} 168 169 SUBROUTINE EXX_pseudoVx(Grid,Orbit,PAWin,trvx) 170 Type(Gridinfo), TARGET, INTENT(IN) :: Grid 171 Type(Orbitinfo), INTENT(IN) :: Orbit 172 Type(PseudoInfo), TARGET, INTENT(IN) :: PAWin 173 REAL(8), INTENT(OUT) :: trvx(:) 174 175 LOGICAL :: success 176 INTEGER :: i,j,k,l,m,n,nbase,io,ip,iq,ib,ic,nocc,kmax,kmin,lp,lq,irc,ind 177 INTEGER :: ii,jj,kk,loop,outer,p 178 INTEGER :: nequations,nterms 179 INTEGER, POINTER :: map(:) 180 INTEGER, PARAMETER :: nshift=4,nmatch=3,lead=4 181 REAL(8) :: occp,occq,wgt,xxx,yyy,r1 182 REAL(8), ALLOCATABLE :: trvx0(:),d(:),d1(:),d2(:),d3(:),v(:),arg(:),F2(:,:) 183 REAL(8), ALLOCATABLE :: F(:),rho(:) 184 INTEGER, PARAMETER :: mxloop=100 185 REAL(8), PARAMETER :: threshold=1.d-8,CondNo=1.d12 186 REAL(8) :: svdcut=1.d-8 187 LOGICAL :: parameterfileexists 188 CHARACTER(120) :: inputline 189 190 Gridwk=>Grid 191 PAW=>PAWin 192 n=Grid%n; S_n=n; irc=PAW%irc; S_irc=irc 193 194 Inquire(file='in_parameters',exist=parameterfileexists) 195 196 If(parameterfileexists) then 197 Open(1002,file='in_parameters') 198 do 199 READ(1002,'(a)') inputline 200 Call Uppercase (inputline) 201 If (TRIM(inputline)=='<SVDCUTOFF>') exit 202 enddo 203 READ(1002,*) svdcut 204 WRITE(STD_OUT,*) 'Modified SVD Cutoff ', svdcut 205 close(1002) 206 endif 207 208 trvx=0 209 210! nocc=0; nbase=PAW%nbase 211! do io=1,nbase 212! if (PAW%occ(io)>threshold) nocc=nocc+1 213! enddo 214! 215! ALLOCATE(S_bmap(nocc),S_omap(nocc),S_E(nocc),S_l(nocc),S_U(nocc),& 216!& S_W(n,nocc),S_X(nbase,nocc),S_Vx(nbase,nocc), S_shift(n), & 217!& S_dExdtphi(n,nocc),S_UOphi(n,nocc),S_Vxvale(n),S_tVxvale(n)) 218! ALLOCATE(trvx0(n),d(n),d1(n),d2(n),d3(n),v(n),arg(irc),F2(n,nocc),& 219!& F(n),rho(n)) 220! 221! map=>S_bmap 222! S_nocc=nocc; S_nbase=nbase; 223! S_bmap=0; S_omap=0 224! 225! io=0 226! do ib=1,nbase 227! if (PAW%occ(ib)>threshold) then 228! io=io+1 229! map(io)=ib 230! endif 231! enddo 232! 233! Do ip=1,nocc 234! S_E(ip)=PAW%eig(map(ip)) 235! S_l(ip)=PAW%l(map(ip)) 236! i=0 237! do io=1,Orbit%norbit 238! if (ABS(Orbit%eig(io)-S_E(ip))<1.d-8) then 239! i=io; S_omap(ip)=io; write(std_out,*) 'ip io map', ip,io ; exit 240! endif 241! enddo 242! if (i==0) then 243! write(std_out,*) 'Error: No match found for ',ip,S_E(ip) 244! stop 245! endif 246! Enddo 247! 248! allocate(S_gvv(n,nocc),S_tgvv(n,nocc)) 249! do ip=1,nocc 250! S_gvv(:,ip)=HSZ%psivale(:,S_omap(ip)) 251! S_tgvv(:,ip)=HSZ%psivale(:,S_omap(ip)) 252! enddo 253! 254! 255! S_W=0;S_X=0 256! F2=0;F=0;rho=0; 257! Do ip=1,nocc 258! occp=PAW%occ(map(ip)) 259! lp=PAW%l(map(ip)) 260! Do iq=1,nocc 261! occq=PAW%occ(map(iq)) 262! lq=PAW%l(map(iq)) 263! kmax=lp+lq 264! kmin=ABS(lp-lq) 265! d1=PAW%otphi(:,map(ip))*PAW%otphi(:,map(iq)) 266! d3=PAW%ophi(:,map(ip))*PAW%ophi(:,map(iq)) 267! Do k=kmin,kmax,2 268! CALL EXXwgt(occp,occq,map(ip),lp,map(iq),lq,k,wgt) 269! IF (wgt>threshold) THEN 270! wgt=(wgt*(2*k+1))/occp ! because of apoisson convention 271! CALL CompensationRHO(Grid,PAW,k,map(ip),map(iq),d2) 272! d=d1+d2 273! CALL apoisson(Grid,k,n,d,trvx0) 274! S_W(:,ip)=S_W(:,ip)-wgt*trvx0(:)*PAW%otphi(:,map(iq)) 275! CALL apoisson(Grid,k,n,d3,v) 276! F2(:,ip)=F2(:,ip)-wgt*v(:)*PAW%ophi(:,map(iq)) 277! do ib=1,nbase 278! if (PAW%l(ib)==lp) then 279! d=PAW%ophi(:,ib)*PAW%ophi(:,map(iq))*v - & 280!& PAW%otphi(:,ib)*PAW%otphi(:,map(iq))*trvx0 281! d(1)=0; d(2:n)=d(2:n)/Grid%r(2:n) 282! S_X(ib,ip)=S_X(ib,ip)-wgt*integrator(Grid,d,1,irc) 283! write(std_out,*) 'Check ',ib,ip,integrator(Grid,d,1,irc),& 284!& integrator(Grid,d) 285! endif 286! enddo 287! ENDIF 288! ENDDO !k 289! ENDDO !iq 290! 291! 292! 293! ! (re)calculate U factor 294! d=(F2(:,ip)-HSZ%rVxvale(:)*PAW%ophi(:,map(ip)))*PAW%ophi(:,map(ip)) 295! d(1)=0;d(2:n)=d(2:n)/Grid%r(2:n) 296! S_U(ip)=integrator(Grid,d) 297! write(std_out,*) 'Checking U ', ip,map(ip),integrator(Grid,d),& 298!& HSZ%Uvale(S_omap(ip)) 299! d=d-HSZ%Uvale(S_omap(ip))*PAW%ophi(:,map(ip))**2 300! write(std_out,*) 'FC rhs',ip,integrator(Grid,d),integrator(Grid,d,1,irc) 301! d=S_W(:,ip)*PAW%otphi(:,map(ip)) 302! d(1)=0;d(2:n)=d(2:n)/Grid%r(2:n) 303! write(std_out,*) 'Checking dtEx ', ip,map(ip),& 304!& integrator(Grid,d)+S_X(map(ip),ip),S_X(map(ip),ip) 305! d=F2(:,ip)*PAW%ophi(:,map(ip)) 306! d(1)=0;d(2:n)=d(2:n)/Grid%r(2:n) 307! write(std_out,*) 'Checking dEx ', ip,map(ip),integrator(Grid,d) 308! d=S_W(:,ip)*PAW%otphi(:,map(ip)) 309! d(1)=0;d(2:n)=d(2:n)/Grid%r(2:n) 310! F(:)=F(:)+occp*d(:) 311! rho(:)=rho(:)+occp*PAW%otphi(:,map(ip))**2 312! do iq=1,nocc 313! if (iq==ip) then 314! F(:)=F(:)-occp*HSZ%Uvale(S_omap(ip))*PAW%otphi(:,map(ip))**2 315! else 316! !if (S_l(iq)==S_l(ip)) then 317! !F(:)=F(:)-occp*HSZ%LMBDvale(S_omap(ip),S_omap(iq))& 318!& ! *PAW%otphi(:,map(ip))*PAW%otphi(:,map(iq)) 319! !endif 320! endif 321! enddo 322! ENDDO 323! ! accumulate dExdtphi 324! S_dExdtphi=0 ; S_UOphi=0 325! do ip=1,nocc 326! S_dExdtphi(2:n,ip)=S_W(2:n,ip)/Grid%r(2:n) 327! S_UOphi(:,ip)=PAW%otphi(:,map(ip)) 328! do ib=1,nbase 329! if (PAW%l(ib)==S_l(ip)) then 330! S_dExdtphi(:,ip)=S_dExdtphi(:,ip)+S_X(ib,ip)*PAW%otp(:,ib) 331! S_UOphi(:,ip)=S_UOphi(:,ip)+PAW%Oij(ib,map(ip))*PAW%otp(:,ib) 332! endif 333! enddo 334! S_UOphi(:,ip)=HSZ%Uvale(S_omap(ip))*S_UOphi(:,ip) 335! enddo 336! 337! ! Calculated matrix element <\phi(ib)|VxVale|\phi(ip)> 338! S_Vx=0 339! do ip=1,nocc 340! lp=PAW%l(S_bmap(ip)) ; 341! do ib=1,nbase 342! if (lp==PAW%l(ib)) THEN 343! d=HSZ%rVxVale(:)*PAW%ophi(:,ib)*PAW%ophi(:,S_bmap(ip)) 344! d(1)=0; d(2:n)=d(2:n)/Gridwk%r(2:n) 345! S_Vx(ib,ip)=integrator(Gridwk,d,1,irc) 346! endif 347! write(std_out,'(" vx ",2i5,1p,1e15.7)') ib,ip,S_Vx(ib,ip) 348! enddo 349! enddo 350! 351! ! Note that S_W is r*(desired function) 352! 353! Open(1001,file='testF1',form='formatted') 354! do i=1,n 355! write(1001,'(1p,50e15.7)') Grid%r(i),HSZ%rVxvale(i),& 356!& (S_W(i,ip),F2(i,ip),ip=1,nocc) 357! enddo 358! close(1001) 359! 360! Do ip=1,nocc 361! do ib=1,nbase 362! write(std_out,*) 'X ', ip,ib,S_X(ib,ip) 363! enddo 364! enddo 365! 366! call Init_trvx(Grid,irc,HSZ%rVxvale,trvx) 367! S_Vxvale=HSZ%rVxvale; S_tVxvale=trvx 368! 369! Open(1001,file='trvx0',form='formatted') 370! do i=1,n 371! write(1001,'(1p,50e15.7)') Grid%r(i),trvx(i),HSZ%rVxvale(i),& 372!& HSZ%rVxcore(i) 373! enddo 374! close(1001) 375! 376! arg(1:irc)=trvx(1:irc) 377! CALL InitAnderson_dr(AC,6,5,irc,0.1d0,1.d3,100,& 378!& 1.d-8,1.d-16,.TRUE.) 379! CALL DoAndersonMix(AC,arg,xxx,tVXsub1,success) 380! trvx(1:irc)=arg(1:irc) 381! WRITE(STD_OUT,*) 'Finished trvx iter # ',AC%CurIter,AC%res 382! CALL FreeAnderson(AC) 383! 384! Open(1001,file='trvx',form='formatted') 385! do i=1,n 386! write(1001,'(1p,50e15.7)') Grid%r(i),trvx(i),HSZ%rVxvale(i),& 387!& HSZ%rVxcore(i) 388! enddo 389! close(1001) 390! 391! Deallocate(trvx0,d,d1,d2,d3,v,F2,arg,rho,F) 392! 393 END SUBROUTINE EXX_pseudoVx 394 395!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 396! ib is the PAW basis index of the outer 397! most orbital for initializing tVx 398!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 399 SUBROUTINE InitialVx_pseudo(Grid,ib,trvx0) 400 TYPE(GridInfo), INTENT(IN) :: Grid 401 INTEGER, INTENT(IN) :: ib 402 REAL(8), INTENT(OUT) :: trvx0(:) 403 404 INTEGER :: lp,kmax,kmin,k,n 405 REAL(8) :: occp,wgt 406 REAL(8), ALLOCATABLE :: d(:),d1(:),d2(:) 407 REAL(8), PARAMETER :: threshold=1.d-8 408 409 n=Grid%n 410 allocate(d(n),d1(n),d2(n)) 411 412 trvx0=0 413 occp=PAW%occ(ib) 414 lp=PAW%l(ib) 415 kmax=2*lp 416 kmin=0 417 d1=(PAW%otphi(:,ib)**2) 418 Do k=kmin,kmax,2 419 CALL EXXwgt(occp,occp,ib,lp,ib,lp,k,wgt) 420 IF (wgt>threshold) THEN 421 wgt=(wgt*(2*k+1))/occp ! because of apoisson convention 422 CALL CompensationRHO(Grid,PAW,k,ib,ib,d) 423 d=d1+d 424 CALL apoisson(Grid,k,n,d,d2) 425 trvx0=trvx0-wgt*d2 426 ENDIF 427 ENDDO !k 428 429 deallocate(d,d1,d2) 430 END SUBROUTINE InitialVx_pseudo 431 432 SUBROUTINE Init_trvx(Grid,p,trvxin,trvxout) 433 Type(GridInfo),INTENT(IN) :: Grid 434 INTEGER, INTENT(IN) :: p ! smooth potential for r<r(irc+1) 435 REAL(8), INTENT(IN) ::trvxin(:) 436 REAL(8), INTENT(OUT) ::trvxout(:) 437 438 INTEGER :: n,i 439 REAL(8) :: r1u1,r2u2,r3u3,d0rc, d1rc,d2rc,rc,x 440 REAL(8) :: u1,u2,u3 441 REAL(8), ALLOCATABLE :: d1(:),d2(:) 442 443 n=Grid%n 444 ALLOCATE(d1(n),d2(n)) 445 446 CALL derivative(Grid,trvxin,d1) 447 CALL derivative(Grid,d1,d2) 448 449 d0rc=trvxin(p) 450 d1rc=d1(p) 451 d2rc=d2(p) 452 rc=Grid%r(p) 453 454 r2u2 = (-d0rc + rc*d1rc -(1.0d0/3.0d0)*rc*rc*d2rc)*3 455 r3u3 = (rc*d1rc -d0rc - r2u2)*(1.0d0/2.0d0) 456 r1u1 = d0rc - r2u2 - r3u3 457 458 u1 = r1u1/rc; 459 u2 = r2u2/(rc*rc) 460 u3 = r3u3/(rc**3) 461 462 trvxout=trvxin 463 464 do i=1,p 465 x = Grid%r(i) 466 trvxout(i) = x*u1 + x*x*u2 + x*x*x*u3 467 enddo 468 469 open(unit=2001,file='rvxsmooth.txt') 470 n=Grid%n 471 do i=1,n 472 write(2001,'(1p,6e15.7)') Grid%r(i), trvxin(i),trvxout(i) 473 enddo 474 close(2001) 475 476 END SUBROUTINE Init_trvx 477 478 SUBROUTINE CompensationRHO(Grid,PAW,k,alpha,beta,rho) 479 TYPE(GridInfo), INTENT(IN) :: Grid 480 TYPE(PseudoInfo), INTENT(IN) :: PAW 481 INTEGER,INTENT(IN) :: k,alpha,beta 482 REAL(8),INTENT(INOUT) :: rho(:) 483 484 REAL(8) , ALLOCATABLE :: Temp1(:),Temp2(:) 485 INTEGER :: n,irc 486 REAL(8) :: n1 487 REAL(8) ,POINTER :: r(:) 488 489 n=Grid%n 490 ALLOCATE(Temp1(n),Temp2(n)) 491 r=>Grid%r 492 irc = PAW%irc 493 494 Temp1 = PAW%ophi(:,alpha)*PAW%ophi(:,beta) - & 495& PAW%otphi(:,alpha)*PAW%otphi(:,beta) 496 Temp1 = Temp1 * (r(:)**k) 497 n1 = integrator(Grid,Temp1,1,irc) 498 rho(:)=(r(:)**(k+2))*PAW%hatshape(:) 499 Temp1=(r(:)**(k))*rho(:) 500 rho(:) =n1*rho(:)/integrator(Grid,Temp1) 501 502 END SUBROUTINE CompensationRHO 503 504 505!!!!!!! 506! tVXsub0 -- version without projector constraints 507!!!!!!! 508 SUBROUTINE tVXsub0(w,en,grad,err,success,update) 509 REAL(8), INTENT(INOUT) :: w(:) 510 REAL(8), INTENT(OUT) :: en 511 REAL(8), INTENT(OUT) :: grad(:) 512 REAL(8), INTENT(OUT) :: err 513 LOGICAL, INTENT(OUT) :: success 514 LOGICAL, INTENT(IN) :: update 515 516 INTEGER :: i,j,k,n,io,l,iter,last,irc,ip,nocc,nbase,many,ib,ic,jo,p 517 REAL(8) :: x,y,energy,occ,boundaryv 518 REAL(8), POINTER :: rv(:) 519 REAL(8),ALLOCATABLE :: shift(:),dvx(:,:),tg(:,:),tu(:),rhs(:),tw(:) 520 INTEGER:: fcount=0 521 REAL(8), PARAMETER :: tol=1.d-10, wgt0=0.01d0 522 CHARACTER(4) :: stuff 523 524 success=.TRUE. 525 irc=S_irc; nocc=S_nocc; n=S_n; nbase=PAW%nbase; p=irc+1 526 527 rv=>PAW%rVeff 528 529 allocate(shift(n),tg(n,nocc),tu(n),rhs(n),tw(n),& 530& dvx(nbase,nocc)) 531 shift=0;grad=0 532 533 ! Matrix element of \tilde{V}_x = w (input) 534 do ip=1,nocc 535 dvx(:,ip)=0; l=PAW%l(S_bmap(ip)) ; 536 do ib=1,nbase 537 if (l==PAW%l(ib)) THEN 538 tu(1:irc)=w(1:irc)*PAW%otphi(1:irc,ib)*PAW%otphi(1:irc,S_bmap(ip)) 539 tu(1)=0; tu(2:irc)=tu(2:irc)/Gridwk%r(2:irc) 540 dvx(ib,ip)=S_Vx(ib,ip)-integrator(Gridwk,tu,1,irc) 541 endif 542 write(std_out,'(" dvx ",2i5,1p,1e15.7)') ib,ip,dvx(ib,ip) 543 enddo 544 enddo 545 call flush_unit(std_out) 546 tw=S_tVxvale 547 tw(1:irc)=w(1:irc) 548 shift=0 ; 549 do ip=1,nocc 550 occ=PAW%occ(S_bmap(ip)) 551 energy=S_E(ip) 552 l=S_l(ip) 553 rhs=0; 554 rhs(2:n)=S_dExdtphi(2:n,ip) & 555& -tw(2:n)*PAW%otphi(2:n,S_bmap(ip))/Gridwk%r(2:n) 556 rhs(2:n)=rhs(2:n)-S_UOphi(2:n,ip) 557 558 do ib=1,nbase 559 if (l==PAW%l(ib)) then 560 rhs(:)=rhs(:)-dvx(ib,ip)*PAW%otp(:,ib) 561 endif 562 enddo 563 564 write(std_out,*)' Check smooth rhs ip ', ip, & 565& overlap(Gridwk,rhs,PAW%otphi(:,S_bmap(ip)),1,irc),& 566& overlap(Gridwk,rhs,PAW%otphi(:,S_bmap(ip))) 567 568 569 call inhomo_numerov_SVD_bv(Gridwk,l,irc+1,energy,tol,rv,rhs,S_tgvv(:,ip)) 570 571 ! check orthogonalities 572 573 do ib=1,nbase 574 if (l==PAW%l(ib)) then 575 write(std_out,*) ' ip ib <g|p> = ', ip,ib ,& 576& overlap(Gridwk,S_tgvv(:,ip),PAW%otp(:,ib)) 577 write(std_out,*) ' ip ib <g|tphi> = ', & 578& overlap(Gridwk,S_tgvv(:,ip),PAW%otphi(:,ib)) 579 endif 580 enddo 581 enddo !ip 582 583 shift=0;tu=0 584 do ip=1,nocc 585 occ=PAW%occ(S_bmap(ip)) 586 587 shift(:)=shift(:)+2*occ*S_tgvv(:,ip)*PAW%otphi(:,S_bmap(ip)) 588 enddo 589 590 grad=0 591 grad(2:irc)=-shift(2:irc)/Gridwk%r(2:irc) 592 err=DOT_PRODUCT(grad(1:irc),grad(1:irc)) 593 en=err 594 S_shift=shift 595 596 WRITE(STD_OUT,'("PAWiter",i5,1p,1e20.12,1p,2e15.7)')fcount,en 597 598 call mkname(fcount,stuff) 599 open(1001,file='pseudo.'//TRIM(stuff),form='formatted') 600 do i=1,n 601 write(1001,'(1p,15e15.7)') Gridwk%r(i),tw(i),shift(i),& 602& (S_tgvv(i,ip),ip=1,nocc),& 603& (S_tgvv(i,ip)*PAW%tphi(i,S_bmap(ip)),ip=1,nocc) 604 enddo 605 close(1001) 606 607 fcount=fcount+1 608 deallocate(shift,tg,tu,rhs,tw,dvx) 609 610 END SUBROUTINE tVXsub0 611 612!!!!!!! 613! tVXsub1 614!!!!!!! 615 SUBROUTINE tVXsub1(w,en,grad,err,success,update) 616 REAL(8), INTENT(INOUT) :: w(:) 617 REAL(8), INTENT(OUT) :: en 618 REAL(8), INTENT(OUT) :: grad(:) 619 REAL(8), INTENT(OUT) :: err 620 LOGICAL, INTENT(OUT) :: success 621 LOGICAL, INTENT(IN) :: update 622 623 INTEGER :: i,j,k,n,io,l,iter,last,irc,ip,nocc,nbase,many,ib,ic,jo,p 624 REAL(8) :: x,y,energy,occ,boundaryv 625 REAL(8), POINTER :: rv(:) 626 REAL(8),ALLOCATABLE :: shift(:),dvx(:,:),tg(:,:),tu(:),rhs(:,:),tw(:) 627 INTEGER:: fcount=0 628 REAL(8), PARAMETER :: tol=1.d-10, wgt0=0.01d0 629 CHARACTER(4) :: stuff 630 631 success=.TRUE. 632 irc=S_irc; nocc=S_nocc; n=S_n; nbase=PAW%nbase; p=irc+1 633 634 rv=>PAW%rVeff 635 636 allocate(shift(n),tg(n,nocc),tu(n),rhs(n,nbase+1),tw(n),& 637& dvx(nbase,nocc)) 638 shift=0;grad=0 639 640 ! Matrix element of \tilde{V}_x = w (input) 641 do ip=1,nocc 642 dvx(:,ip)=0; l=PAW%l(S_bmap(ip)) ; 643 do ib=1,nbase 644 if (l==PAW%l(ib)) THEN 645 tu(1:irc)=w(1:irc)*PAW%otphi(1:irc,ib)*PAW%otphi(1:irc,S_bmap(ip)) 646 tu(1)=0; tu(2:irc)=tu(2:irc)/Gridwk%r(2:irc) 647 dvx(ib,ip)=S_Vx(ib,ip)-integrator(Gridwk,tu,1,irc) 648 endif 649 write(std_out,'(" dvx ",2i5,1p,1e15.7)') ib,ip,dvx(ib,ip) 650 enddo 651 enddo 652 call flush_unit(std_out) 653 tw=S_tVxvale 654 tw(1:irc)=w(1:irc) 655 shift=0 ; 656 do ip=1,nocc 657 occ=PAW%occ(S_bmap(ip)) 658 energy=S_E(ip) 659 l=S_l(ip) 660 rhs=0; 661 rhs(2:n,1)=S_dExdtphi(2:n,ip) & 662& -tw(2:n)*PAW%otphi(2:n,S_bmap(ip))/Gridwk%r(2:n) 663 rhs(2:n,1)=rhs(2:n,1)-S_UOphi(2:n,ip) 664 665 many=1 666 do ib=1,nbase 667 if (l==PAW%l(ib)) then 668 rhs(:,1)=rhs(:,1)-dvx(ib,ip)*PAW%otp(:,ib) 669 many=many+1 670 rhs(:,many)=PAW%otp(:,ib) 671 endif 672 enddo 673 674 S_tgvv(:,ip)=0 675 If (many<=2) then 676 call inhomo_numerov_SVD_bv(Gridwk,l,n-1,energy,tol,rv,rhs(:,1),S_tgvv(:,ip)) 677 else 678 call multisolv(S_bmap(ip),l,many,energy,rv,rhs(:,1:many),S_tgvv(:,ip)) 679 endif 680 681 ! Remove homogeneous solution 682 x=0; y=0 683 write(std_out,*) 'Fix tail ', ip,S_bmap(ip) 684 do i=irc+1,n 685 x=x+PAW%otphi(i,S_bmap(ip))**2 686 y=y+(S_tgvv(i,ip)-S_gvv(i,ip))*PAW%otphi(i,S_bmap(ip)) 687 enddo 688 write(std_out,*) 'Adjust tail value ', ip ,y,x 689 S_tgvv(:,ip)=S_tgvv(:,ip)-(y/x)*PAW%otphi(:,S_bmap(ip)) 690 691 do ib=1,nbase 692 if (l==PAW%l(ib)) then 693 write(std_out,*) ' ip ib <g|p> = ', ip,ib ,& 694& overlap(Gridwk,S_tgvv(:,ip),PAW%otp(:,ib)) 695 write(std_out,*) ' ip ib <g|tphi> = ', & 696& overlap(Gridwk,S_tgvv(:,ip),PAW%otphi(:,ib)) 697 endif 698 enddo 699 enddo !ip 700 701 shift=0;tu=0 702 do ip=1,nocc 703 occ=PAW%occ(S_bmap(ip)) 704 shift(:)=shift(:)+2*occ*S_tgvv(:,ip)*PAW%otphi(:,S_bmap(ip)) 705 enddo 706 707 grad=0 708 grad(2:irc)=-shift(2:irc)/Gridwk%r(2:irc) 709 err=DOT_PRODUCT(grad(1:irc),grad(1:irc)) 710 en=err 711 S_shift=shift 712 713 WRITE(STD_OUT,'("PAWiter",i5,1p,1e20.12,1p,2e15.7)')fcount,en 714 715 call mkname(fcount,stuff) 716 open(1001,file='pseudo.'//TRIM(stuff),form='formatted') 717 do i=1,n 718 write(1001,'(1p,15e15.7)') Gridwk%r(i),tw(i),shift(i),& 719& (S_gvv(i,ip),S_tgvv(i,ip),PAW%otphi(i,S_bmap(ip)),ip=1,nocc) 720 enddo 721 close(1001) 722 723 fcount=fcount+1 724 deallocate(shift,tg,tu,rhs,tw,dvx) 725 726 END SUBROUTINE tVXsub1 727 728 SUBROUTINE FindVX(mxloop,n,threshold,arg,trvx) 729 Integer, INTENT(IN) :: mxloop,n 730 REAL(8), INTENT(IN) :: threshold 731 REAL(8), INTENT(INOUT) :: arg(:) 732 REAL(8), INTENT(INOUT) :: trvx(:) 733 734 LOGICAL :: success 735 REAL(8) :: xxx 736 INTEGER :: loop,i 737 INTEGER, parameter :: mxl=50 738 REAL(8), parameter :: mix=0.2d0 739 CHARACTER(4) :: nm 740 trvx=arg 741 Do loop=1,mxloop*10 742 CALL InitAnderson_dr(AC,6,2,n,0.001d0,1.d1,mxl,1.d-6,1.d-6,.true.) 743 CALL DoAndersonMix(AC,arg,xxx,tVXsub0,success) 744 write(std_out,*) 'Finished EXX_PseudoVx',AC%CurIter,AC%res 745 CALL FreeAnderson(AC) 746 747 trvx=arg 748 arg=0 749 arg(2:n)=-S_shift(2:n)/Gridwk%r(2:n) 750 xxx=DOT_PRODUCT(arg(1:n),arg(1:n)) 751 write(std_out,*) 'tVX loop ', loop,xxx 752 call mkname(loop,nm) 753 Open(1001,file='anal'//TRIM(nm),form='formatted') 754 do i=1,n 755 write(1001,'(1p,20e15.7)') Gridwk%r(i),trvx(i),S_shift(i),& 756& trvx(i)+mix*arg(i) 757 enddo 758 close(1001) 759 if (xxx<threshold) then 760 write(std_out,*) ' tVX loop converged with ', loop,xxx 761 exit 762 endif 763 trvx=trvx+mix*arg 764 arg=trvx 765 enddo 766 767 END SUBROUTINE FindVX 768 769 770 SUBROUTINE multisolv(in,l,mult,energy,rv,rhs,res) 771 INTEGER, INTENT(IN):: in,l,mult 772 REAL(8), INTENT(IN):: energy,rv(:),rhs(:,:) 773 REAL(8), INTENT(INOUT):: res(:) 774 775 integer :: i,io,jo,ib,many,n,nbase 776 integer, ALLOCATABLE :: map(:) 777 REAL(8), ALLOCATABLE :: tu(:),tw(:,:),M(:,:),MM(:,:),MMM(:,:),c(:),cc(:) 778 REAL(8), PARAMETER :: tol=1.d-10 779 780 n=S_n; res=0 781 nbase=PAW%nbase 782 allocate(map(nbase),M(nbase,nbase),MM(nbase,nbase),MMM(nbase,nbase),& 783& tw(n,nbase+1),tu(n),c(nbase),cc(nbase)) 784 785 many=0; map=0 786 do ib=1,nbase 787 if (l==PAW%l(ib)) then 788 many=many+1; map(many)=ib; 789 endif 790 enddo 791 792 if (mult-1/=many) then 793 write(std_out,*) 'Error in multisolv ', mult,many 794 stop 795 endif 796 797 call inhomo_numerov_SVD_bvm(Gridwk,l,n-1,mult,energy,tol,rv,rhs(:,1:mult),& 798& tw(:,1:mult)) 799 800 MMM=0; M=0 801 do io=1,many 802 do jo=1,many 803 tu(:)=PAW%otp(:,map(io))*tw(:,jo+1) 804 MMM(io,jo)=integrator(Gridwk,tu) 805 M(io,jo)=PAW%dij(map(io),map(jo))-& 806& energy*PAW%oij(map(io),map(jo)) 807 enddo 808 enddo 809 MM=0 810 do io=1,many 811 do jo=1,many 812 do ib=1,many 813 MM(io,jo)=MM(io,jo)+M(io,ib)*MMM(ib,jo) 814 enddo 815 enddo 816 MM(io,io)=MM(io,io)+1 817 enddo 818 819 c=0;cc=0 820 do io=1,many 821 tu(:)=PAW%otp(:,map(io))*tw(:,1) 822 c(io)=integrator(Gridwk,tu) 823 enddo 824 do io=1,many 825 do jo=1,many 826 cc(io)=cc(io)-M(io,jo)*c(jo) 827 enddo 828 enddo 829 830 MMM=MM 831 call linsol(MMM,cc,many,nbase,nbase,nbase) 832 833 write(std_out,*) ' linsol results ' 834 do io=1,many 835 write(std_out,*) io,cc(io) 836 enddo 837 838 res=tw(:,1) 839 do io=1,many 840 res(:)=res(:)+cc(io)*tw(:,io+1) 841 enddo 842 843 deallocate(map,M,MM,MMM,tw,tu,c,cc) 844 845 END SUBROUTINE multisolv 846 847 SUBROUTINE Calc_w(Grid,PAW,io,w) 848 TYPE(GridInfo), INTENT(IN) :: Grid 849 TYPE(PseudoInfo), INTENT(IN) :: PAW 850 INTEGER, INTENT(IN) :: io 851 REAL(8), INTENT(OUT) :: w(:,:) ! w(1:n,1:nbase) 852 853 INTEGER :: i,j,l,n,nbase 854 REAL(8) :: eig,zeroval,term 855 REAL(8),POINTER :: rv(:),r(:) 856 REAL(8) ,ALLOCATABLE :: rhs(:),temp(:) 857 858 eig=PAW%eig(io) 859 l=PAW%l(io) !<-- lp 860 rv=>PAW%rveff 861 n=Grid%n 862 r=>Grid%r 863 ALLOCATE(rhs(n),temp(n)) 864 rhs =0 865 w=0 866 temp=0 867 868 nbase=PAW%nbase 869 DO i=1,nbase 870 IF (l==PAW%l(i))THEN ! <-- lp = lk 871 rhs = PAW%otp(:,i) 872 CALL inhomo_bound_numerov(Grid,l,n,eig,rv,rhs,temp) 873 w(:,i)=temp(:) 874 ENDIF 875 ENDDO 876 877 OPEN(unit=2001,file='w.txt') 878 DO i=1,n 879 WRITE(2001,'(1p,50e15.7)') r(i), (w(i,j),j=1,nbase) 880 ENDDO 881 CLOSE(2001) 882 883 DEALLOCATE(rhs,temp) 884 END SUBROUTINE Calc_w 885 886 SUBROUTINE Calc_u(Grid,PAW,io,rhs,u) 887 TYPE(GridInfo), INTENT(IN) :: Grid 888 TYPE(PseudoInfo), INTENT(IN) :: PAW 889 INTEGER, INTENT(IN) :: io 890 REAL(8), INTENT(IN) :: rhs(:) 891 REAL(8), INTENT(OUT) :: u(:) 892 893 INTEGER :: i,j,l,n 894 REAL(8) :: eig,zeroval,term 895 REAL(8),POINTER :: rv(:),r(:) 896 REAL(8) ,ALLOCATABLE :: temp(:) 897 898 eig=PAW%eig(io) 899 l=PAW%l(io) 900 rv=>PAW%rveff 901 n=Grid%n 902 r=>Grid%r 903 904 ALLOCATE(temp(n)) 905 temp=0;u=0 906 907 CALL inhomo_bound_numerov(Grid,l,n,eig,rv,rhs,u) 908 909 OPEN(unit=2001,file='uio.txt') 910 DO i=1,n 911 WRITE(2001,'(1p,6e15.7)') r(i), u(i) 912 ENDDO 913 CLOSE(2001) 914 915 Deallocate(temp) 916 END SUBROUTINE Calc_u 917 918END Module exx_pseudo 919