1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2! This module contains the following active subroutines: 3! SetPAWOptions1, SetPAWOptions2, StoreTOCCWFN, Troullier, kerker, 4! nonncps, checkghosts, sethat, coretailselfenergy, setcoretail, 5! fixtcorewfn, selfhatpot, setbasis, makebasis_bloechl, 6! makebasis_custom, makebasis_modrrkj, readmatchradius, 7! makebasis_V_setvloc, formprojectors, bsolv, ftprod, fthatpot, 8! ftkin, ftvloc, unboundsep, boundsep, PStoAE, Set_PAW_MatrixElements, 9! logderiv, FindVlocfromVeff, SCFPAW, PAWIter_LDA, exploreparms, 10! EXPLORElogderiv, Report_PseudobasisRP, phase_unwrap 11! Check_overlap_of_projectors 12! smoothcore, settau 13! 14! 5/2018 phase_unwrap contributed by Casey Brock from Vanderbilt U. 15!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 16 17#if defined HAVE_CONFIG_H 18#include "config.h" 19#endif 20 21MODULE pseudo 22 23 USE io_tools 24 USE GlobalMath 25 USE atomdata 26 USE aeatom 27 USE excor 28 USE exx_pseudo 29 USE hf_pseudo 30 USE numerov_mod 31 USE paw_sub 32 USE pseudodata 33 USE pseudo_sub 34 USE radialDirac 35 USE radialsr 36 USE input_dataset_mod 37 38 IMPLICIT NONE 39 40 Type(Pseudoinfo), TARGET :: PAW 41 42 ! Parameters controlling PAW options 43 INTEGER,PRIVATE,PARAMETER :: BLOECHL=1, VANDERBILT=2, CUSTOM=3, MODRRKJ=7, HFPROJ=8 44 INTEGER,PRIVATE,PARAMETER :: BLOECHLPS=0, POLYNOM=1, POLYNOM2=2, RRKJ=3 45 INTEGER,PRIVATE,PARAMETER :: VANDERBILTORTHO=0, GRAMSCHMIDTORTHO=1 46 INTEGER,PRIVATE,PARAMETER :: SVDORTHO=2 47 INTEGER,PRIVATE,PARAMETER :: MTROULLIER=1, ULTRASOFT=2, BESSEL=3, KERKER_E=4, KERKER_P=5 48 INTEGER,PRIVATE,PARAMETER :: HARTREE_FOCK=4, SETVLOC=5 49 50 REAL(8),PRIVATE, PARAMETER :: coretailtol=1.d-12 51 INTEGER, PRIVATE :: coretailpoints=-1,besselopt=-1 52 INTEGER, PRIVATE :: Projectorindex=-1,PSindex=-1,Orthoindex=-1,Vlocalindex=-1 53 REAL(8), PRIVATE :: gaussparam=-1 54 55CONTAINS 56 57 SUBROUTINE SetPAWOptions1(ifen,Grid) 58 INTEGER, INTENT(IN) :: ifen 59 Type(Gridinfo), INTENT(IN) :: Grid 60 61 INTEGER :: n,i,j,l,irc 62 REAL(8):: h,rc 63 LOGICAL :: multi_rc 64 65 n=Grid%n 66 67! Maximum l for basis functions (from input dataset) 68 PAW%lmax=input_dataset%lmax 69 70! Cut-off radii (from input dataset) 71 multi_rc=(input_dataset%rc/=input_dataset%rc_shap.or.& 72& input_dataset%rc/=input_dataset%rc_vloc.or.& 73& input_dataset%rc/=input_dataset%rc_core.or.& 74& input_dataset%rc_shap/=input_dataset%rc_vloc.or.& 75& input_dataset%rc_shap/=input_dataset%rc_core.or.& 76& input_dataset%rc_vloc/=input_dataset%rc_core) 77 PAW%irc =FindGridIndex(Grid,input_dataset%rc) 78 PAW%irc_shap=FindGridIndex(Grid,input_dataset%rc_shap) 79 PAW%irc_vloc=FindGridIndex(Grid,input_dataset%rc_vloc) 80 PAW%irc_core=FindGridIndex(Grid,input_dataset%rc_core) 81 82 PAW%rc =Grid%r(PAW%irc) 83 PAW%rc_shap=Grid%r(PAW%irc_shap) 84 PAW%rc_vloc=Grid%r(PAW%irc_vloc) 85 PAW%rc_core=Grid%r(PAW%irc_core) 86 irc=PAW%irc;rc=PAW%rc 87 88 if (irc>n-Grid%ishift) stop 'error -- rc is too big !' 89 WRITE(STD_OUT,*) ' adjusted rc ',rc, Grid%r(irc) 90 WRITE(STD_OUT,*) ' irc,rc = ',irc,rc 91 if (multi_rc) then 92 WRITE(STD_OUT,*) ' adjusted rc_shape ',PAW%rc_shap 93 WRITE(STD_OUT,*) ' adjusted rc_vloc ',PAW%rc_vloc 94 WRITE(STD_OUT,*) ' adjusted rc_core ',PAW%rc_core 95 endif 96 WRITE(ifen,*) ' paw parameters: ' 97 WRITE(ifen,*) ' lmax = ',PAW%lmax 98 WRITE(ifen,*) ' rc = ',PAW%rc 99 WRITE(ifen,*) ' irc = ',PAW%irc 100 if (multi_rc) then 101 WRITE(ifen,*) ' rc_shape = ',PAW%rc_shap 102 WRITE(ifen,*) ' rc_vloc = ',PAW%rc_vloc 103 WRITE(ifen,*) ' rc_core = ',PAW%rc_core 104 endif 105 End Subroutine SetPAWOptions1 106 107 SUBROUTINE SetPAWOptions2(ifen,Grid,Orbit,Pot,success) 108 INTEGER, INTENT(IN) :: ifen 109 Type(Gridinfo), INTENT(IN) :: Grid 110 Type(OrbitInfo), INTENT(IN) :: Orbit 111 Type(PotentialInfo), INTENT(IN) :: Pot 112 LOGICAL , INTENT(OUT) :: success 113 114 INTEGER :: i,pdeg,l 115 REAL(8) :: qcut,x,y,e 116 117 success=.true. 118 119 !Pseudization and orthogonalization parameters (from input dataset) 120 pdeg=input_dataset%pseudo_polynom2_pdeg 121 qcut=input_dataset%pseudo_polynom2_qcut 122 IF (input_dataset%projector_type==PROJECTOR_TYPE_BLOECHL) Projectorindex=BLOECHL 123 IF (input_dataset%projector_type==PROJECTOR_TYPE_VANDERBILT) Projectorindex=VANDERBILT 124 IF (input_dataset%projector_type==PROJECTOR_TYPE_MODRRKJ) Projectorindex=MODRRKJ 125 IF (input_dataset%projector_type==PROJECTOR_TYPE_CUSTOM) Projectorindex=CUSTOM 126 IF (input_dataset%projector_type==PROJECTOR_TYPE_HF) Projectorindex=HFPROJ 127 IF (input_dataset%pseudo_type==PSEUDO_TYPE_BLOECHL) PSindex=BLOECHLPS 128 IF (input_dataset%pseudo_type==PSEUDO_TYPE_POLYNOM) PSindex=POLYNOM 129 IF (input_dataset%pseudo_type==PSEUDO_TYPE_POLYNOM2) PSindex=POLYNOM2 130 IF (input_dataset%pseudo_type==PSEUDO_TYPE_RRKJ) PSindex=RRKJ 131 IF (input_dataset%pseudo_type==PSEUDO_TYPE_HF) PSindex=HARTREE_FOCK 132 IF (input_dataset%ortho_type==ORTHO_TYPE_GRAMSCHMIDT) Orthoindex=GRAMSCHMIDTORTHO 133 IF (input_dataset%ortho_type==ORTHO_TYPE_VANDERBILT) Orthoindex=VANDERBILTORTHO 134 IF (input_dataset%ortho_type==ORTHO_TYPE_SVD) Orthoindex=SVDORTHO 135 IF (input_dataset%ortho_type==ORTHO_TYPE_HF) Orthoindex=-13 136 137 write(PAW%Proj_description,'("Projector type:")') 138 if (PSindex==BLOECHLPS) then 139 write(PAW%Proj_description,'(a," Bloechl")') trim(PAW%Proj_description) 140 else if (Projectorindex==MODRRKJ) then 141 write(PAW%Proj_description,'(a," modified RKKJ projectors")') & 142& trim(PAW%Proj_description) 143 else if (PSindex==POLYNOM) then 144 write(PAW%Proj_description,'(a," polynomial pseudization")') & 145& trim(PAW%Proj_description) 146 else if (PSindex==POLYNOM2) then 147 write(PAW%Proj_description,'(a," improved polynomial pseudization")') & 148& trim(PAW%Proj_description) 149 else if (PSindex==RRKJ) then 150 write(PAW%Proj_description,'(a," RRKJ pseudization")') & 151& trim(PAW%Proj_description) 152 else if (PSindex==HARTREE_FOCK) then 153 write(PAW%Proj_description,'(" HF projectors using Vanderbilt-like scheme")') 154 endif 155 156 if (Orthoindex==VANDERBILTORTHO) then 157 PAW%orthogonalization_scheme='vanderbilt' 158 write(PAW%Proj_description,'(a," + Vanderbilt ortho.")') & 159& trim(PAW%Proj_description) 160 else if (Orthoindex==GRAMSCHMIDTORTHO) then 161 PAW%orthogonalization_scheme='gramschmidt' 162 write(PAW%Proj_description,'(a," + Gram-Schmidt ortho.")') & 163& trim(PAW%Proj_description) 164 else if (Orthoindex==SVDORTHO) then 165 PAW%orthogonalization_scheme='svd' 166 write(PAW%Proj_description,'(a," + SVD ortho.")') & 167& trim(PAW%Proj_description) 168 end if 169 170 write(std_out,*) PAW%Proj_description 171 172 !Shape function parameters (from input dataset) 173 gaussianshapefunction=.false.;besselshapefunction=.false. 174 if (input_dataset%shapefunc_type==SHAPEFUNC_TYPE_GAUSSIAN) then 175 gaussianshapefunction=.true. 176 gaussparam=input_dataset%shapefunc_gaussian_param 177 CALL sethat(Grid,PAW,gaussparam=gaussparam) ! Gaussian shape function 178 write(PAW%Comp_description,& 179& '("Gaussian compensation charge shape with gausstol = ",1p,1e12.4)') gaussparam 180 else if (input_dataset%shapefunc_type==SHAPEFUNC_TYPE_BESSEL) then 181 besselshapefunction=.true. 182 CALL sethat(Grid,PAW,besselopt=i) ! Bessel shape function 183 if (PAW%irc_shap/=PAW%irc) then 184 write(PAW%Comp_description,& 185& '("Bessel compensation charge shape zeroed at ",1p,1e12.4)') PAW%rc_shap 186 else 187 write(PAW%Comp_description,& 188& '("Bessel compensation charge shape zeroed at rc")') 189 endif 190 else 191 CALL sethat(Grid,PAW) ! sinc^2 shape function 192 if (PAW%irc_shap/=PAW%irc) then 193 write(PAW%Comp_description,& 194& '("Sinc^2 compensation charge shape zeroed at ",1p,1e12.4)') PAW%rc_shap 195 else 196 write(PAW%Comp_description,& 197& '("Sinc^2 compensation charge shape zeroed at rc")') 198 endif 199 endif 200 201 !Core tolerance for HF (from input dataset) 202 PAW%coretol=MAX(input_dataset%hf_coretol,0.d0) 203 if (Projectorindex==HFPROJ) write(std_out,*) 'Resetting coretol to ', PAW%coretol 204 205 !Vlocal parameters (from input dataset) 206 pdeg=input_dataset%pseudo_polynom2_pdeg 207 qcut=input_dataset%pseudo_polynom2_qcut 208 IF (input_dataset%vloc_type==VLOC_TYPE_MTROULLIER) Vlocalindex=MTROULLIER 209 IF (input_dataset%vloc_type==VLOC_TYPE_ULTRASOFT) Vlocalindex=ULTRASOFT 210 IF (input_dataset%vloc_type==VLOC_TYPE_BESSEL) Vlocalindex=BESSEL 211 IF (input_dataset%vloc_type==VLOC_TYPE_SETVLOC) Vlocalindex=SETVLOC 212 IF (input_dataset%vloc_type==VLOC_TYPE_KERKER_EXPF) Vlocalindex=KERKER_E 213 IF (input_dataset%vloc_type==VLOC_TYPE_KERKER_POLY) Vlocalindex=KERKER_P 214 215 if (Vlocalindex==MTROULLIER) then 216 l=input_dataset%vloc_l ; e=input_dataset%vloc_ene 217 WRITE(PAW%Vloc_description,& 218& '("Vloc: Norm-conserving Troullier-Martins with l= ",i1,";e= ",1p,1e12.4)')l,e 219 endif 220 if (Vlocalindex==ULTRASOFT) then 221 l=input_dataset%vloc_l ; e=input_dataset%vloc_ene 222 WRITE(PAW%Vloc_description,& 223& '("Vloc: Non norm-conserving form with l= ",i1,";e= ",1p,1e12.4)')l,e 224 endif 225 if (Vlocalindex==BESSEL) then 226 WRITE(PAW%Vloc_description,& 227& '("Vloc: truncated form - Vps(r)=A.sin(qr)/r for r<rc")') 228 endif 229 if (Vlocalindex==SETVLOC) then 230 PAW%VlocCoef=input_dataset%vloc_setvloc_coef 231 i=FindGridIndex(Grid,input_dataset%vloc_setvloc_rad) 232 PAW%VlocRad=Grid%r(i) 233 WRITE(PAW%Vloc_description, & 234& '("Vloc == VlocCoef*shapfunc , VlocCoef,Rad = ",1p,2e15.7)') PAW%VlocCoef,PAW%VlocRad 235 PAW%vloc= 0.d0;y=PAW%VlocRad 236 PAW%vloc(1)=PAW%VlocCoef 237 do i=2,Grid%n 238 if (Grid%r(i)<PAW%VlocRad) then 239 PAW%vloc(i)=PAW%VlocCoef*(SIN(pi*Grid%r(i)/y)/(pi*Grid%r(i)/y))**2 240 endif 241 enddo 242 endif 243 if (Vlocalindex==KERKER_E.or.Vlocalindex==KERKER_P) then 244 l=input_dataset%vloc_l ; e=input_dataset%vloc_ene 245 if (Vlocalindex==KERKER_E) PAW%Vloc_description="Norm-conserving Exp Vloc" 246 if (Vlocalindex==KERKER_P) PAW%Vloc_description="Norm-conserving Poly Vloc" 247 WRITE(PAW%Vloc_description,'(a,"; l = ",i1,"; powers = ",4i3,"; e = ",1p,e12.3)')& 248& trim(PAW%Vloc_description),l,input_dataset%vloc_kerker_power(1:4),e 249 ENDIF 250 251 WRITE(STD_OUT,*) PAW%Vloc_description 252 253 IF (Vlocalindex==MTROULLIER.and.Projectorindex/=HFPROJ) & 254& CALL troullier(Grid,Pot,PAW,l,e) 255 IF (Vlocalindex==MTROULLIER.and.Projectorindex==HFPROJ) then 256 CALL make_hf_basis_only(Grid,Pot,PAW) 257 CALL troullier_HF(Grid,Pot,PAW,l,e) 258 ENDIF 259 IF (Vlocalindex==ULTRASOFT) CALL nonncps(Grid,Pot,PAW,l,e) 260 IF (Vlocalindex==BESSEL) CALL besselps(Grid,Pot,PAW) 261 IF (Vlocalindex==KERKER_E.or.Vlocalindex==KERKER_P) CALL kerker(Grid,Pot,PAW) 262 !! Note: if SETVLOC only HF or VANDERBILT schemes work 263 IF (Projectorindex==BLOECHL) THEN 264 CALL makebasis_bloechl(Grid,Pot,0) 265 ELSE IF (Projectorindex==CUSTOM.AND.PSindex==BLOECHLPS) THEN 266 CALL makebasis_bloechl(Grid,Pot,1) 267 ELSE IF (Projectorindex==VANDERBILT.OR.Projectorindex==CUSTOM) THEN 268 if (Vlocalindex==SETVLOC) then 269 Call makebasis_V_setvloc(Grid,Pot,PAW) 270 else 271 CALL makebasis_custom(Grid,Pot,PSindex,Orthoindex,pdeg,qcut) 272 endif 273 ELSE IF (Projectorindex==HFPROJ) THEN 274 CALL make_hf_tp_only(Grid,Pot,PAW) 275 ELSE IF (Projectorindex==MODRRKJ) THEN 276 CALL makebasis_modrrkj(Grid,Pot,Orthoindex,success) 277 ENDIF 278 279 WRITE(ifen,*) TRIM(PAW%Vloc_description) 280 WRITE(ifen,*) TRIM(PAW%Proj_description) 281 WRITE(ifen,*) TRIM(PAW%Comp_description) 282 283 CALL StoreTOCCWFN(PAW) 284 285 END SUBROUTINE SetPAWOptions2 286 287 SUBROUTINE StoreTOCCWFN(PAW) 288 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 289 290 INTEGER :: io,ib 291 292 do io=1,PAW%TOCCWFN%norbit 293 if (PAW%valencemap(io)>0) then 294 ib=PAW%valencemap(io) 295 PAW%TOCCWFN%wfn(:,io)=PAW%tphi(:,ib) 296 endif 297 enddo 298 END SUBROUTINE StoreTOCCWFN 299 !*************************************************************** 300 ! SUBROUTINE troullier(lmax,Grid,Pot) 301 ! Creates screened norm-conserving pseudopotential following 302 ! approach of N. Troullier and J. L. Martins, PRB 43, 1993 (1991) 303 ! Uses p(r)=a0+f(r); f(r)=SUMm(Coef(m)*r^(2*m), where 304 ! m=1,2..6 305 ! Psi(r) = r^(l+1)*exp(p(r)) 306 !*************************************************************** 307 SUBROUTINE Troullier(Grid,Pot,PAW,l,e) 308 TYPE(Gridinfo), INTENT(IN) :: Grid 309 TYPE(Potentialinfo), INTENT(IN) :: Pot 310 TYPE(Pseudoinfo), INTENT(INOUT) :: PAW 311 INTEGER,INTENT(IN) :: l 312 REAL(8),INTENT(IN) :: e 313 314 REAL(8), ALLOCATABLE :: VNC(:) 315 REAL(8) :: A0,A,B,B0,C,C0,D,F,S 316 REAL(8) :: Coef(6),Coef0,Coef0old 317 REAL(8) :: h,rc,delta,x,pp,dpp,ddpp,dddpp,ddddpp 318 REAL(8) :: gam,bet 319 INTEGER :: i,j,k,n,iter,nr,nodes,irc,ok,m,wavetype 320 INTEGER, PARAMETER :: niter=5000 321 REAL(8), PARAMETER :: small=1.0d-9 322 REAL(8), ALLOCATABLE :: wfn(:),p(:),dum(:) 323 REAL(8), POINTER :: r(:),rv(:) 324 CHARACTER(132) :: line 325 326 n=Grid%n 327 h=Grid%h 328 r=>Grid%r 329 rv=>Pot%rv 330 nr=min(PAW%irc_vloc+10,n) 331 irc=PAW%irc_vloc 332 rc=PAW%rc_vloc 333 334 ALLOCATE(VNC(n),wfn(nr),p(nr),dum(nr),stat=ok) 335 IF (ok /=0) THEN 336 WRITE(STD_OUT,*) 'Error in troullier -- in allocating wfn,p', nr,ok 337 STOP 338 ENDIF 339 340 !write(std_out,*) ' Troullier ', n,nr,irc 341 !call flush_unit(std_out) 342 if (scalarrelativistic) then 343 CALL unboundsr(Grid,Pot,nr,l,e,wfn,nodes) 344 else 345 CALL unboundsch(Grid,Pot%rv,Pot%v0,Pot%v0p,nr,l,e,wfn,nodes) 346 endif 347 348 IF (wfn(irc)<0) wfn=-wfn 349 dum(1:irc)=(wfn(1:irc)**2) 350 S=integrator(Grid,dum(1:irc),1,irc) 351 A0=LOG(wfn(irc)/(rc**(l+1))) 352 B0=(rc*Gfirstderiv(Grid,irc,wfn)/wfn(irc)-(l+1)) 353 C0=rc*(rv(irc)-rc*e)-B0*(B0+2*l+2) 354 D=-rc*(rv(irc)-rc*Gfirstderiv(Grid,irc,rv))-2*B0*C0-2*(l+1)*(C0-B0) 355 F=rc*(2*rv(irc)-rc*(2*Gfirstderiv(Grid,irc,rv) & 356& -rc*Gsecondderiv(Grid,irc,rv)))+& 357& 4*(l+1)*(C0-B0)-2*(l+1)*D-2*C0**2-2*B0*D 358 359 WRITE(STD_OUT,*) 'In troullier -- matching parameters',S,A0,B0,C0,D,F 360 361 delta=1.d10 362 iter=0 363 Coef0=0 364 365 DO WHILE(delta>small.AND.iter<=niter) 366 iter=iter+1 367 A=A0-Coef0 368 B=B0 369 C=C0 370 CALL EvaluateTp(l,A,B,C,D,F,coef) 371 372 dum=0 373 DO i=1,irc 374 x=(r(i)/rc)**2 375 p(i)=x*(Coef(1)+x*(Coef(2)+x*(Coef(3)+& 376& x*(Coef(4)+x*(Coef(5)+x*Coef(6)))))) 377 dum(i)=((r(i)**(l+1))*EXP(p(i)))**2 378 ENDDO 379 Coef0old=Coef0 380 381 x=integrator(Grid,dum(1:irc),1,irc) 382 Coef0=(LOG(S/x))/2 383 384 delta=ABS(Coef0-Coef0old) 385 !WRITE(STD_OUT,'(" VNC: iter Coef0 delta",i5,1p,2e15.7)') iter,Coef0,delta 386 ENDDO 387 388 WRITE(STD_OUT,*) ' VNC converged in ', iter,' iterations' 389 WRITE(STD_OUT,*) ' Coefficients -- ', Coef0,Coef(1:6) 390 ! 391 ! Now calculate VNC 392 OPEN(88,file='NC',form='formatted') 393 ! 394 VNC=0 395 DO i=2,nr 396 x=(r(i)/rc)**2 397 p(i)=Coef0+x*(Coef(1)+x*(Coef(2)+& 398& x*(Coef(3)+x*(Coef(4)+x*(Coef(5)+x*Coef(6)))))) 399 dpp=2*r(i)/(rc**2)*(Coef(1)+x*(2*Coef(2)+x*(3*Coef(3)+& 400& x*(4*Coef(4)+x*(5*Coef(5)+x*6*Coef(6)))))) 401 ddpp=(1/(rc**2))*(2*Coef(1)+x*(12*Coef(2)+x*(30*Coef(3)+& 402& x*(56*Coef(4)+x*(90*Coef(5)+x*132*Coef(6)))))) 403 dddpp=(r(i)/rc**4)*(24*Coef(2)+x*(120*Coef(3)+x*(336*Coef(4)+& 404& x*(720*Coef(5)+x*1320*Coef(6))))) 405 ddddpp=(1/(rc**4)*(24*Coef(2)+x*(360*Coef(3)+x*(1680*Coef(4)+& 406& x*(5040*Coef(5)+x*11880*Coef(6)))))) 407 IF (i==irc) THEN 408 WRITE(STD_OUT,*) 'check dp ', dpp, B0/rc 409 WRITE(STD_OUT,*) 'check ddp ', ddpp, C0/rc**2 410 WRITE(STD_OUT,*) 'check dddp', dddpp, D/rc**3 411 WRITE(STD_OUT,*) 'check ddddp', ddddpp, F/rc**4 412 ENDIF 413 VNC(i)=e+ddpp+dpp*(dpp+2*(l+1)/r(i)) 414 dum(i)=(r(i)**(l+1))*EXP(p(i)) 415 WRITE(88,'(1p,5e15.7)') r(i),wfn(i),dum(i),VNC(i)*r(i),rv(i) 416 ENDDO 417 CLOSE(88) 418 x=overlap(Grid,dum(1:irc),dum(1:irc),1,irc) 419 WRITE(STD_OUT,*) 'check norm ',x,S 420 421 VNC(irc:n)=rv(irc:n)/r(irc:n) 422 PAW%rveff(1:n)=VNC(1:n)*r(1:n) 423 424 DEALLOCATE(VNC,wfn,p,dum) 425 END SUBROUTINE troullier 426 427 !*************************************************************** 428 ! SUBROUTINE kerker(lmax,Grid,Pot) 429 ! Creates screened norm-conserving pseudopotential following 430 ! approach of G. P. Kerker, J. Phys. C. 13,L189-L194 (1980) 431 ! Uses p(r)=a0+f(r); f(r)=SUMi(Coef(i)*r^m(i)), where m(i) 432 ! are input powers 433 ! Psi(r) = r^(l+1)*exp(p(r)) if PStype = EXPF 434 ! Psi(r) = r^(l+1)*(p(r)) if PStype = POLY 435 !*************************************************************** 436 SUBROUTINE kerker(Grid,Pot,PAW) 437 TYPE(Gridinfo), INTENT(IN) :: Grid 438 TYPE(Potentialinfo), INTENT(IN) :: Pot 439 TYPE(Pseudoinfo), INTENT(INOUT) :: PAW 440 441 REAL(8), ALLOCATABLE :: VNC(:) 442 REAL(8) :: A0,A,B,C,D,S,Coef(4),Coef0,Coef0old 443 REAL(8) :: h,e,rc,delta,x,pp,dpp,ddpp,dddpp 444 REAL(8) :: gam,bet 445 INTEGER :: i,j,k,n,iter,nr,nodes,irc,l,ok,m(4),wavetype 446 INTEGER, PARAMETER :: EXPF=1, POLY=2 447 INTEGER, PARAMETER :: niter=5000 448 REAL(8), PARAMETER :: small=1.0d-12 449 CHARACTER(10) :: vtype 450 REAL(8), ALLOCATABLE :: wfn(:),p(:),dum(:) 451 REAL(8), POINTER :: r(:),rv(:) 452 453 !Read data from input dataset 454 IF (input_dataset%vloc_type==VLOC_TYPE_KERKER_EXPF) wavetype=EXPF 455 IF (input_dataset%vloc_type==VLOC_TYPE_KERKER_POLY) wavetype=POLY 456 l=input_dataset%vloc_l ; e=input_dataset%vloc_ene 457 m(1:4)=input_dataset%vloc_kerker_power(1:4) 458 459 n=Grid%n 460 n=Grid%n 461 h=Grid%h 462 r=>Grid%r 463 rv=>Pot%rv 464 nr=min(PAW%irc_vloc+10,n) 465 irc=PAW%irc_vloc 466 rc=PAW%rc_vloc 467 ALLOCATE(VNC(n),wfn(nr),p(nr),dum(nr),stat=ok) 468 IF (ok /=0) THEN 469 WRITE(STD_OUT,*) 'Error in kerker -- in allocating wfn,p', nr,ok 470 STOP 471 ENDIF 472 473 if (scalarrelativistic) then 474 CALL unboundsr(Grid,Pot,nr,l,e,wfn,nodes) 475 else 476 CALL unboundsch(Grid,Pot%rv,Pot%v0,Pot%v0p,nr,l,e,wfn,nodes) 477 endif 478 479 480 481 IF (wfn(irc)<0) wfn=-wfn 482 dum(1:irc)=(wfn(1:irc)**2) 483 S=integrator(Grid,dum(1:irc),1,irc) 484 IF (wavetype==EXPF) THEN 485 A0=LOG(wfn(irc)/(rc**(l+1))) 486 B=(rc*Gfirstderiv(Grid,irc,wfn)/wfn(irc)-(l+1)) 487 C=rc*(rv(irc)-rc*e)-B*(B+2*l+2) 488 D=-rc*(rv(irc)-rc*Gfirstderiv(Grid,irc,rv))-2*B*C-2*(l+1)*(C-B) 489 ENDIF 490 491 IF (wavetype==POLY) THEN 492 A0=(wfn(irc)/(rc**(l+1))) 493 B=(rc*Gfirstderiv(Grid,irc,wfn))/(rc**(l+1))-(l+1)*A0 494 C=rc*(rv(irc)-rc*e)*A0-2*(l+1)*B 495 D=-rc*(rv(irc)-rc*Gfirstderiv(Grid,irc,rv))*A0+2*(l+1)*(B-C)+& 496& rc*(rv(irc)-rc*e)*B 497 ENDIF 498 499 500 WRITE(STD_OUT,*) 'In kerker -- matching parameters',S,A0,B,C,D 501 502 delta=1.d10 503 iter=0 504 Coef0=0 505 506 DO WHILE(delta>small.AND.iter<=niter) 507 iter=iter+1 508 A=A0-Coef0 509 CALL EvaluateP(m,A,B,C,D,Coef) 510 511 dum=0 512 DO i=1,irc 513 x=(r(i)/rc) 514 p(i)=(x**m(1))*Coef(1)+(x**m(2))*Coef(2)+(x**m(3))*Coef(3)+(x**m(4))*Coef(4) 515 IF (wavetype==EXPF)dum(i)=((r(i)**(l+1))*EXP(p(i)))**2 516 IF (wavetype==POLY)dum(i)=(wfn(i))**2-((r(i)**(l+1))*(p(i)))**2 517 ENDDO 518 Coef0old=Coef0 519 IF (wavetype==EXPF) THEN 520 x=integrator(Grid,dum(1:irc),1,irc) 521 Coef0=(LOG(S/x))/2 522 ENDIF 523 IF (wavetype==POLY) THEN 524 gam=(2*l+3)*integrator(Grid,dum(1:irc),1,irc)/(rc**(2*l+3)) 525 bet=(2*l+3)*(Coef(1)/(2*l+3+m(1))+Coef(2)/(2*l+3+m(2))+& 526& Coef(3)/(2*l+3+m(3))+Coef(4)/(2*l+3+m(4))) 527 !WRITE(STD_OUT,'("VNC: iter -- bet,gam = ",i5,1p,4e15.7)') iter,bet,gam 528 x=bet**2+gam 529 Coef0old=Coef0 530 IF (x<0.d0) THEN 531 WRITE(STD_OUT,*) 'Warning in Kerker subroutine x = ',x 532 Coef0=Coef0+0.1*A0 533 ELSE 534 Coef0=SQRT(x)-bet 535 ENDIF 536 ENDIF 537 delta=ABS(Coef0-Coef0old) 538 !WRITE(STD_OUT,*) ' VNC: iter Coef0 delta', iter,Coef0,delta 539 ENDDO 540 541 WRITE(STD_OUT,*) ' VNC converged in ', iter,' iterations' 542 WRITE(STD_OUT,*) ' Coefficients -- ', Coef0,Coef(1:4) 543 ! 544 ! Now calculate VNC 545 OPEN(88,file='NC',form='formatted') 546 ! 547 VNC=0 548 DO i=1,nr 549 x=(r(i)/rc) 550 p(i)=Coef0+(x**m(1))*Coef(1)+(x**m(2))*Coef(2)+& 551& (x**m(3))*Coef(3)+(x**m(4))*Coef(4) 552 dpp=(m(1)*(x**(m(1)-1))*Coef(1)+m(2)*(x**(m(2)-1))*Coef(2)+& 553& m(3)*(x**(m(3)-1))*Coef(3)+m(4)*(x**(m(4)-1))*Coef(4))/rc 554 ddpp=(m(1)*(m(1)-1)*(x**(m(1)-2))*Coef(1)+& 555& m(2)*(m(2)-1)*(x**(m(2)-2))*Coef(2)+& 556& m(3)*(m(3)-1)*(x**(m(3)-2))*Coef(3)+& 557& m(4)*(m(4)-1)*(x**(m(4)-2))*Coef(4))/(rc**2) 558 dddpp=(m(1)*(m(1)-1)*(m(1)-2)*(x**(m(1)-3))*Coef(1)+& 559& m(2)*(m(2)-1)*(m(2)-2)*(x**(m(2)-3))*Coef(2)+& 560& m(3)*(m(3)-1)*(m(3)-2)*(x**(m(3)-3))*Coef(3)+& 561& m(4)*(m(4)-1)*(m(4)-2)*(x**(m(4)-3))*Coef(4))/(rc**3) 562 IF (i==irc) THEN 563 WRITE(STD_OUT,*) 'check dp ', dpp, B/rc 564 WRITE(STD_OUT,*) 'check ddp ', ddpp, C/rc**2 565 WRITE(STD_OUT,*) 'check dddp', dddpp, D/rc**3 566 ENDIF 567 IF (wavetype==EXPF) THEN 568 VNC(i)=e+ddpp+dpp*(dpp+2*(l+1)/r(i)) 569 dum(i)=(r(i)**(l+1))*EXP(p(i)) 570 ENDIF 571 IF (wavetype==POLY) THEN 572 VNC(i)=e+(ddpp+2*(l+1)*dpp/r(i))/p(i) 573 dum(i)=(r(i)**(l+1))*(p(i)) 574 ENDIF 575 WRITE(88,'(1p,5e15.7)') r(i),wfn(i),dum(i),VNC(i)*r(i),rv(i) 576 ENDDO 577 CLOSE(88) 578 x=overlap(Grid,dum(1:irc),dum(1:irc),1,irc) 579 WRITE(STD_OUT,*) 'check norm ',x,S 580 581 VNC(irc:n)=rv(irc:n)/r(irc:n) 582 PAW%rveff(1:n)=VNC(1:n)*r(1:n) 583 584 DEALLOCATE(VNC,wfn,p,dum) 585 END SUBROUTINE kerker 586 587 !*************************************************************** 588 ! SUBROUTINE nonncps(lmax,Grid,Pot) 589 ! Creates screened pseudopotential by inverting Schroedinger 590 ! equation from a pseudized radial wave function of the form: 591 ! Psi(r) = r**(l+1) * exp (a + b*r**2 + c*r**4 + d*r**6) 592 ! No norm-conserving condition is imposed on Psi 593 !*************************************************************** 594 SUBROUTINE nonncps(Grid,Pot,PAW,l,e) 595 TYPE(Gridinfo), INTENT(IN) :: Grid 596 TYPE(Potentialinfo), INTENT(IN) :: Pot 597 TYPE(Pseudoinfo), INTENT(INOUT) :: PAW 598 INTEGER,INTENT(IN) :: l 599 REAL(8),INTENT(IN) :: e 600 601 INTEGER :: i,irc,n,nr,ok,nodes,i1,i2,i3,i4 602 REAL(8) :: rc,x,y1,y2,y3,p0,p1,p2,p3,sgn 603 REAL(8) :: b(4),c(4),d(4),amat(4,4) 604 REAL(8),ALLOCATABLE :: VNC(:),wfn(:) 605 REAL(8),POINTER :: r(:),rv(:) 606 CHARACTER(132) :: line 607 608 !Polynomial definitions 609 p0(x,y1,y2,y3)=(x-y1)*(x-y2)*(x-y3) 610 p1(x,y1,y2,y3)=(x-y2)*(x-y3)+(x-y1)*(x-y3)+(x-y1)*(x-y2) 611 p2(x,y1,y2,y3)=2.0d0*((x-y1)+(x-y2)+(x-y3)) 612 p3(x,y1,y2,y3)=6.0d0 613 614 n=Grid%n 615 r=>Grid%r 616 rv=>Pot%rv 617 nr=min(PAW%irc_vloc+10,n) 618 irc=PAW%irc_vloc 619 rc=PAW%rc_vloc 620 621 ALLOCATE(VNC(n),wfn(nr),stat=ok) 622 IF (ok/=0) stop 'Error in uspseudo -- allocating arrays' 623 624 if (scalarrelativistic) then 625 CALL unboundsr(Grid,Pot,nr,l,e,wfn,nodes) 626 else 627 CALL unboundsch(Grid,Pot%rv,Pot%v0,Pot%v0p,nr,l,e,wfn,nodes) 628 endif 629 IF (wfn(irc)<0) wfn=-wfn 630 631 DO i=2,nr 632 wfn(i)=wfn(i)/r(i)**dble(l+1) 633 ENDDO 634 635 i1=irc-1;i2=i1+1;i3=i2+1;i4=i3+1 636 c(1)=wfn(i1)/p0(r(i1),r(i2),r(i3),r(i4)) 637 c(2)=wfn(i2)/p0(r(i2),r(i3),r(i4),r(i1)) 638 c(3)=wfn(i3)/p0(r(i3),r(i4),r(i1),r(i2)) 639 c(4)=wfn(i4)/p0(r(i4),r(i1),r(i2),r(i3)) 640 d(1)=c(1)*p0(rc,r(i2),r(i3),r(i4)) + c(2)*p0(rc,r(i3),r(i4),r(i1)) + & 641& c(3)*p0(rc,r(i4),r(i1),r(i2)) + c(4)*p0(rc,r(i1),r(i2),r(i3)) 642 d(2)=c(1)*p1(rc,r(i2),r(i3),r(i4)) + c(2)*p1(rc,r(i3),r(i4),r(i1)) + & 643& c(3)*p1(rc,r(i4),r(i1),r(i2)) + c(4)*p1(rc,r(i1),r(i2),r(i3)) 644 d(3)=c(1)*p2(rc,r(i2),r(i3),r(i4)) + c(2)*p2(rc,r(i3),r(i4),r(i1)) + & 645& c(3)*p2(rc,r(i4),r(i1),r(i2)) + c(4)*p2(rc,r(i1),r(i2),r(i3)) 646 d(4)=c(1)*p3(rc,r(i2),r(i3),r(i4)) + c(2)*p3(rc,r(i3),r(i4),r(i1)) + & 647& c(3)*p3(rc,r(i4),r(i1),r(i2)) + c(4)*p3(rc,r(i1),r(i2),r(i3)) 648 649 sgn=d(1)/abs(d(1));d(1:4)=d(1:4)*sgn 650 b(1)=log(d(1));b(2:4)=d(2:4) 651 amat(1,1)= 1.0d0 652 amat(2:4,1)= 0.0d0 653 amat(1,2)= rc**2 654 amat(2,2)= 2.0d0*d(1)*rc 655 amat(3,2)= 2.0d0*d(1) +2.0d0*d(2)*rc 656 amat(4,2)= 4.0d0*d(2) +2.0d0*d(3)*rc 657 amat(1,3)= rc**4 658 amat(2,3)= 4.0d0*d(1)*rc**3 659 amat(3,3)= 12.0d0*d(1)*rc**2+ 4.0d0*d(2)*rc**3 660 amat(4,3)= 24.0d0*d(1)*rc +24.0d0*d(2)*rc**2+4.0d0*d(3)*rc**3 661 amat(1,4)= rc**6 662 amat(2,4)= 6.0d0*d(1)*rc**5 663 amat(3,4)= 30.0d0*d(1)*rc**4+ 6.0d0*d(2)*rc**5 664 amat(4,4)= 120.0d0*d(1)*rc**3+60.0d0*d(2)*rc**4+6.0d0*d(3)*rc**5 665 666 CALL linsol(amat,b,4,4,4,4) 667 !write(std_out,*) 'Completed linsol with coefficients' 668 !write(std_out,'(1p,10e15.7)') (b(i),i=1,4) 669 670 PAW%rveff(1)=0.d0 671 DO i=2,irc-1 672 c(1)=2.0d0*b(2)*r(i)+ 4.0d0*b(3)*r(i)**3+ 6.0d0*b(4)*r(i)**5 673 c(2)=2.0d0*b(2) +12.0d0*b(3)*r(i)**2+30.0d0*b(4)*r(i)**4 674 PAW%rveff(i)=r(i)*(e+dble(2*l+2)*c(1)/r(i)+c(1)**2+c(2)) 675 ENDDO 676 PAW%rveff(irc:n)=rv(irc:n) 677 678 DEALLOCATE(VNC,wfn) 679 680 END SUBROUTINE nonncps 681 682 683 684 SUBROUTINE checkghosts(Grid,Orbit,FC,PAW) 685 TYPE(GridInfo), INTENT(IN) :: Grid 686 TYPE(OrbitInfo), INTENT(IN) :: Orbit 687 TYPE(FCInfo), INTENT(IN) :: FC 688 TYPE(PseudoInfo), INTENT(in) :: PAW 689 INTEGER :: l,nr,nodes,i,io 690 REAL(8) :: energy,h 691 REAL(8), POINTER :: r(:) 692 REAL(8), ALLOCATABLE :: wfn(:),VNC(:) 693 TYPE(PotentialInfo) :: Pot 694 695 h=Grid%h 696 r=>Grid%r 697 nr=min(PAW%irc+5,Grid%n) 698 call InitPot(Pot,nr) 699 ALLOCATE(VNC(nr),wfn(nr),stat=i) 700 IF (i/=0) STOP 'Error in checkghosts allocation' 701 POT%rv(1:nr)=PAW%rveff(1:nr) 702 call zeropot(Grid,POT%rv,POT%v0,POT%v0p) 703 704 DO l=0,PAW%lmax 705 DO io=1,Orbit%norbit 706 IF((.NOT.Orbit%iscore(io)).AND.(Orbit%l(io)==l)) THEN 707 energy=Orbit%eig(io) 708 WRITE(STD_OUT,*) 'Check for ghosts with l', l,energy 709 CALL unboundsch(Grid,Pot%rv,Pot%v0,Pot%v0p,nr,l,energy,wfn,nodes) 710 !DO i=1,nr 711 ! WRITE(l+17,'(1p,2e15.7)') Grid%r(i),wfn(i) 712 !ENDDO 713 WRITE(STD_OUT,*) ' Found # nodes = ', nodes 714 EXIT 715 ENDIF 716 ENDDO 717 ENDDO 718 719 CALL DestroyPot(Pot) 720 DEALLOCATE(VNC,wfn) 721 722 END SUBROUTINE checkghosts 723 724 SUBROUTINE sethat(Grid,PAW,gaussparam,besselopt) 725 TYPE(GridInfo), INTENT(IN) :: Grid 726 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 727 INTEGER,INTENT(IN), OPTIONAL :: besselopt 728 REAL(8),INTENT(IN), OPTIONAL :: gaussparam 729 730 INTEGER :: n,irc,irc_shap,i 731 REAL(8), POINTER :: r(:) 732 REAL(8) :: h,con,rc,rc_shap,selfen,d,dd,jbes1,jbes2,qr 733 REAL(8) :: al(2),ql(2) 734 735 n=Grid%n 736 h=Grid%h 737 irc=PAW%irc 738 rc=PAW%rc 739 irc_shap=PAW%irc_shap 740 rc_shap=PAW%rc_shap 741 r=>Grid%r 742 743 PAW%hatden=0 744 PAW%projshape=0 745 PAW%hatshape=0 746 PAW%projshape(1)=1 747 PAW%hatshape(1)=1 748 DO i=2,irc-1 749 PAW%projshape(i)=(SIN(pi*r(i)/rc)/(pi*r(i)/rc))**2 750 ENDDO 751 if(present(gaussparam)) then 752 d=rc_shap/SQRT(LOG(1.d0/gaussparam)) 753 PAW%gausslength=d 754 DO i=2,irc 755 PAW%hatshape(i)=EXP(-(r(i)/d)**2) 756 ENDDO 757 PAW%irc_shap=PAW%irc 758 PAW%rc_shap=PAW%rc 759 else if(present(besselopt)) then 760 call shapebes(al,ql,0,rc_shap) 761 DO i=1,irc_shap-1 762 qr=ql(1)*r(i);CALL jbessel(jbes1,d,dd,0,0,qr) 763 qr=ql(2)*r(i);CALL jbessel(jbes2,d,dd,0,0,qr) 764 PAW%hatshape(i)=al(1)*jbes1+al(2)*jbes2 765 ENDDO 766 else 767 DO i=2,irc_shap-1 768 PAW%hatshape(i)=(SIN(pi*r(i)/rc_shap)/(pi*r(i)/rc_shap))**2 769 ENDDO 770 endif 771 PAW%hatden(1:irc)=PAW%hatshape(1:irc)*(r(1:irc)**2) 772 773 ! normalize 774 if (.not.besselshapefunction) then 775 con=integrator(Grid,PAW%hatden,1,PAW%irc_shap) 776 WRITE(STD_OUT,*) ' check hatden normalization', con 777 PAW%hatden=PAW%hatden/con 778 endif 779 780 CALL poisson(Grid,con,PAW%hatden,PAW%hatpot,selfen) 781 WRITE(STD_OUT,*) 'Self energy for L=0 hat density ', selfen 782 783 END SUBROUTINE sethat 784 785 SUBROUTINE coretailselfenergy(Grid,PAW,ctctse,cthatse) 786 TYPE(GridInfo), INTENT(IN) :: Grid 787 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 788 REAL(8), INTENT(OUT) :: ctctse,cthatse 789 790 INTEGER :: i,irc,n 791 REAL(8) :: rc,h,x,y,z 792 REAL(8), allocatable :: d1(:),d2(:) 793 794 n=Grid%n 795 h=Grid%h 796 irc=PAW%irc 797 798 allocate(d1(n),d2(n),stat=i) 799 if (i /= 0) then 800 write(std_out,*) 'coretailselfenergy: allocation error -- ', n,i 801 stop 802 endif 803 804 x=integrator(Grid,PAW%tcore) 805 write(std_out,*) 'tcore charge ' , x 806 CALL poisson(Grid,x,PAW%tcore,d1,ctctse) 807 d2(2:n)=PAW%hatden(2:n)*d1(2:n)/Grid%r(2:n) 808 d2(1)=0 809 cthatse=integrator(Grid,d2(1:irc),1,PAW%irc_shap) 810 write(std_out,*) 'ctctse,cthatse = ', ctctse,cthatse 811 812 deallocate(d1,d2) 813 814 END SUBROUTINE coretailselfenergy 815 816 817! SUBROUTINE setcoretail(Grid,coreden) 818! TYPE(GridInfo), INTENT(IN) :: Grid 819! REAL(8), INTENT(IN) :: coreden(:) 820! 821! REAL(8) :: rc,h,x,y,z,u0,u2,u4 822! REAL(8), allocatable :: d1(:),d2(:) 823! INTEGER :: i,j,k,n,irc 824! 825! n=Grid%n 826! h=Grid%h 827! irc=PAW%irc_core 828! rc=PAW%rc_core 829! 830! allocate(d1(n),d2(n),stat=i) 831! if (i /= 0) then 832! write(std_out,*) 'setcoretail: allocation error -- ', n,i 833! stop 834! endif 835! CALL derivative(Grid,coreden,d1) 836! CALL derivative(Grid,d1,d2) 837! 838! x=coreden(irc) 839! y=d1(irc)*rc 840! z=d2(irc)*(rc*rc) 841! write(std_out,*) 'setcoretail: x,y,z = ', x,y,z 842! 843! u0=3*x - 9*y/8 + z/8 844! u2=-3*x + 7*y/4 - z/4 845! u4=x - 5*y/8 + z/8 846! 847! write(std_out,*) 'setcoretail: u0,u2,u4 = ', u0,u2,u4 848! 849! PAW%core=coreden 850! PAW%tcore=coreden 851! 852! do i=1,irc 853! x=(Grid%r(i)/rc)**2 854! PAW%tcore(i)= x*(u0+x*(u2+x*u4)) 855! enddo 856! 857! ! Find coretailpoints 858! z = integrator(Grid,coreden) 859! 860! coretailpoints=PAW%irc+Grid%ishift !! coretailpoints should be>=PAW%irc 861! do i=PAW%irc+Grid%ishift,n 862! if(ABS(z-integrator(Grid,coreden,1,i))<coretailtol) then 863! coretailpoints=i 864! exit 865! endif 866! enddo 867! write(std_out,*) 'coretailpoints = ',coretailpoints 868! PAW%coretailpoints=coretailpoints 869! 870! deallocate(d1,d2) 871! 872! END SUBROUTINE setcoretail 873 874 875!+++++++++++++++++++++++++++++++++++++++++++++++++++++++ 876! SUBROUTINE smoothcore 877! program to take array orig (all electron den or tau 878! and return smooth polynomial function for 0 \le r \le rc_core 879! with first and second derivatives continuous 880!+++++++++++++++++++++++++++++++++++++++++++++++++++++++ 881 SUBROUTINE smoothcore(Grid,orig,smooth) 882 TYPE(GridInfo), INTENT(IN) :: Grid 883 REAL(8), INTENT(IN) :: orig(:) 884 REAL(8), INTENT(INOUT) :: smooth(:) 885 886 REAL(8) :: rc,h,x,y,z,u0,u2,u4 887 REAL(8), allocatable :: d1(:),d2(:) 888 INTEGER :: i,j,k,n,irc 889 890 n=Grid%n 891 h=Grid%h 892 irc=PAW%irc_core 893 rc=PAW%rc_core 894 895 allocate(d1(n),d2(n),stat=i) 896 if (i /= 0) then 897 write(std_out,*) 'smoothcore: allocation error -- ', n,i 898 stop 899 endif 900 CALL derivative(Grid,orig,d1) 901 CALL derivative(Grid,d1,d2) 902 903 x=orig(irc) 904 y=d1(irc)*rc 905 z=d2(irc)*(rc*rc) 906 write(std_out,*) 'smoothcore: x,y,z = ', x,y,z 907 908 u0=3*x - 9*y/8 + z/8 909 u2=-3*x + 7*y/4 - z/4 910 u4=x - 5*y/8 + z/8 911 912 write(std_out,*) 'smoothcore: u0,u2,u4 = ', u0,u2,u4 913 914 smooth=orig 915 do i=1,irc 916 x=(Grid%r(i)/rc)**2 917 smooth(i)= x*(u0+x*(u2+x*u4)) 918 enddo 919 920 deallocate(d1,d2) 921 END SUBROUTINE smoothcore 922 923 924 SUBROUTINE setcoretail(Grid,coreden) 925 TYPE(GridInfo), INTENT(IN) :: Grid 926 REAL(8), INTENT(IN) :: coreden(:) 927 928 REAL(8) :: rc,h,x,y,z,u0,u2,u4 929 INTEGER :: i,j,k,n,irc 930 931 n=Grid%n 932 h=Grid%h 933 irc=PAW%irc_core 934 rc=PAW%rc_core 935 936 CALL smoothcore(Grid,coreden,PAW%tcore) 937 938 PAW%core=coreden 939 940 ! Find coretailpoints 941 z = integrator(Grid,coreden) 942 943 coretailpoints=PAW%irc+Grid%ishift !! coretailpoints should be>=PAW%irc 944 do i=PAW%irc+Grid%ishift,n 945 if(ABS(z-integrator(Grid,coreden,1,i))<coretailtol) then 946 coretailpoints=i 947 exit 948 endif 949 enddo 950 write(std_out,*) 'coretailpoints = ',coretailpoints 951 PAW%coretailpoints=coretailpoints 952 953 END SUBROUTINE setcoretail 954 955 SUBROUTINE setttau(Grid,coretau) 956 TYPE(GridInfo), INTENT(IN) :: Grid 957 REAL(8), INTENT(IN) :: coretau(:) 958 REAL(8) :: sqr4pi 959 960 write(std_out,*) 'in setttau ' 961 CALL smoothcore(Grid,coretau,PAW%tcoretau) 962 963 PAW%coretau=coretau 964 write(std_out,*) 'completed setttau' 965 966 sqr4pi=sqrt(4*pi)*1.d-10 967 PAW%itau=max(PAW%coretailpoints,1) 968 do while (PAW%itau<Grid%n.and. & 969& abs(PAW%tcoretau(PAW%itau))>sqr4pi*Grid%r(PAW%itau)**2) 970 PAW%itau=PAW%itau+1 971 end do 972 973 END SUBROUTINE setttau 974 975!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 976! Set smooth core functions for EXXKLI case using proceedure 977! identical to HF 978! 979!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 980 SUBROUTINE fixtcorewfn(Grid,PAW) 981 TYPE(GridInfo), INTENT(IN) :: Grid 982 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 983 984 INTEGER :: i,j,k,l,n,io,jo,ib,ic,irc 985 LOGICAL :: last 986 INTEGER , PARAMETER :: np=5 987 988 n=Grid%n 989 irc=PAW%irc 990 991 !Do i=1,n 992 ! write(444,'(1p,10e15.7)') Grid%r(i),PAW%tcore(i) 993 !enddo 994 Do io=1,PAW%TOCCWFN%norbit 995 write(std_out,*) 'fixtcore', io,PAW%valencemap(io),PAW%TOCCWFN%iscore(io) 996 IF(PAW%valencemap(io)<0) THEN ! core state 997 PAW%TOCCWFN%wfn(:,io)=0 998 If (PAW%TOCCWFN%iscore(io)) then 999 last=.true. ; l=PAW%TOCCWFN%l(io) 1000 ! check to find if this is the outermost core state for this l 1001 if (io<PAW%TOCCWFN%norbit) then 1002 do jo=io+1,PAW%TOCCWFN%norbit 1003 if (PAW%TOCCWFN%iscore(jo).and.PAW%TOCCWFN%l(jo)==l) & 1004& last=.false. 1005 enddo 1006 endif 1007 if (last) then 1008 IF (MAXVAL(ABS(PAW%OCCWFN%wfn(irc-np/2:irc+np/2,io))) & 1009& >PAW%coretol) THEN 1010 CALL Smoothfunc(Grid%r,l,np,irc,n,PAW%OCCWFN%wfn(:,io), & 1011& PAW%TOCCWFN%wfn(:,io)) 1012 write(std_out,*) 'last',io 1013 ENDIF 1014 endif 1015 else 1016 write(std_out,*) 'Something wrong -- should be core state ',& 1017& io,PAW%TOCCWFN%l(io),PAW%TOCCWFN%eig(io),PAW%TOCCWFN%occ(io) 1018 stop 1019 endif 1020 ENDIF 1021 Enddo 1022 1023 ! reset PAW%tcore to be consistent 1024 PAW%tcore=0 1025 do io=1,PAW%TOCCWFN%norbit 1026 if (PAW%TOCCWFN%iscore(io)) then 1027 PAW%tcore=PAW%tcore+PAW%TOCCWFN%occ(io)*(PAW%TOCCWFN%wfn(:,io)**2) 1028 endif 1029 enddo 1030 !Do i=1,n 1031 ! write(445,'(1p,10e15.7)') Grid%r(i),PAW%tcore(i) 1032 !enddo 1033 END SUBROUTINE fixtcorewfn 1034 1035 SUBROUTINE selfhatpot(Grid,PAW,l,eself) 1036 TYPE(GridInfo), INTENT(IN) :: Grid 1037 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 1038 INTEGER, INTENT(IN) :: l 1039 REAL(8), INTENT(OUT) :: eself 1040 1041 1042 INTEGER :: n,irc,i 1043 REAL(8), POINTER :: r(:) 1044 REAL(8), ALLOCATABLE :: den(:),a(:) 1045 REAL(8) :: h,con 1046 REAL(8) :: qr,jbes1,jbes2,dum1,dum2,al(2),ql(2) 1047 1048 n=Grid%n 1049 h=Grid%h 1050 r=>Grid%r 1051 irc=PAW%irc 1052 1053 ALLOCATE(den(n),a(n),stat=i) 1054 IF (i/=0) THEN 1055 WRITE(STD_OUT,*) 'Error in selfhatpot allocation',irc,i 1056 STOP 1057 ENDIF 1058 1059 if (besselshapefunction) then 1060 call shapebes(al,ql,l,PAW%rc_shap) 1061 DO i=1,PAW%irc_shap 1062 qr=ql(1)*r(i);CALL jbessel(jbes1,dum1,dum2,l,0,qr) 1063 qr=ql(2)*r(i);CALL jbessel(jbes2,dum1,dum2,l,0,qr) 1064 den(i)=(al(1)*jbes1+al(2)*jbes2)*r(i)**2 1065 ENDDO 1066 if (n>PAW%irc_shap) den(PAW%irc_shap+1:n)=0.d0 1067 else 1068 DO i=1,n 1069 den(i)=(r(i)**l)*PAW%hatden(i) 1070 ENDDO 1071 a(1:n)=den(1:n)*(r(1:n)**l) 1072 con=integrator(Grid,a,1,PAW%irc_shap) 1073 den=den/con 1074 endif 1075 1076 a=0 1077 1078 CALL apoisson(Grid,l,n,den,a) 1079 1080 ! apoisson returns a*r 1081 a(2:n)=a(2:n)/r(2:n) 1082 a(1)=0 1083 1084 eself=0.5d0*overlap(Grid,a,den) 1085 1086 DEALLOCATE(den,a) 1087 1088 END SUBROUTINE selfhatpot 1089 1090 1091!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1092! For diracrelativistic case, this subroutine is not yet ready 1093! Suggested future changes are added as comments, such as 1094! Upper component of wave function is loaded with renormalization 1095!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1096 SUBROUTINE setbasis(Grid,Pot,Orbit) 1097 TYPE(GridInfo), INTENT(IN) :: Grid 1098 TYPE(PotentialInfo), INTENT(IN) :: Pot 1099 TYPE(OrbitInfo), INTENT(INOUT) :: Orbit 1100 1101 1102 INTEGER :: n,irc,nbase,l,lmax,mxbase,lng,currentnode 1103 INTEGER :: i,j,k,io,ok,nbl,nr,nodes,ib,loop,niter,iter,ibasis_add 1104 REAL(8) :: h,rc,q00,energy,rat,delta,thisconv,qeff,tq,range 1105 REAL(8) :: ecoul,etxc,eexc 1106 REAL(8), POINTER :: r(:) 1107 1108 n=Grid%n 1109 h=Grid%h 1110 r=>Grid%r 1111 irc=PAW%irc 1112 nr=MIN(irc+100,n-100) 1113 rc=PAW%rc 1114 lmax=PAW%lmax 1115 1116 ! set AErefrv 1117 PAW%AErefrv=Pot%rv 1118 PAW%rvx=Pot%rvx 1119 1120 PAW%exctype=Orbit%exctype 1121 1122 nbase=PAW%nbase 1123 mxbase=nbase+5*max(1,PAW%lmax) 1124 1125 ! "filter" occupied states for long-range noise 1126 DO io=1,Orbit%norbit 1127 Call Filter(n,Orbit%wfn(:,io),machine_zero) 1128 ENDDO 1129 1130 call CopyOrbit(Orbit,PAW%OCCwfn) 1131 call CopyOrbit(Orbit,PAW%TOCCwfn) 1132 PAW%valencemap=-13 1133 IF (Orbit%exctype=='HF') THEN 1134 range=r(n) 1135 rat=-1.d30; k=1 1136 do io=1,Orbit%norbit 1137 if ((Orbit%occ(io)>1.d-5).and.Orbit%eig(io)>rat) then 1138 rat=Orbit%eig(io) 1139 k=io 1140 endif 1141 enddo 1142 do i=n,irc+1,-1 1143 if (ABS(Orbit%wfn(i,k))>1.d-4) then 1144 range=r(i) 1145 exit 1146 endif 1147 enddo 1148 write(std_out,*) 'range ', k,range; call flush_unit(std_out) 1149 ENDIF 1150 1151 WRITE(STD_OUT,*) ' basis functions:' 1152 WRITE(STD_OUT,*)' No. n l energy occ ' 1153 1154 nbase=0 ; ibasis_add=1 1155 DO l=0,lmax 1156 currentnode=-1 1157 nbl=0 1158 DO io=1,Orbit%norbit ! cycle through all configuration 1159 IF (Orbit%l(io).EQ.l) THEN 1160 currentnode=Orbit%np(io)-l-1 1161 IF (.NOT.Orbit%iscore(io)) THEN 1162 nbl=nbl+1 1163 nbase=nbase+1 1164 PAW%np(nbase)=Orbit%np(io) 1165 PAW%l(nbase)=l 1166 PAW%nodes(nbase)=PAW%np(nbase)-l-1 1167 write(std_out,*) 'l,nbase,node',l,nbase,currentnode 1168 PAW%eig(nbase)=Orbit%eig(io) 1169 PAW%occ(nbase)=Orbit%occ(io) 1170 PAW%phi(:,nbase)=Orbit%wfn(:,io) 1171 1172 if(diracrelativistic) then 1173 STOP 'Error -- setbasis subroutine not ready for diracrelativistic!' 1174 !rat=1.d0/sqrt(overlap(Grid,PAW%phi(:,nbase),PAW%phi(:,nbase))) 1175 !PAW%phi(:,nbase)=rat*PAW%phi(:,nbase) 1176 !PAW%kappa(nbase)=kappa 1177 endif 1178 1179 PAW%valencemap(io)=nbase 1180 IF(Orbit%exctype=='HF') then 1181 PAW%eig(nbase)=HF%lmbd(io,io) 1182 PAW%lmbd(:,nbase)=HF%lmbd(:,io) 1183 write(std_out,*) 'lmbd for nbase ', nbase 1184 write(std_out,'(1p,20e15.7)') PAW%lmbd(1:Orbit%norbit,nbase) 1185 PAW%lmbd(io,nbase)=0.d0 ! to avoid double counting 1186 endif 1187 WRITE(STD_OUT,'(3i6,1p,2e15.6)') nbase,PAW%np(nbase),l, & 1188& PAW%eig(nbase),PAW%occ(nbase); call flush_unit(std_out) 1189 ENDIF 1190 ENDIF 1191 ENDDO 1192 1193 generalizedloop: DO 1194 IF (ibasis_add>input_dataset%nbasis_add) EXIT generalizedloop 1195 IF (input_dataset%basis_add_l(ibasis_add)/=l) EXIT generalizedloop 1196 energy=input_dataset%basis_add_energy(ibasis_add) 1197 IF (energy<0.d0) & 1198& WRITE(STD_OUT,*) 'energy is negative',energy,' -- WARNING WARNING !!!' 1199 nbase=nbase+1 1200 IF (nbase > mxbase ) THEN 1201 WRITE(STD_OUT,*) 'Error in setbasis -- too many functions ', nbase,mxbase 1202 STOP 1203 ENDIF 1204 PAW%l(nbase)=l 1205 PAW%np(nbase)=999 1206 PAW%nodes(nbase)=currentnode+1 1207 currentnode=PAW%nodes(nbase) 1208 write(std_out,*) 'l,nbase,node',l,nbase,currentnode 1209 PAW%eig(nbase)=energy 1210 PAW%occ(nbase)=0.d0 1211 PAW%phi(1:n,nbase)=0.d0 1212 if (scalarrelativistic) then 1213 CALL unboundsr(Grid,Pot,n,l,energy,PAW%phi(:,nbase),nodes) 1214 elseif ( Orbit%exctype=='HF') then 1215 CALL HFunocc(Grid,Orbit,l,energy,Pot%rv,Pot%v0,Pot%v0p, & 1216& PAW%phi(:,nbase),PAW%rng(nbase)) 1217 else 1218 CALL unboundsch(Grid,Pot%rv,Pot%v0,Pot%v0p,& 1219& nr,l,energy,PAW%phi(:,nbase),nodes) 1220 endif 1221 rat=MAX(ABS(PAW%phi(irc,nbase)),ABS(PAW%phi(irc+1,nbase))) 1222 rat=DSIGN(rat,PAW%phi(irc,nbase)) 1223 PAW%phi(1:n,nbase)=PAW%phi(1:n,nbase)/rat 1224 !IF(Orbit%exctype=='HF') PAW%lmbd(:,nbase)=PAW%lmbd(:,nbase)/rat 1225 WRITE(STD_OUT,'(3i6,1p,2e15.6)') nbase,PAW%np(nbase),l, & 1226& PAW%eig(nbase),PAW%occ(nbase) 1227 nbl=nbl+1 1228 ibasis_add=ibasis_add+1 1229 ENDDO generalizedloop 1230 1231 1232 ENDDO ! end lmax loop 1233 1234 WRITE(STD_OUT,*) 'completed phi basis with ',nbase,' functions ' 1235 PAW%nbase=nbase ! reset nbase 1236 1237 END SUBROUTINE setbasis 1238 1239 !************************************************************************** 1240 ! Program to generate atomic basis functions 1241 ! Version using Bloechl's form of projector and orthogonalization procedure 1242 ! At the end of this subroutine, the basis functions and projectors are 1243 ! orthogonalized with a Gram-Schmidt like scheme 1244 !************************************************************************** 1245 SUBROUTINE makebasis_bloechl(Grid,Pot,option) 1246 TYPE(GridInfo), INTENT(IN) :: Grid 1247 TYPE(PotentialInfo), INTENT(IN) :: POT 1248 INTEGER,INTENT(IN) :: option 1249 1250 INTEGER :: n,irc,nr 1251 INTEGER :: i,j,k,io,ok,lmax 1252 REAL(8) :: h,rc 1253 REAL(8), ALLOCATABLE :: tmp(:),VNC(:) 1254 TYPE(PotentialInfo), TARGET:: PS 1255 REAL(8), POINTER :: r(:) 1256 1257 n=Grid%n 1258 h=Grid%h 1259 r=>Grid%r 1260 irc=PAW%irc 1261 nr=irc+20 1262 rc=PAW%rc 1263 lmax=PAW%lmax 1264 PS%nz=0.d0 1265 1266 ALLOCATE(tmp(n),VNC(n),PS%rv(n),stat=ok) 1267 IF (ok /= 0 ) THEN 1268 WRITE(STD_OUT,*) 'Error in allocating denout in makebasis ',n,ok 1269 STOP 1270 ENDIF 1271 1272 ! find basis functions 1273 PS%rv=PAW%rveff 1274 call zeropot(Grid,PS%rv,PS%v0,PS%v0p) 1275 !write(std_out,*) 'VNC v0 ', PS%v0,PS%v0p,PS%rv(5) 1276 CALL formprojectors(Grid,Pot,PS,option) 1277 1278 DEALLOCATE(tmp,VNC,PS%rv) 1279 do io=1,PAW%nbase 1280 PAW%rcio(io)=PAW%rc 1281 end do 1282 1283 END SUBROUTINE makebasis_bloechl 1284 1285 1286 !************************************************************************** 1287 ! Program to generate atomic basis functions 1288 ! 1289 ! 1) Pseudization of partial waves: 1290 ! - simple polynom scheme [optps=1] 1291 ! r^(l+1).Sum[Ci.r^2i] 0<=i<=4 1292 ! OR - ultrasoft polynom scheme [optps=2] 1293 ! r^(l+1).{Sum[Ci.r^2i]+Sum[Cj.r^2j]} 0<=i<=3 1294 ! 3<j adjusted using Fourier filtering 1295 ! OR - RRKJ scheme with 2 Bessel functions (PHYS REV B 41,1227 (1990)) [optps=3] 1296 ! 1297 ! 2) Build and orthogonalization of projectors 1298 ! - Vanderbilt generation method (PHYS REV B 41,7892 (1990)) [optorth=0] 1299 ! OR - Gram-Schmidt like sheme [optorth=1] 1300 !************************************************************************** 1301 SUBROUTINE makebasis_custom(Grid,Pot,optps,optorth,pdeg,qcut) 1302 TYPE(GridInfo), INTENT(IN) :: Grid 1303 TYPE(PotentialInfo), INTENT(IN) :: Pot 1304 INTEGER, INTENT(IN) :: optps,optorth,pdeg 1305 REAL(8), INTENT(IN) :: qcut 1306 1307 INTEGER :: i,j,k,l,io,jo,ok,lmax,nbase,n,irc,irc_vloc,nr,np,thisrc 1308 INTEGER :: icount,jcount,istart,ifinish,ibase,jbase 1309 REAL(8) :: choice,rc,xx,yy,gg,g,gp,gpp,al(2),ql(2) 1310 INTEGER, ALLOCATABLE :: omap(:) 1311 REAL(8), ALLOCATABLE :: VNC(:),Ci(:),aa(:,:),ai(:,:) 1312 REAL(8), POINTER :: r(:) 1313 1314 if (optps<1.or.optps>3) stop 'bug: error calling makebasis_custom routine' 1315 if (optorth<0.or.optorth>1) stop 'bug: error calling makebasis_custom routine' 1316 1317 n=Grid%n 1318 r=>Grid%r 1319 irc=PAW%irc 1320 irc_vloc=PAW%irc_vloc 1321 nbase=PAW%nbase 1322 lmax=PAW%lmax 1323 1324 np=5;if (optps==2) np=pdeg+1 ; 1325 if (optps==1.or.optps==2) allocate(Ci(np)) 1326 1327 ! Set screened local pseudopotential 1328 allocate(VNC(n),stat=i) 1329 if (i/=0) stop 'allocation error in makebasis_vanderbilt' 1330 VNC(2:n)=PAW%rveff(2:n)/r(2:n) 1331 call extrapolate(Grid,VNC) 1332 1333 ! Loop on basis elements 1334 do io=1,nbase 1335 l=PAW%l(io) 1336 1337 ! Read matching radius (from input dataset) 1338 thisrc=FindGridIndex(Grid,input_dataset%basis_func_rc(io)) 1339 thisrc=MIN(thisrc,irc) ! make sure rc<total rc 1340 rc=r(thisrc);PAW%rcio(io)=rc 1341 WRITE(STD_OUT,'(a,3i5,1p,1e15.7)') ' For this wfn: ',io,PAW%np(io),PAW%l(io),PAW%eig(io) 1342 WRITE(STD_OUT,'(a,f10.7)') ' >>> rc =', rc 1343 if (thisrc<3.or.thisrc>irc.or. & 1344& (optps==1.and.thisrc>n-3).or. & 1345& (optps==2.and.thisrc>n-6)) then 1346 write(std_out,*) 'rc out of range', thisrc,n,irc 1347 stop 1348 endif 1349 1350 ! Find partial wave pseudization 1351 if (optps==1) then 1352 call pspolyn(PAW%phi(:,io),Ci,r,l,np,thisrc,n) 1353 else if (optps==2) then 1354 call psuspolyn(PAW%phi(:,io),Ci,r,l,np,thisrc,n,qcut) 1355 else if (optps==3) then 1356 call psbes(PAW%phi(:,io),al,ql,Grid,l,thisrc,n) 1357 endif 1358 1359 ! Compute pseudized partial wave and unnormalized projector 1360 PAW%tphi(:,io)=PAW%phi(:,io) 1361 PAW%tp(:,io)=0.d0 1362 if (optps==1.or.optps==2) then 1363 do i=1,thisrc-1 1364 xx=r(i)*r(i) 1365 PAW%tphi(i,io)=Ci(1)+Ci(2)*xx 1366 gp=2.d0*Ci(2) 1367 gpp=2.d0*Ci(2) 1368 do j=3,np 1369 PAW%tphi(i,io)=PAW%tphi(i,io)+Ci(j)*xx**(j-1) 1370 gp=gp+dble(2*j-2)*Ci( j)*xx**(j-2) 1371 gpp=gpp+dble((2*j-2)*(2*j-3))*Ci(j)*xx**(j-2) 1372 enddo 1373 PAW%tphi(i,io)=PAW%tphi(i,io)*r(i)**(l+1) 1374 PAW%tp(i,io)=(dble(2*(l+1))*gp+gpp)*r(i)**(l+1)+(PAW%eig(io)-VNC(i))*PAW%tphi(i,io) 1375 enddo 1376 else if (optps==3) then 1377 PAW%tphi(1,io)=0.d0 1378 do i=2,thisrc-1 1379 xx=ql(1)*r(i) 1380 call jbessel(g,gp,gpp,l,2,xx) 1381 PAW%tphi(i,io)=al(1)*g*r(i) 1382 gg=al(1)*(2.d0*ql(1)*gp+ql(1)*xx*gpp) 1383 xx=ql(2)*r(i) 1384 call jbessel(g,gp,gpp,l,2,xx) 1385 PAW%tphi(i,io)=PAW%tphi(i,io)+al(2)*g*r(i) 1386 gg=gg+al(2)*(2.d0*ql(2)*gp+ql(2)*xx*gpp) 1387 PAW%tp(i,io)=(PAW%eig(io)-VNC(i)-dble(l*(l+1))/(r(i)**2))*PAW%tphi(i,io)+gg 1388 enddo 1389 endif 1390 1391 if (thisrc<irc_vloc) then 1392 do i=thisrc,irc_vloc-1 1393 gpp=Gsecondderiv(Grid,i,PAW%tphi(:,io)) 1394 PAW%tp(i,io)=(PAW%eig(io)-VNC(i)-dble(l*(l+1))/(r(i)**2))*PAW%tphi(i,io)+gpp 1395 enddo 1396 endif 1397 1398 ! If Gram-Schmidt orthogonalization, form normalized projector 1399 if (optorth==1) then 1400 xx=overlap(Grid,PAW%tp(:,io),PAW%tphi(:,io),1,max(thisrc,irc_vloc)) 1401 PAW%tp(:,io)=PAW%tp(:,io)/xx 1402 endif 1403 1404 enddo !nbase 1405 1406 deallocate(VNC);if (optps==1.or.optps==2) deallocate(Ci) 1407 1408 ! Form orthogonalized projector functions 1409 do io=1,nbase 1410 PAW%ophi(:,io)=PAW%phi(:,io) 1411 PAW%otphi(:,io)=PAW%tphi(:,io) 1412 PAW%Kop(1,io)=0 1413 PAW%Kop(2:n,io)=(PAW%eig(io)-Pot%rv(2:n)/Grid%r(2:n))*PAW%phi(2:n,io) 1414 if (optorth==1) PAW%otp(:,io)=PAW%tp(:,io) 1415 enddo 1416 1417 ! First option : VANDERBILT SCHEME 1418 if (optorth==0) then 1419 do l=0,lmax 1420 icount=0 1421 do io=1,nbase 1422 if (PAW%l(io)==l) icount=icount+1 1423 enddo 1424 write(std_out,*) 'For l = ', l,icount,' basis functions' 1425 if (icount==0) cycle 1426 allocate(aa(icount,icount),ai(icount,icount),omap(icount)) 1427 aa=0;icount=0 1428 do io=1,nbase 1429 if (PAW%l(io)==l) then 1430 icount=icount+1 1431 omap(icount)=io 1432 endif 1433 enddo 1434 do i=1,icount 1435 io=omap(i) 1436 do j=1,icount 1437 jo=omap(j) 1438 aa(i,j)=overlap(Grid,PAW%otphi(:,io),PAW%tp(:,jo),1,irc) 1439 enddo 1440 enddo 1441 ai=aa;call minverse(ai,icount,icount,icount) 1442 1443 do i=1,icount 1444 io=omap(i) 1445 PAW%ck(io)=ai(i,i) 1446 PAW%otp(:,io)=0 1447 do j=1,icount 1448 jo=omap(j) 1449 PAW%otp(:,io)=PAW%otp(:,io)+PAW%tp(:,jo)*ai(j,i) 1450 enddo 1451 enddo 1452 deallocate(aa,ai,omap) 1453 enddo 1454 1455 ! Second option : GRAM-SCHMIDT SCHEME 1456 else if (optorth==1) then 1457 DO l=0,lmax 1458 icount=0 1459 do io=1,nbase 1460 if (PAW%l(io)==l) icount=icount+1 1461 enddo 1462 if (icount==0) cycle 1463 allocate(aa(icount,icount),ai(icount,icount),omap(icount)) 1464 icount=0 1465 DO io=1,nbase 1466 IF (PAW%l(io)==l) THEN 1467 IF (icount==0) istart=io 1468 IF (icount>=0) ifinish=io 1469 icount=icount+1;omap(icount)=io 1470 ENDIF 1471 ENDDO 1472 DO ibase=istart,ifinish 1473 DO jbase=istart,ibase 1474 IF (jbase.LT.ibase) THEN 1475 xx=overlap(Grid,PAW%otp(:,jbase),PAW%otphi(:,ibase),1,irc) 1476 yy=overlap(Grid,PAW%otphi(:,jbase),PAW%otp(:,ibase),1,irc) 1477 PAW%ophi(1:n,ibase)=PAW%ophi(1:n,ibase)-PAW%ophi(1:n,jbase)*xx 1478 PAW%Kop(1:n,ibase)=PAW%Kop(1:n,ibase)-PAW%Kop(1:n,jbase)*xx 1479 PAW%otphi(1:n,ibase)=PAW%otphi(1:n,ibase)-PAW%otphi(1:n,jbase)*xx 1480 PAW%otp(1:n,ibase)=PAW%otp(1:n,ibase)-PAW%otp(1:n,jbase)*yy 1481 aa(ibase-istart+1,jbase-istart+1)=xx 1482 aa(jbase-istart+1,ibase-istart+1)=xx 1483 ELSE IF (jbase.EQ.ibase) THEN 1484 xx=overlap(Grid,PAW%otp(:,jbase),PAW%otphi(:,ibase),1,irc) 1485 choice=1.d0/SQRT(ABS(xx)) 1486 PAW%otp(1:n,ibase)=PAW%otp(1:n,ibase)*DSIGN(choice,xx) 1487 PAW%otphi(1:n,ibase)=PAW%otphi(1:n,ibase)*choice 1488 PAW%ophi(1:n,ibase)=PAW%ophi(1:n,ibase)*choice 1489 PAW%Kop(1:n,ibase)=PAW%Kop(1:n,ibase)*choice 1490 aa(ibase-istart+1,ibase-istart+1)=xx 1491 ENDIF 1492 ENDDO 1493 ENDDO 1494 ai=aa; !call minverse(aa,icount,icount,icount) 1495 do i=1,icount 1496 io=omap(i);PAW%ck(io)=ai(i,i) 1497 enddo 1498 deallocate(aa,ai,omap) 1499 ENDDO 1500 endif 1501 1502 END SUBROUTINE makebasis_custom 1503 1504!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1505!!! makebasis_modrrkj 1506!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1507 SUBROUTINE makebasis_modrrkj(Grid,Pot,Orthoindex,success) 1508 TYPE(GridInfo), INTENT(IN) :: Grid 1509 TYPE(PotentialInfo), INTENT(IN) :: Pot 1510 INTEGER, INTENT(IN) :: Orthoindex 1511 1512 INTEGER :: i,j,k,l,io,jo,ok,lmax,nbase,n,irc,irc_vloc,nr,np,thisrc 1513 INTEGER :: icount,jcount,istart,ifinish,ibase,jbase,lprev 1514 INTEGER :: match=5 1515 REAL(8) :: dk 1516 REAL(8) :: choice,rc,xx,yy,gg,g,gp,gpp,al(2),ql(2),root1,root2,logderiv 1517 INTEGER, ALLOCATABLE :: omap(:),rcindex(:) 1518 REAL(8), ALLOCATABLE :: VNC(:),Ci(:),aa(:,:),ai(:,:),jl(:,:),kv(:) 1519 REAL(8), ALLOCATABLE :: f(:),rcval(:) 1520 REAL(8), ALLOCATABLE :: U(:,:),VT(:,:),WORK(:),S(:),X(:,:) 1521 INTEGER :: LWORK 1522 REAL(8), POINTER :: r(:) 1523 LOGICAL :: success 1524 1525 success=.false. 1526 n=Grid%n 1527 r=>Grid%r 1528 irc=PAW%irc 1529 irc_vloc=PAW%irc_vloc 1530 nbase=PAW%nbase 1531 lmax=PAW%lmax 1532 1533 ALLOCATE(rcindex(nbase),rcval(nbase)) 1534 1535 ! Input matching radii for each basis function 1536 1537 call readmatchradius(Grid,rcindex,rcval) 1538 1539 1540 ! Set screened local pseudopotential 1541 allocate(VNC(n),stat=i) 1542 if (i/=0) stop 'allocation error in make_modrrkj' 1543 VNC(2:n)=PAW%rveff(2:n)/r(2:n) 1544 call extrapolate(Grid,VNC) 1545 1546 1547 allocate(kv(match),jl(n,match),f(n),Ci(match),aa(match,match)) 1548 if (i/=0) stop 'allocation error in make_nrecipe' 1549 ! Loop on basis elements 1550 lprev=-1;np=0 1551 do io=1,nbase 1552 1553 PAW%Kop(1,io)=0 1554 PAW%Kop(2:n,io)=(PAW%eig(io)-Pot%rv(2:n)/Grid%r(2:n))*PAW%phi(2:n,io) 1555 1556 l=PAW%l(io) 1557 if (l==lprev) then 1558 np=np+1 1559 else 1560 np=1 1561 lprev=l 1562 endif 1563 1564 thisrc=rcindex(io) 1565 rc=rcval(io);PAW%rcio(io)=rc 1566 1567 ! Set normalization for PAW%eig(io)>0 1568 if (PAW%eig(io)>0) then 1569 PAW%phi(:,io)=PAW%phi(:,io)/PAW%phi(thisrc,io) 1570 endif 1571 ! Find logderiv 1572 logderiv=Gfirstderiv(Grid,thisrc,PAW%phi(:,io))/PAW%phi(thisrc,io) 1573 write(std_out,*) 'logderiv ', io, l, logderiv 1574 1575 ! Find bessel function zeros 1576 call solvbes(kv,1.d0,0.d0,l,np) 1577 write(std_out,*) 'Searching for solution ',kv(np) 1578 if (np==1) then 1579 root1=0.001d0; root2=kv(1)-0.001d0; xx=0.5d0*(root1+root2) 1580 else 1581 root1=kv(np-1)+0.001d0; root2=kv(np)-0.001d0; xx=0.5d0*(root1+root2) 1582 endif 1583 icount=0 1584 Do 1585 call jbessel(g,gp,gpp,l,2,xx) 1586 yy=(logderiv-(1.d0/xx+gp/g))/(-1.d0/xx**2+gpp/g-(gp/g)**2) 1587 if (abs(yy)<1.d-6) then 1588 write(std_out,*) 'exiting Bessel loop', icount,xx,yy 1589 exit 1590 else 1591 if (xx+yy>root2) then 1592 xx=0.5*(xx+root2) 1593 else if (xx+yy<root1) then 1594 xx=0.5*(xx+root1) 1595 else 1596 xx=xx+yy 1597 endif 1598 endif 1599 !write(std_out,*) 'iter Bessel', xx,yy 1600 icount=icount+1 1601 if (icount>1000) then 1602 write(std_out,*) 'Giving up on Bessel root' 1603 return 1604 endif 1605 Enddo 1606 1607 success=.true. 1608 write(std_out,*) 'Found Bessel root at xx' , xx 1609 dk=(pi/40)/rc ! could be made adjustable 1610 write(std_out,*) 'dk value for this case ', dk 1611 kv=0 1612 kv(1)=xx/rc-dk*(match-1)/2 1613 write(std_out,*) '# kv 1', kv(1) 1614 Do i=2,match 1615 kv(i)=kv(i-1)+dk 1616 write(std_out,*) '# kv ', i, kv(i) 1617 enddo 1618 1619 Do i=1,match 1620 f(:)=kv(i)*r(:) 1621 call sphbes(l,n,f) 1622 jl(:,i)=f(:)*r(:) 1623 Enddo 1624 1625 ! Match bessel functions with wavefunctions 1626 1627 Ci=0; aa=0 1628 Do j=1,match 1629 i=match/2-j+1 1630 Ci(j)=PAW%phi(thisrc-i*1,io) ! 1 step should vary 1631 aa(j,1:match)=jl(thisrc-i*1,1:match) 1632 enddo 1633 1634 call SolveAXeqBM(match,aa,Ci,match-1) 1635 write(std_out,*) 'Completed SolveAXeqB with coefficients' 1636 write(std_out,'(1p,10e15.7)') (Ci(i),i=1,match) 1637 1638 ! Compute pseudized partial wave and unnormalized projector 1639 PAW%tphi(:,io)=PAW%phi(:,io); 1640 PAW%tp(:,io)=0.d0 1641 do i=1,thisrc-1 1642 PAW%tphi(i,io)=sum(Ci(1:match)*jl(i,1:match)) 1643 do j=1,match 1644 PAW%tp(i,io)=PAW%tp(i,io)+& 1645& Ci(j)*(-PAW%eig(io)+(kv(j)**2)+VNC(i))*jl(i,j) 1646 enddo 1647 enddo 1648 1649 if (thisrc<=irc_vloc) then 1650 do i=thisrc,irc_vloc-1 1651 gpp=Gsecondderiv(Grid,i,PAW%tphi(:,io)) 1652 PAW%tp(i,io)=(-PAW%eig(io)& 1653& +VNC(i)+dble(l*(l+1))/(r(i)**2))*PAW%tphi(i,io) -gpp 1654 enddo 1655 endif 1656 1657 1658 enddo !nbase 1659 1660 Deallocate(aa) 1661! 1662 Selectcase(Orthoindex) 1663 Case default ! Orthoindex=VANDERBILTORTHO 1664 !! Form orthogonalized projectors -- VANDERBILTORTHO 1665 write(std_out,*) ' Vanderbilt ortho' 1666 do l=0,lmax 1667 icount=0 1668 do io=1,nbase 1669 if (PAW%l(io)==l) icount=icount+1 1670 enddo 1671 if (icount==0) cycle 1672 write(std_out,*) 'For l = ', l,icount,' basis functions' 1673 allocate(aa(icount,icount),ai(icount,icount),omap(icount)) 1674 aa=0;icount=0 1675 do io=1,nbase 1676 if (PAW%l(io)==l) then 1677 icount=icount+1 1678 omap(icount)=io 1679 endif 1680 enddo 1681 do i=1,icount 1682 io=omap(i) 1683 PAW%otphi(:,io)=PAW%tphi(:,io) 1684 PAW%ophi(:,io)=PAW%phi(:,io) 1685 enddo 1686 do i=1,icount 1687 io=omap(i) 1688 do j=1,icount 1689 jo=omap(j) 1690 aa(i,j)=overlap(Grid,PAW%otphi(:,io),PAW%tp(:,jo),1,irc) 1691 enddo 1692 enddo 1693 ai=aa;call minverse(ai,icount,icount,icount) 1694 1695 do i=1,icount 1696 io=omap(i) 1697 PAW%ck(io)=ai(i,i) 1698 PAW%otp(:,io)=0 1699 do j=1,icount 1700 jo=omap(j) 1701 PAW%otp(:,io)=PAW%otp(:,io)+PAW%tp(:,jo)*ai(j,i) 1702 enddo 1703 enddo 1704 1705 write(std_out,*) 'Check otp for l = ', l 1706 do i = 1, icount 1707 io=omap(i) 1708 do j = 1, icount 1709 jo=omap(j) 1710 write(std_out,*) 'Overlap i j ', i,j, & 1711& overlap(Grid,PAW%otphi(:,io),PAW%otp(:,jo),1,irc) 1712 enddo 1713 enddo 1714 deallocate(aa,ai,omap) 1715 enddo 1716 1717 Case(GRAMSCHMIDTORTHO) ! Form orthonormal projectors -- GRAM-SCHMIDT SCHEME 1718 write(std_out,*) 'Gramschmidt ortho' 1719 DO l=0,lmax 1720 icount=0 1721 do io=1,nbase 1722 if (PAW%l(io)==l) icount=icount+1 1723 enddo 1724 if (icount==0) cycle 1725 allocate(aa(icount,icount),ai(icount,icount),omap(icount)) 1726 icount=0 1727 DO io=1,nbase 1728 IF (PAW%l(io)==l) THEN 1729 IF (icount==0) istart=io 1730 IF (icount>=0) ifinish=io 1731 icount=icount+1;omap(icount)=io 1732 ENDIF 1733 ENDDO 1734 DO ibase=istart,ifinish 1735 PAW%otp(:,ibase)=PAW%tp(:,ibase) 1736 PAW%otphi(:,ibase)=PAW%tphi(:,ibase) 1737 PAW%ophi(:,ibase)=PAW%phi(:,ibase) 1738 DO jbase=istart,ibase 1739 IF (jbase.LT.ibase) THEN 1740 xx=overlap(Grid,PAW%otp(:,jbase),PAW%otphi(:,ibase),1,irc) 1741 yy=overlap(Grid,PAW%otphi(:,jbase),PAW%otp(:,ibase),1,irc) 1742 PAW%ophi(1:n,ibase)=PAW%ophi(1:n,ibase)-PAW%ophi(1:n,jbase)*xx 1743 PAW%Kop(1:n,ibase)=PAW%Kop(1:n,ibase)-PAW%Kop(1:n,jbase)*xx 1744 PAW%otphi(1:n,ibase)=PAW%otphi(1:n,ibase)-PAW%otphi(1:n,jbase)*xx 1745 PAW%otp(1:n,ibase)=PAW%otp(1:n,ibase)-PAW%otp(1:n,jbase)*yy 1746 aa(ibase-istart+1,jbase-istart+1)=xx 1747 aa(jbase-istart+1,ibase-istart+1)=xx 1748 ELSE IF (jbase.EQ.ibase) THEN 1749 xx=overlap(Grid,PAW%otp(:,jbase),PAW%otphi(:,ibase),1,irc) 1750 !choice=1.d0/SQRT(ABS(xx)) 1751 !PAW%otp(1:n,ibase)=PAW%otp(1:n,ibase)*DSIGN(choice,xx) 1752 !PAW%otphi(1:n,ibase)=PAW%otphi(1:n,ibase)*choice 1753 !PAW%ophi(1:n,ibase)=PAW%ophi(1:n,ibase)*choice 1754 !PAW%Kop(1:n,ibase)=PAW%Kop(1:n,ibase)*choice 1755 PAW%otp(1:n,ibase)=PAW%otp(1:n,ibase)/xx 1756 aa(ibase-istart+1,ibase-istart+1)=xx 1757 ENDIF 1758 ENDDO 1759 ENDDO 1760 ai=aa; 1761 do i=1,icount 1762 io=omap(i);PAW%ck(io)=ai(i,i) 1763 enddo 1764 deallocate(aa,ai,omap) 1765 ENDDO 1766 1767 Case (SVDORTHO) ! SVD construction 1768 write(std_out,*) 'SVD ortho' 1769 DO l=0,lmax 1770 icount=0 1771 do io=1,nbase 1772 if (PAW%l(io)==l) icount=icount+1 1773 enddo 1774 if (icount==0) cycle 1775 allocate(aa(icount,icount),ai(icount,icount),omap(icount),X(n,icount)) 1776 icount=0 1777 DO io=1,nbase 1778 IF (PAW%l(io)==l) THEN 1779 IF (icount==0) istart=io 1780 IF (icount>=0) ifinish=io 1781 icount=icount+1;omap(icount)=io 1782 ENDIF 1783 ENDDO 1784 do i=1,icount 1785 io=omap(i) 1786 PAW%otphi(:,io)=0 1787 PAW%ophi(:,io)=0 1788 PAW%otp(:,io)=0 1789 do j=1,icount 1790 jo=omap(j) 1791 aa(i,j)=overlap(Grid,PAW%tphi(:,io),PAW%tp(:,jo),1,irc) 1792 enddo 1793 enddo 1794 LWORK=MAX(200,icount,icount) 1795 ALLOCATE(U(icount,icount),VT(icount,icount),WORK(LWORK),S(icount)) 1796 ai=aa; 1797 CALL DGESVD('A','A',icount,icount,ai,icount,S,& 1798& U,icount,VT,icount,WORK,LWORK,k) 1799 1800 write(std_out,*) 'For l = ', l , 'Completed SVD' 1801 write(std_out,*) 'S = ', S(:) 1802 1803 Do i=1,icount 1804 io=omap(i) 1805 X(:,i)=PAW%Kop(:,io); 1806 Enddo 1807 do j=1,icount 1808 jo=omap(j) 1809 PAW%ophi(:,jo)=0 1810 PAW%otphi(:,jo)=0 1811 PAW%otp(:,jo)=0 1812 PAW%Kop(:,jo)=0 1813 do i=1,icount 1814 io=omap(i) 1815 PAW%ophi(:,jo)=PAW%ophi(:,jo)+PAW%phi(:,io)*U(i,j) 1816 PAW%otphi(:,jo)=PAW%otphi(:,jo)+PAW%tphi(:,io)*U(i,j) 1817 PAW%Kop(:,jo)=PAW%Kop(:,jo)+X(:,i)*U(i,j) 1818 PAW%otp(:,jo)=PAW%otp(:,jo)+PAW%tp(:,io)*VT(j,i)/S(j) 1819 enddo 1820 enddo 1821 1822 do i=1,icount 1823 io=omap(i) 1824 do j=1,icount 1825 jo=omap(j) 1826 write(std_out,*) 'Overlap ', i,j, & 1827& overlap(Grid,PAW%otphi(:,io),PAW%otp(:,jo),1,irc) 1828 enddo 1829 enddo 1830 1831 deallocate(aa,ai,omap,U,VT,S,WORK,X) 1832 1833 ENDDO 1834 1835 End Select 1836 1837 Deallocate(VNC,kv,jl,f,Ci,rcindex,rcval) 1838 END SUBROUTINE makebasis_modrrkj 1839 1840 SUBROUTINE readmatchradius(Grid,rcindex,rcval) 1841 TYPE(GridInfo), INTENT(IN) :: Grid 1842 INTEGER, INTENT(INOUT) :: rcindex(:) 1843 REAL(8), INTENT(INOUT) :: rcval(:) 1844 1845 INTEGER :: io,nbase,n,irc,lmax,thisrc 1846 REAL(8) :: rc 1847 1848 nbase=PAW%nbase 1849 irc=PAW%irc 1850 1851 ! Loop on basis elements 1852 rcindex=0;rcval=0 1853 do io=1,nbase 1854 1855 ! Read matching radius (from input dataset) 1856 thisrc=FindGridIndex(Grid,input_dataset%basis_func_rc(io)) 1857 thisrc=MIN(thisrc,irc) ! make sure rc<total rc 1858 rc=Grid%r(thisrc) 1859 WRITE(STD_OUT,'(a,3i5,1p,1e15.7)') ' For this wfn: ',io,PAW%np(io),PAW%l(io),PAW%eig(io) 1860 WRITE(STD_OUT,'(a,f10.7)') ' >>> rc =', rc 1861 if (thisrc<3.or.thisrc>irc) then 1862 write(std_out,*) 'rc out of range', thisrc,n,irc 1863 stop 1864 endif 1865 rcindex(io)=thisrc;rcval(io)=rc 1866 write(std_out,'(" For io = ", i5," rc = ", i5,f10.5)') io,rcindex(io),rcval(io) 1867 1868 ENDDO 1869 END SUBROUTINE 1870 1871 1872!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1873 !modified version of SUBROUTINE makebasis_custom which was 1874 ! written by Marc Torrent from earlier version by NAWH 1875 ! for the case the vloc==0 (KS only; not HF) 1876 ! PAW%core and PAW%tcore already loaded 1877!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1878 SUBROUTINE makebasis_V_setvloc(Grid,Pot,PAW) 1879 Type(GridInfo), INTENT(IN) :: Grid 1880 Type(PotentialInfo), INTENT(IN) :: Pot 1881 Type(PseudoInfo), INTENT(INOUT) :: PAW 1882 1883 INTEGER :: i,j,k,l,n,ib,jb,io,jo,nbase,irc,np,thisrc,lmax,ni 1884 INTEGER :: icount 1885 INTEGER, ALLOCATABLE :: omap(:) 1886 REAL(8) :: rc,gp,gpp,xx,occ,stuff 1887 REAL(8), POINTER :: r(:) 1888 REAL(8), ALLOCATABLE :: Ci(:),dum(:),tdum(:),hat(:),aa(:,:),ai(:,:) 1889 REAL(8), ALLOCATABLE :: dpdr(:),pdr(:) 1890 REAL(8), PARAMETER :: threshold=1.d-6 1891 1892 if (TRIM(PAW%exctype)=='HF'.or.TRIM(PAW%exctype)=='EXX') then 1893 write(std_out,*) 'makebasis_V_setvloc is not designed for ', PAW%exctype 1894 stop 1895 endif 1896 n=Grid%n 1897 r=>Grid%r 1898 irc=PAW%irc 1899 lmax=PAW%lmax 1900 nbase=PAW%nbase 1901 np=5 1902 allocate(Ci(np),dum(n),tdum(n),hat(n)) 1903 PAW%rveff=PAW%vloc*Grid%r 1904 1905 ! Loop on basis elements 1906 do ib=1,nbase 1907 l=PAW%l(ib) 1908 ! Read matching radius (from input dataset) 1909 thisrc=FindGridIndex(Grid,input_dataset%basis_func_rc(ib)) 1910 thisrc=MIN(thisrc,irc) ! make sure rc<total rc 1911 rc=r(thisrc);PAW%rcio(io)=rc 1912 WRITE(STD_OUT,'(a,3i5,1p,1e15.7)') ' For this wfn: ',ib,PAW%np(ib),PAW%l(ib),PAW%eig(ib) 1913 WRITE(STD_OUT,'(a,f10.7)') ' >>> rc =', rc 1914 if (thisrc<3.or.thisrc>irc.or.thisrc>n-3) then 1915 write(std_out,*) 'rc out of range', thisrc,n,irc 1916 stop 1917 endif 1918 call pspolyn(PAW%phi(:,ib),Ci,r,l,np,thisrc,n) 1919 ! Compute pseudized partial wave and part of projector 1920 PAW%tphi(:,ib)=PAW%phi(:,ib) 1921 PAW%tp(:,ib)=0.d0 ; PAW%Kop(:,ib)=0 1922 do i=1,thisrc-1 1923 xx=r(i)*r(i) 1924 PAW%tphi(i,ib)=Ci(1)+Ci(2)*xx 1925 gp=2.d0*Ci(2) 1926 gpp=2.d0*Ci(2) 1927 do j=3,np 1928 PAW%tphi(i,ib)=PAW%tphi(i,ib)+Ci(j)*xx**(j-1) 1929 gp=gp+dble(2*j-2)*Ci( j)*xx**(j-2) 1930 gpp=gpp+dble((2*j-2)*(2*j-3))*Ci(j)*xx**(j-2) 1931 enddo 1932 PAW%tphi(i,ib)=PAW%tphi(i,ib)*r(i)**(l+1) 1933 PAW%tp(i,ib)=-(dble(2*(l+1))*gp+gpp)*r(i)**(l+1) !kinetic pt. 1934 PAW%Kop(i,ib)=PAW%tp(i,ib) 1935 enddo 1936 do i=thisrc,n 1937 gpp=Gsecondderiv(Grid,i,PAW%tphi(:,ib)) 1938 PAW%Kop(i,ib)=-gpp+dble(l*(l+1))/(r(i)**2)*PAW%tphi(i,ib) 1939 if (i<irc) PAW%tp(i,ib)=PAW%Kop(i,ib) 1940 enddo 1941 enddo !nbase 1942 1943 do io=1,PAW%OCCWFN%norbit 1944 if(PAW%valencemap(io)>0) then 1945 ib=PAW%valencemap(io) 1946 PAW%TOCCWFN%wfn(:,io)=PAW%tphi(:,ib) 1947 !PAW%OCCWFN%wfn(:,io)=PAW%phi(:,ib) !presumably also true 1948 endif 1949 enddo 1950 1951 ! Find rveff 1952 PAW%den=0.d0; PAW%tden=0.d0 1953 allocate(dpdr(n),pdr(n)) 1954 do io=1,PAW%OCCWFN%norbit 1955 if (.not.PAW%OCCwfn%iscore(io)) then 1956 occ=PAW%OCCwfn%occ(io) 1957 PAW%den=PAW%den+occ*PAW%OCCwfn%wfn(:,io)**2 1958 PAW%tden=PAW%tden+occ*PAW%TOCCwfn%wfn(:,io)**2 1959 endif 1960 enddo 1961 dum=PAW%core+PAW%den-PAW%tcore-PAW%tden 1962 xx=-Pot%nz+integrator(Grid,dum) 1963 write(std_out,*) 'Checking charge for pseudo scheme ', xx,integrator(Grid,tdum) 1964 PAW%rveff=PAW%rveff+xx*PAW%hatpot 1965 tdum=PAW%tcore+PAW%tden 1966 call poisson(Grid,gp,tdum,dum,stuff) ! Coulomb from tcore and tden 1967 write(std_out,*) 'After Poisson ', gp,stuff 1968 PAW%rveff=PAW%rveff+dum 1969 CALL exch(Grid,tdum,hat,gp,stuff) 1970 PAW%rveff=PAW%rveff+hat 1971 do i=1,n 1972 write(901,'(1p,20e15.7)') Grid%r(i),PAW%rveff(i),PAW%AErefrv(i),dum(i),& 1973& xx*PAW%hatpot(i),hat(i) 1974 enddo 1975 tdum=0; tdum(2:n)=PAW%rveff(2:n)/r(2:n) 1976 ! complete projector functions ; at this point PAW%tp stores -K*tilde{wfn} 1977 do ib=1,nbase 1978 PAW%tp(1:irc-1,ib)=PAW%tp(1:irc-1,ib)+& 1979& (tdum(1:irc-1)-PAW%eig(ib))*PAW%tphi(1:irc-1,ib) 1980 enddo 1981 1982 ! Form orthogonalized projector functions 1983 do ib=1,nbase 1984 PAW%ophi(:,ib)=PAW%phi(:,ib) 1985 PAW%otphi(:,ib)=PAW%tphi(:,ib) 1986 enddo 1987 1988 do l=0,lmax 1989 icount=0 1990 do ib=1,nbase 1991 if (PAW%l(ib)==l) icount=icount+1 1992 enddo 1993 if (icount==0) cycle 1994 write(std_out,*) 'For l = ', l,icount,' basis functions' 1995 allocate(aa(icount,icount),ai(icount,icount),omap(icount)) 1996 aa=0;icount=0 1997 do ib=1,nbase 1998 if (PAW%l(ib)==l) then 1999 icount=icount+1 2000 omap(icount)=ib 2001 endif 2002 enddo 2003 do i=1,icount 2004 ib=omap(i) 2005 do j=1,icount 2006 jb=omap(j) 2007 aa(i,j)=overlap(Grid,PAW%otphi(:,ib),PAW%tp(:,jb),1,irc) 2008 enddo 2009 enddo 2010 ai=aa;call minverse(ai,icount,icount,icount) 2011 2012 do i=1,icount 2013 ib=omap(i) 2014 PAW%ck(ib)=ai(i,i) 2015 PAW%otp(:,ib)=0 2016 do j=1,icount 2017 jb=omap(j) 2018 PAW%otp(:,ib)=PAW%otp(:,ib)+PAW%tp(:,jb)*ai(j,i) 2019 enddo 2020 enddo 2021 deallocate(aa,ai,omap) 2022 enddo 2023 2024 do i=1,n 2025 write(801,'(1p,50e15.7)') r(i),(PAW%otp(i,ib),PAW%tp(i,ib),ib=1,nbase) 2026 enddo 2027 2028 deallocate(Ci,dum,tdum,hat) 2029 END SUBROUTINE makebasis_V_setvloc 2030 2031 !************************************************************************* 2032 ! program to generate projector functions for Blochl's paw formalism 2033 ! starting with smooth functions 2034 ! for every basis function phi, choose smooth function with 2035 ! form for r<rc: (r**(l+1))*sum(n)*(cn*(r**(2*n))) 2036 ! where, rc==r(iiirc) 2037 ! phi(i,io), eval(io), tocc(io) original basis function, energy, occupancy 2038 ! tphi(i,io) smooth basis function corresponding to phi(i,io) 2039 ! tp(i,io)=constant*normalized Harmonic Oscillator function 2040 ! tp == 0 for r>r(iiirc) 2041 ! functions defined to be identically zero for r>r(iiirc)) 2042 !************************************************************************* 2043 SUBROUTINE formprojectors(Grid,Pot,PS,option) 2044 TYPE(GridInfo), INTENT(IN) :: Grid 2045 TYPE(PotentialInfo), INTENT(IN) :: Pot,PS 2046 INTEGER,INTENT(IN) :: option 2047 2048 INTEGER :: nbase,lmax,l,io,irc,wantednodes,nb,n,icount 2049 INTEGER :: istart,ifinish,ibase,jbase 2050 REAL(8) :: h,xx,yy,choice,rc 2051 INTEGER,allocatable :: irc_io(:) 2052 2053 n=Grid%n 2054 h=Grid%h 2055 lmax=PAW%lmax 2056 nbase=PAW%nbase 2057 irc=PAW%irc 2058 ! 2059 ! form projector functions for each l 2060 ! 2061 2062 !Read rc from input dataset 2063 if (option==1) then 2064 allocate(irc_io(nbase)) 2065 do io=1,nbase 2066 irc_io(io)=FindGridIndex(Grid,input_dataset%basis_func_rc(io)) 2067 rc=Grid%r(irc_io(io)) 2068 if(irc_io(io)>PAW%irc) then 2069 write(std_out,*) 'rc out of range', irc_io(io),n,PAW%irc 2070 stop 2071 endif 2072 enddo 2073 endif 2074 2075 DO l=0,lmax 2076 icount=0 2077 DO io=1,nbase 2078 IF (PAW%l(io)==l) THEN 2079 IF (icount==0) istart=io 2080 IF (icount >= 0) ifinish=io 2081 icount=icount+1 2082 wantednodes=icount-1 2083 ! form unorthonormalized projector functions tp 2084 WRITE(STD_OUT,*) '******* projector for l = ',l 2085 if (option==1) then 2086 CALL bsolv(Grid,PS,io,wantednodes,irc_io(io)) 2087 else 2088 CALL bsolv(Grid,PS,io,wantednodes) 2089 endif 2090 PAW%ophi(:,io)=PAW%phi(:,io) 2091 PAW%otphi(:,io)=PAW%tphi(:,io) 2092 PAW%otp(:,io)=PAW%tp(:,io) 2093 PAW%Kop(1,io)=0 2094 PAW%Kop(2:n,io)=(PAW%eig(io)-Pot%rv(2:n)/Grid%r(2:n))& 2095& *PAW%phi(2:n,io) 2096 ENDIF 2097 ENDDO 2098 !write(std_out,*) 'orthnormalization' 2099 !write(std_out,*) 'start orthogonalization',istart,ifinish 2100 DO ibase=istart,ifinish 2101 DO jbase=istart,ibase 2102 IF (jbase.LT.ibase) THEN 2103 xx=overlap(Grid,PAW%otp(:,jbase),PAW%otphi(:,ibase),1,irc) 2104 yy=overlap(Grid,PAW%otphi(:,jbase),PAW%otp(:,ibase),1,irc) 2105 !write(std_out,*) 'before',jbase,ibase,xx,yy 2106 PAW%ophi(1:n,ibase)=PAW%ophi(1:n,ibase)-PAW%ophi(1:n,jbase)*xx 2107 PAW%Kop(1:n,ibase)=PAW%Kop(1:n,ibase)-PAW%Kop(1:n,jbase)*xx 2108 PAW%otphi(1:n,ibase)=PAW%otphi(1:n,ibase)-PAW%otphi(1:n,jbase)*xx 2109 PAW%otp(1:n,ibase)=PAW%otp(1:n,ibase)-PAW%otp(1:n,jbase)*yy 2110 ELSE IF (jbase.EQ.ibase) THEN 2111 xx=overlap(Grid,PAW%otp(:,ibase),PAW%otphi(:,ibase),1,irc) 2112 !write(std_out,*) 'before',jbase,ibase,xx 2113 choice=1.d0/SQRT(ABS(xx)) 2114 PAW%otp(1:n,ibase)=PAW%otp(1:n,ibase)*DSIGN(choice,xx) 2115 PAW%otphi(1:n,ibase)=PAW%otphi(1:n,ibase)*choice 2116 PAW%ophi(1:n,ibase)=PAW%ophi(1:n,ibase)*choice 2117 PAW%Kop(1:n,ibase)=PAW%Kop(1:n,ibase)*choice 2118 ENDIF 2119 ENDDO 2120 ENDDO 2121 2122 ENDDO 2123 2124 if (option==1) deallocate(irc_io) 2125 2126 END SUBROUTINE formprojectors 2127 !************************************************************************* 2128 ! on input tphi=phi 2129 ! on output tphi recalculated for r<nrc*h 2130 ! l = angular momentum quantum number 2131 ! nodes = the number of desired nodes in tp and tphi 2132 ! shape function 2133 !************************************************************************* 2134 SUBROUTINE bsolv(Grid,Pot,io,wantednodes,irc_io) 2135 TYPE(GridInfo), INTENT(IN) :: Grid 2136 TYPE(PotentialInfo), INTENT(IN) :: Pot 2137 INTEGER, INTENT(IN) :: io,wantednodes 2138 INTEGER, OPTIONAL :: irc_io 2139 2140 REAL(8), ALLOCATABLE :: f(:),chi(:),fakerv(:) 2141 REAL(8), ALLOCATABLE :: tphi(:),tp(:) 2142 REAL(8), POINTER :: r(:),rv(:) 2143 REAL(8) :: h,en,v0,v0p,ol,h2,eh2,veh2,xnorm,rc_io 2144 INTEGER :: l,n,num,nrc,i,j,k,ii,jj,iter,irc,nodes,ok 2145 REAL(8) :: cpmin,cpmax,zeroval 2146 REAL(8) :: cp,cph2,del,deriv,tderiv 2147 REAL(8), PARAMETER :: small=1.d-10,step=1.d0 2148 INTEGER, PARAMETER :: mxiter=1000 2149 2150 n=Grid%n 2151 h=Grid%h 2152 r=>Grid%r 2153 rv=>Pot%rv 2154 if (present(irc_io)) then 2155 irc=irc_io 2156 else 2157 irc=PAW%irc 2158 endif 2159 ALLOCATE(fakerv(n),f(n),chi(n),tphi(n),tp(n),stat=ok) 2160 IF (ok /= 0) THEN 2161 WRITE(STD_OUT,*) 'Error in bsolv allocation ', n,ok 2162 STOP 2163 ENDIF 2164 tphi=PAW%phi(:,io) 2165 en=PAW%eig(io) 2166 l=PAW%l(io) 2167 deriv=(Gfirstderiv(Grid,irc,tphi))/tphi(irc) ! log derivative 2168 tderiv=0 ! initial log derive 2169 cp=0 ! initial projector constant 2170 del=1000 ! initial error 2171 cpmin=-100 2172 cpmax=100 2173 if (present(irc_io)) then 2174 rc_io=r(irc_io) 2175 f(1)=1.d0;f(irc_io:n)=0.d0 2176 DO i=2,irc_io-1 2177 f(i)=(SIN(pi*r(i)/rc_io)/(pi*r(i)/rc_io))**2 2178 ENDDO 2179 else 2180 DO i=1,n 2181 f(i)=PAW%projshape(i) 2182 ENDDO 2183 endif 2184 2185 WRITE(STD_OUT,*) 'in bsolv -- l, en, n',l,en,wantednodes 2186 2187 iter=0 2188 2189 DO 2190 2191 iter=iter+1 2192 2193 WRITE(STD_OUT,*) 'bsolv iter cp',iter,cp 2194 2195 fakerv=rv-cp*f*Grid%r 2196 call zeropot(Grid,fakerv,v0,v0p) 2197 ! initialize chi 2198 chi=0 2199 chi(2)=wfninit(0.d0,l,v0,v0p,en,Grid%r(2)) 2200 zeroval=0 2201 if (l==1) zeroval=2 2202 2203 call forward_numerov(Grid,l,irc+5,en,fakerv,zeroval,chi,nodes) 2204 2205 WRITE(STD_OUT,'("iter nodes cp cpmin cpmax",2i5,1p,3e15.7)')& 2206& iter,nodes,cp,cpmin,cpmax 2207 IF (nodes.EQ.wantednodes) THEN 2208 tderiv=(Gfirstderiv(Grid,irc,chi))/chi(irc) ! log derivative 2209 tp=0 2210 tp(1:irc)=chi(1:irc)/chi(irc) 2211 chi(1:irc)=f(1:irc)*(tp(1:irc))**2 2212 xnorm=integrator(Grid,chi(1:irc),1,irc) 2213 del=(tderiv-deriv)/xnorm 2214 WRITE(STD_OUT,*) 'iter nodes del', iter,nodes,del 2215 2216 IF (ABS(del).LT.small) EXIT 2217 2218 IF (iter.GE.mxiter) THEN 2219 WRITE(STD_OUT,*)' terminating projector',iter 2220 STOP 2221 ENDIF 2222 2223 IF (ABS(del).GT.step) del=DSIGN(step,del) 2224 cp=cp+del 2225 IF (cp.GT.cpmax) cp=cpmax-ranx()*step 2226 IF (cp.LT.cpmin) cp=cpmin+ranx()*step 2227 2228 ELSE IF(nodes.GT.wantednodes) THEN 2229 cpmax=cp 2230 cp=cpmax-ranx()*step 2231 ELSE IF(nodes.LT.wantednodes) THEN 2232 cpmin=cp 2233 cp=cpmin+ranx()*step 2234 ENDIF 2235 2236 ENDDO 2237 2238 tphi(1:irc-1)=tphi(irc)*tp(1:irc-1) 2239 tphi(irc:n)=PAW%phi(irc:n,io) 2240 tp(1:irc)=f(1:irc)*tphi(1:irc) 2241 chi(1:irc)=tphi(1:irc)*tp(1:irc) 2242 xnorm=integrator(Grid,chi(1:irc),1,irc) 2243 WRITE(STD_OUT,*) 'normalization for projector l,n=',l,nodes,xnorm 2244 tp(1:irc)=tp(1:irc)/xnorm 2245 tp(irc+1:n)=0 2246 2247 PAW%tphi(:,io)=tphi 2248 PAW%tp(:,io)=tp 2249 PAW%ck(io)=cp 2250 2251 WRITE(STD_OUT,*) 'completed bsolv',io,cp 2252 2253 DEALLOCATE(fakerv,f,chi,tphi,tp) 2254 END SUBROUTINE bsolv 2255 2256 !*********************************************************************** 2257 ! program to calculate Fourier transform of paw product tp*tphi 2258 ! in order to get an idea of the convergence 2259 ! and output them to files tprod.l 2260 !*********************************************************************** 2261 SUBROUTINE ftprod(Grid) 2262 TYPE(GridInfo), INTENT(IN) :: Grid 2263 2264 INTEGER, PARAMETER :: nq=200 2265 REAL(8), PARAMETER :: qmax=15.d0 2266 2267 REAL(8), ALLOCATABLE :: q(:),dum(:),dum1(:) 2268 REAL(8), POINTER :: r(:) 2269 REAL(8) :: h,dq,tphij 2270 INTEGER :: i,ib,l,iq,n,irc 2271 CHARACTER(4) flnm 2272 ! 2273 WRITE(STD_OUT,*) 'calculating Fourier transforms of tp*tphi products ',& 2274& 'For bound states only ' 2275 2276 n=Grid%n 2277 h=Grid%h 2278 r=>Grid%r 2279 irc=PAW%irc 2280 ALLOCATE(q(nq),dum(n),dum1(n),stat=i) 2281 IF (i/=0) THEN 2282 WRITE(STD_OUT,*) 'Error in allocating space in ftprod',n,nq,i 2283 STOP 2284 ENDIF 2285 2286 dq=qmax/nq 2287 DO ib=1,PAW%nbase 2288 IF (PAW%eig(ib)< 0.d0) THEN 2289 l=PAW%l(ib) 2290 CALL mkname(ib,flnm) 2291 OPEN(55,file='tprod.'//TRIM(flnm),form='formatted') 2292 DO iq=1,nq 2293 q(iq)=iq*dq 2294 DO i=1,n 2295 dum(i)=q(iq)*r(i) 2296 ENDDO 2297 CALL sphbes(l,n,dum) 2298 DO i=1,n 2299 dum1(i)=dum(i)*r(i)*PAW%otphi(i,ib) 2300 IF (i.LE.irc) dum(i)=dum(i)*r(i)*PAW%otp(i,ib) 2301 ENDDO 2302 tphij=integrator(Grid,dum1(1:n))*& 2303& integrator(Grid,dum(1:irc),1,irc) 2304 2305 WRITE(55,'(1p,2e16.7)') q(iq),tphij 2306 ENDDO 2307 CLOSE(55) 2308 ENDIF !(e<=0) 2309 ENDDO ! ib 2310 DEALLOCATE(q,dum,dum1) 2311 END SUBROUTINE ftprod 2312 !*********************************************************************** 2313 ! program to calculate Fourier transform of hatpot functions 2314 ! in order to get an idea of the convergence 2315 ! and output them to files hatpot.l 2316 !*********************************************************************** 2317 SUBROUTINE fthatpot(Grid) 2318 TYPE(GridInfo), INTENT(IN) :: Grid 2319 2320 INTEGER, PARAMETER :: nq=200 2321 REAL(8), PARAMETER :: qmax=15.d0 2322 2323 REAL(8), ALLOCATABLE :: q(:),dum(:),dum1(:),dum2(:),dum3(:),arg(:) 2324 REAL(8), POINTER :: r(:) 2325 REAL(8) :: h,dq,fthatp,fthatd,fthattd 2326 INTEGER :: i,ib,l,iq,n,irc,ll 2327 CHARACTER(4) flnm 2328 ! 2329 WRITE(STD_OUT,*) 'calculating Fourier transforms of hatpot for each l ' 2330 2331 n=Grid%n 2332 h=Grid%h 2333 r=>Grid%r 2334 irc=PAW%irc 2335 ALLOCATE(q(nq),dum(n),dum1(n),dum2(n),dum3(n),arg(n),stat=i) 2336 IF (i/=0) THEN 2337 WRITE(STD_OUT,*) 'Error in allocating space in fthatpot',n,nq,i 2338 STOP 2339 ENDIF 2340 2341 dq=qmax/nq 2342 ll=2*MAXVAL(PAW%l(:)) 2343 DO l=0,ll 2344 CALL mkname(l,flnm) 2345 OPEN(55,file='fthatpot.'//TRIM(flnm),form='formatted') 2346 CALL hatL(Grid,PAW,l,dum2) 2347 DO iq=1,nq 2348 q(iq)=iq*dq 2349 DO i=1,n 2350 dum(i)=q(iq)*r(i) 2351 ENDDO 2352 CALL sphbes(l,n,dum) 2353 DO i=1,n 2354 dum3(i)=dum(i)*dum2(i) 2355 ENDDO 2356 fthatd=integrator(Grid,dum3(1:n)) 2357 fthatp=8*PI*fthatd/q(iq)**2 2358 2359 fthattd=0 2360 IF(l==0) THEN 2361 arg(1:n)=dum(1:n)*PAW%tden(1:n) 2362 fthattd=integrator(Grid,arg) 2363 ENDIF 2364 2365 WRITE(55,'(1p,5e16.7)') q(iq),fthatp,fthatd,fthatp*fthatd,fthatp*fthattd 2366 ENDDO 2367 CLOSE(55) 2368 ENDDO ! l 2369 DEALLOCATE(q,dum,dum1,dum2,dum3,arg) 2370 END SUBROUTINE fthatpot 2371 2372 !*********************************************************************** 2373 ! program to calculate Fourier transform of tphi*q^2 2374 ! in order to get an idea of the convergence of kinetic energy 2375 ! and output them to files ftkin.l 2376 !*********************************************************************** 2377 SUBROUTINE ftkin(Grid) 2378 TYPE(GridInfo), INTENT(IN) :: Grid 2379 2380 INTEGER, PARAMETER :: nq=200 2381 REAL(8), PARAMETER :: qmax=15.d0 2382 2383 REAL(8), ALLOCATABLE :: q(:),dum(:),dum1(:) 2384 REAL(8), POINTER :: r(:) 2385 REAL(8) :: h,dq,kin 2386 INTEGER :: i,ib,l,iq,n,irc 2387 CHARACTER(4) flnm 2388 ! 2389 WRITE(STD_OUT,*) 'calculating Fourier transforms of tphi ',& 2390& 'For bound states only ' 2391 2392 n=Grid%n 2393 h=Grid%h 2394 r=>Grid%r 2395 ALLOCATE(q(nq),dum(n),dum1(n),stat=i) 2396 IF (i/=0) THEN 2397 WRITE(STD_OUT,*) 'Error in allocating space in ftkin',n,nq,i 2398 STOP 2399 ENDIF 2400 2401 dq=qmax/nq 2402 DO ib=1,PAW%nbase 2403 IF (PAW%eig(ib)<=0.d0) THEN 2404 l=PAW%l(ib) 2405 CALL mkname(ib,flnm) 2406 OPEN(55,file='ftkin.'//TRIM(flnm),form='formatted') 2407 DO iq=1,nq 2408 q(iq)=iq*dq 2409 DO i=1,n 2410 dum(i)=q(iq)*r(i) 2411 ENDDO 2412 CALL sphbes(l,n,dum) 2413 DO i=1,n 2414 dum1(i)=dum(i)*r(i)*PAW%tphi(i,ib) 2415 ENDDO 2416 kin=integrator(Grid,dum1(1:n))*(q(iq)) 2417 kin=kin*kin 2418 2419 WRITE(55,'(1p,2e16.7)') q(iq),kin 2420 ENDDO 2421 CLOSE(55) 2422 ENDIF !(e<=0) 2423 ENDDO ! ib 2424 DEALLOCATE(q,dum,dum1) 2425 END SUBROUTINE ftkin 2426 2427 !*********************************************************************** 2428 ! program to calculate Fourier transform of vloc and tden 2429 ! in order to get an idea of their convergence 2430 ! and output them to files ftvloc 2431 !*********************************************************************** 2432 SUBROUTINE ftvloc(Grid) 2433 TYPE(GridInfo), INTENT(IN) :: Grid 2434 2435 INTEGER, PARAMETER :: nq=200 2436 REAL(8), PARAMETER :: qmax=15.d0 2437 2438 REAL(8), ALLOCATABLE :: q(:),dum(:),dum1(:) 2439 REAL(8), POINTER :: r(:) 2440 REAL(8) :: h,dq,vloc,tden 2441 INTEGER :: i,ib,l,iq,n,irc 2442 CHARACTER(4) flnm 2443 ! 2444 WRITE(STD_OUT,*) 'calculating Fourier transforms of vloc and tden' 2445 2446 n=Grid%n 2447 h=Grid%h 2448 irc=PAW%irc 2449 r=>Grid%r 2450 ALLOCATE(q(nq),dum(n),dum1(n),stat=i) 2451 IF (i/=0) THEN 2452 WRITE(STD_OUT,*) 'Error in allocating space in ftvloc',n,nq,i 2453 STOP 2454 ENDIF 2455 2456 OPEN(55,file='ftvloc',form='formatted') 2457 dq=qmax/nq 2458 l=0 2459 DO iq=1,nq 2460 q(iq)=iq*dq 2461 DO i=1,n 2462 dum(i)=q(iq)*r(i) 2463 ENDDO 2464 CALL sphbes(l,n,dum) 2465 DO i=1,n 2466 dum1(i)=dum(i)*PAW%vloc(i)*(r(i)**2) 2467 dum(i)=dum(i)*PAW%tden(i) 2468 ENDDO 2469 vloc=integrator(Grid,dum1(1:n),1,irc) 2470 tden=integrator(Grid,dum(1:n)) 2471 2472 WRITE(55,'(1p,4e16.7)') q(iq),vloc,tden,vloc*tden 2473 ENDDO 2474 CLOSE(55) 2475 2476 DEALLOCATE(q,dum,dum1) 2477 END SUBROUTINE ftvloc 2478 2479 2480 !*************************************************************************** 2481 ! pgm to solve separable radial schroedinger equation 2482 ! at energy 'energy' and at angular momentum l 2483 ! 2484 ! with smooth potential rveff/r, given in uniform mesh of n points 2485 ! r=i*h, i=1,...n-1 ;assuming p(r)=C*r**(l+1)*polynomial(r) for r==0; 2486 ! p((n+1)*h)=0 2487 ! 2488 ! uses Noumerov algorithm 2489 ! 2490 ! For l=0,1 corrections are needed to approximate wfn(r=0) 2491 ! These depend upon: 2492 ! e0 (current guess of energy eigenvalue) 2493 ! l,nz==0 2494 ! v(0) == v0 electronic potential at r=0 2495 ! v'(0) == v0p derivative of electronic potential at r=0 2496 ! 2497 ! also returns node == number of nodes for calculated state 2498 ! 2499 ! proj == projector functions 2500 ! hij and qij == hamiltonianian and overlap matrix elements 2501 !*************************************************************************** 2502 SUBROUTINE unboundsep(Grid,Pot,PAW,nr,l,energy,wfn,node) 2503 TYPE(GridInfo), INTENT(IN) :: Grid 2504 TYPE(PotentialInfo), INTENT(IN) :: Pot 2505 TYPE(PseudoInfo), INTENT(IN) :: PAW 2506 INTEGER, INTENT(IN) :: l,nr 2507 REAL(8), INTENT(IN) :: energy 2508 REAL(8), INTENT(INOUT) :: wfn(:) 2509 INTEGER, INTENT(INOUT) :: node 2510 2511 INTEGER :: n,ia,ib,ic,nbase,icount,jcount,lcount,irc 2512 REAL(8) :: summ,h,scale,zeroval 2513 REAL(8), ALLOCATABLE :: y(:,:),b(:),a(:,:) 2514 2515 n=Grid%n; h=Grid%h; nbase=PAW%nbase; irc=PAW%irc 2516 !write(std_out,*) 'Entering unboundsep with l energy = ', l, energy 2517 2518 IF (nr<irc) THEN 2519 WRITE(STD_OUT,*) 'error in unboundsep -- nr < irc', nr,irc 2520 STOP 2521 ENDIF 2522 ! 2523 ! initialize wfn 2524 wfn=0 2525 wfn(2)=wfninit(0.d0,l,Pot%v0,Pot%v0p,energy,Grid%r(2)) 2526 zeroval=0 2527 if (l==1) zeroval=2 2528 2529 call forward_numerov(Grid,l,nr,energy,Pot%rv,zeroval,wfn,node) 2530 2531 ALLOCATE(y(nr,nbase),b(nbase),stat=ib) 2532 IF (ib/=0) THEN 2533 WRITE(STD_OUT,*) 'Error in unboundsep allocation', nr,nbase,ib 2534 STOP 2535 ENDIF 2536 2537 lcount=0 2538 DO ib=1,nbase 2539 IF (l==PAW%l(ib)) THEN 2540 lcount=lcount+1 2541 CALL inhomogeneous_numerov(Grid,l,nr,energy,& 2542& PAW%rveff,PAW%otp(:,ib),y(:,lcount)) 2543 ENDIF 2544 ENDDO 2545 2546 IF(lcount>0) THEN 2547 ALLOCATE(a(lcount,lcount),stat=ib) 2548 IF (ib/=0) THEN 2549 WRITE(STD_OUT,*) 'Error in unboundsep allocation', lcount,ib 2550 STOP 2551 ENDIF 2552 2553 icount=0 2554 DO ib=1,nbase 2555 summ=0.d0 2556 IF (l==PAW%l(ib)) THEN 2557 icount=icount+1 2558 DO ic=1,nbase 2559 IF (l==PAW%l(ic)) THEN 2560 summ=summ+(PAW%dij(ib,ic)-energy*PAW%oij(ib,ic))*& 2561& overlap(Grid,PAW%otp(:,ic),wfn,1,irc) 2562 ENDIF 2563 ENDDO 2564 b(icount)=-summ 2565 ENDIF 2566 ENDDO 2567 ! 2568 icount=0 2569 DO ia=1,nbase 2570 IF (l==PAW%l(ia)) THEN 2571 icount=icount+1 2572 jcount=0 2573 DO ib=1,nbase 2574 IF (l==PAW%l(ib)) THEN 2575 jcount=jcount+1 2576 summ=0.d0 2577 IF (ia.EQ.ib) summ=1.d0 2578 DO ic=1,nbase 2579 IF (l==PAW%l(ic)) THEN 2580 summ=summ+(PAW%dij(ia,ic)-energy*PAW%oij(ia,ic))*& 2581& overlap(Grid,PAW%otp(:,ic),y(:,jcount),1,irc) 2582 ENDIF 2583 ENDDO 2584 a(icount,jcount)=summ 2585 ENDIF 2586 ENDDO 2587 ENDIF 2588 ENDDO 2589 ! 2590 CALL linsol(a,b,lcount,lcount,lcount,nbase) 2591 2592 icount=0 2593 DO ib=1,nbase 2594 IF(l==PAw%l(ib)) THEN 2595 icount=icount+1 2596 wfn(1:nr)=wfn(1:nr)+b(icount)*y(1:nr,icount) 2597 ENDIF 2598 ENDDO 2599 DEALLOCATE(a) 2600 ENDIF 2601 ! 2602 ! normalize to unity within integration range 2603 ! 2604 scale=1.d0/sepnorm(Grid,PAW,nr,l,wfn) 2605 IF (scale.LE.0.d0) THEN 2606 WRITE(STD_OUT,*) 'warning -- negative norm for l=',l 2607 scale=-scale 2608 IF (scale.EQ.0.d0) scale=1.d0 2609 ENDIF 2610 scale=DSIGN(SQRT(scale),wfn(nr-2)) 2611 wfn(1:nr)=wfn(1:nr)*scale 2612 DEALLOCATE(b,y) 2613 END SUBROUTINE unboundsep 2614 2615 !*************************************************************************** 2616 ! pgm to solve separable radial schroedinger equation 2617 ! for bound state near energy 'energy' and at angular momentum l 2618 ! 2619 ! with smooth potential rveff/r, given in uniform mesh of n points 2620 ! r=i*h, i=1,...n-1 ;assuming p(r)=C*r**(l+1)*polynomial(r) for r==0; 2621 ! p((n+1)*h)=0 2622 ! 2623 ! uses Noumerov algorithm 2624 ! 2625 ! For l=0,1 corrections are needed to approximate wfn(r=0) 2626 ! These depend upon: 2627 ! e0 (current guess of energy eigenvalue) 2628 ! l,nz==0 2629 ! v(0) == v0 electronic potential at r=0 2630 ! v'(0) == v0p derivative of electronic potential at r=0 2631 ! 2632 ! also returns node == number of nodes for calculated state 2633 ! 2634 ! proj == projector functions 2635 ! hij and qij == hamiltonianian and overlap matrix elements 2636 !*************************************************************************** 2637 SUBROUTINE boundsep(Grid,Pot,PAW,l,node,energy,emin,emax,wfn) 2638 TYPE(GridInfo), INTENT(IN) :: Grid 2639 TYPE(PotentialInfo), INTENT(IN) :: Pot 2640 TYPE(PseudoInfo), INTENT(IN) :: PAW 2641 INTEGER, INTENT(IN) :: l,node 2642 REAL(8), INTENT(INOUT) :: energy,emin,emax 2643 REAL(8), INTENT(INOUT) :: wfn(:) 2644 2645 INTEGER, PARAMETER :: niter=100 2646 INTEGER :: n,i,j,ia,ib,ic,nbase,icount,jcount,lcount 2647 INTEGER :: match,mmatch,irc,node1,iter 2648 REAL(8), PARAMETER :: convre=1.d-10,err=1.d-9 2649 REAL(8) :: summ,h,scale,best,energy0,dele,x,rin,rout,zeroval,ppp 2650 REAL(8), ALLOCATABLE :: p1(:),u(:),y(:,:),b(:),a(:,:) 2651 2652 n=Grid%n; h=Grid%h; nbase=PAW%nbase; irc=PAW%irc 2653 ! 2654 ALLOCATE(p1(n),u(n),y(n,nbase),b(nbase),stat=ib) 2655 IF (ib/=0) THEN 2656 WRITE(STD_OUT,*) 'Error in boundsep allocation', n,nbase,ib 2657 STOP 2658 ENDIF 2659 !WRITE(STD_OUT,*) 'in boundsep with ', l,node,energy,emin,emax 2660 lcount=0 2661 DO ib=1,nbase 2662 IF (l==PAW%l(ib)) THEN 2663 lcount=lcount+1 2664 ENDIF 2665 ENDDO 2666 ! 2667 IF(lcount>0) THEN 2668 ALLOCATE(a(lcount,lcount),stat=ib) 2669 IF (ib/=0) THEN 2670 WRITE(STD_OUT,*) 'Error in boundsep allocation', lcount,ib 2671 STOP 2672 ENDIF 2673 ENDIF 2674 2675 IF (emax.GT.0.d0) emax=0.d0 2676 best=1.d10 2677 IF (energy.LT.emin) energy=emin+err 2678 IF (energy.GT.emax) energy=emax 2679 energy0=energy 2680 2681 iter=0 2682 2683 2684 DO 2685 match=MIN(irc+10,n-10) 2686 x=0.5d0*(Pot%rv(n)/Grid%r(n)+Pot%rv(n-1)/Grid%r(n-1))+& 2687& l*(l+1)/(Grid%r(n)**2) 2688 ppp=SQRT(ABS(x-energy)) 2689 p1=0 2690 p1(n)=1 2691 p1(n-1)=exp(-ppp*(Grid%r(n-1)-Grid%r(n))) 2692 2693 !write(std_out,*) 'before backward', n,p1(n-1),p1(n) 2694 !write(std_out,*) 'before backward', x,ppp,exp(-ppp*(Grid%r(n-1)-Grid%r(n))) 2695 !write(std_out,*) 'x,energy', x,energy,ABS(x-energy),SQRT(ABS(x-energy)) 2696 !call flush_unit(std_out) 2697 2698 CALL backward_numerov(Grid,l,match-5,energy,Pot%rv,p1) 2699 rin=Gfirstderiv(Grid,match,p1)/p1(match) 2700 mmatch=match+1 2701 !WRITE(STD_OUT,*) 'match, rin ' ,match,rin 2702 ! 2703 ! perform outward integration until match point -- it is assumed 2704 ! that projector functions proj are zero for r>r(match) 2705 ! 2706 ! initialize u 2707 u=0 2708 u(2)=wfninit(0.d0,l,Pot%v0,Pot%v0p,energy,Grid%r(2)) 2709 zeroval=0 2710 if (l==1) zeroval=2 2711 2712 call forward_numerov(Grid,l,mmatch+5,energy& 2713& ,Pot%rv,zeroval,u,node1) 2714 ! 2715 lcount=0 2716 DO ib=1,nbase 2717 IF (l==PAW%l(ib)) THEN 2718 lcount=lcount+1 2719 CALL inhomogeneous_numerov(Grid,l,mmatch+5,energy,& 2720& Pot%rv,PAW%otp(:,ib),y(:,lcount)) 2721 ENDIF 2722 ENDDO 2723 ! 2724 IF(lcount>0) THEN 2725 2726 icount=0 2727 DO ib=1,nbase 2728 summ=0.d0 2729 IF (l==PAW%l(ib)) THEN 2730 icount=icount+1 2731 DO ic=1,nbase 2732 IF (l==PAW%l(ic)) THEN 2733 summ=summ+(PAW%dij(ib,ic)-energy*PAW%oij(ib,ic))*& 2734& overlap(Grid,PAW%otp(:,ic),u,1,irc) 2735 ENDIF 2736 ENDDO 2737 b(icount)=-summ 2738 ENDIF 2739 ENDDO 2740 ! 2741 icount=0 2742 DO ia=1,nbase 2743 IF (l==PAW%l(ia)) THEN 2744 icount=icount+1 2745 jcount=0 2746 DO ib=1,nbase 2747 IF (l==PAW%l(ib)) THEN 2748 jcount=jcount+1 2749 summ=0.d0 2750 IF (ia.EQ.ib) summ=1.d0 2751 DO ic=1,nbase 2752 IF (l==PAW%l(ic)) THEN 2753 summ=summ+(PAW%dij(ia,ic)-energy*PAW%oij(ia,ic))*& 2754& overlap(Grid,PAW%otp(:,ic),y(:,jcount),1,irc) 2755 ENDIF 2756 ENDDO 2757 a(icount,jcount)=summ 2758 ENDIF 2759 ENDDO 2760 ENDIF 2761 ENDDO 2762 ! 2763 CALL linsol(a,b,lcount,lcount,lcount,nbase) 2764 2765 wfn=0 2766 wfn(1:mmatch+5)=u(1:mmatch+5) 2767 icount=0 2768 DO ib=1,nbase 2769 IF(l==PAw%l(ib)) THEN 2770 icount=icount+1 2771 wfn(1:mmatch+5)=wfn(1:mmatch+5)+b(icount)*y(1:mmatch+5,icount) 2772 ENDIF 2773 ENDDO 2774 ENDIF 2775 2776 rout=Gfirstderiv(Grid,match,wfn)/wfn(match) 2777 !WRITE(STD_OUT,'("node,match,rin,rout",2i8,1p,2e15.7)') node1,match,rin,rout 2778 ! -- estimate correction 2779 node1=0 2780 wfn(:)=wfn(:)/wfn(match) 2781 DO j=3,match 2782 IF (wfn(j)*wfn(j-1).LT.0.d0) node1=node1+1 2783 ENDDO 2784 !WRITE(STD_OUT,*) 'actual number of nodes', node1 2785 2786!This test is obsolete: pseudo-WFs do not have to be orthogonal 2787! IF (node1<node) THEN 2788! ! too few nodes -- raise energy 2789! emin=MAX(energy+err,emin) 2790! energy=emax-(emax-energy)*ranx() 2791! !WRITE(STD_OUT,*) 'too few nodes -- energy raised', energy,emin,emax 2792! ELSEIF (node1>node) THEN 2793! ! too many nodes -- lower energy 2794! emax=MIN(energy-err,emax) 2795! energy=emin+(energy-emin)*ranx() 2796! !WRITE(STD_OUT,*) 'too many nodes -- energy lowered', energy,emin,emax 2797! !do i=1,mmatch 2798! !write(200+iter,'(1p,7e15.7)')Grid%r(i),wfn(i) 2799! !enddo 2800! ELSEIF (node1==node) THEN 2801 DO j=match,n 2802 wfn(j)=p1(j)/p1(match) 2803 ENDDO 2804 ! normalization 2805 scale=1.d0/sepnorm(Grid,PAW,n,l,wfn) 2806 dele=(rout-rin)*scale 2807 !WRITE(STD_OUT,*) 'dele,scale',dele,scale 2808 scale=SQRT(scale) 2809 wfn=scale*wfn 2810 2811 !do i=1,n 2812 ! write(100+iter,'(1p,2E15.7)') Grid%r(i),wfn(i) 2813 !enddo 2814 2815 x=ABS(dele) 2816 IF (x.LT.best) THEN 2817 energy0=energy 2818 best=x 2819 ENDIF 2820 IF (ABS(dele).LE.convre) EXIT 2821 energy=energy+dele 2822 !WRITE(STD_OUT,*) 'next energy' , energy,dele 2823 ! if energy is out of range, pick random energy in correct range 2824 IF (emin-energy.GT.convre.OR.energy-emax.GT.convre) THEN 2825 energy=emin+(emax-emin)*ranx()+err 2826 ! WRITE(STD_OUT,*) 'energy out of range -- reranged --', energy 2827 ENDIF 2828! ENDIF 2829 iter=iter+1 2830 !WRITE(STD_OUT,*) 'Energy for next iteration ', iter,energy 2831 IF (iter.GT.niter) THEN 2832 WRITE(STD_OUT,*) 'no convergence in boundsep',l,dele,energy 2833 WRITE(STD_OUT,*) ' best guess of eig, dele = ',energy0,best 2834 STOP 2835 ENDIF 2836 ENDDO 2837 ! 2838 ! normalize to unity within integration range 2839 ! 2840 CALL filter(n,wfn,1.d-11) 2841 scale=1.d0/sepnorm(Grid,PAW,n,l,wfn) 2842 IF (scale.LE.0.d0) THEN 2843 WRITE(STD_OUT,*) 'warning -- negative norm for l=',l 2844 scale=-scale 2845 IF (scale.EQ.0.d0) scale=1.d0 2846 ENDIF 2847 scale=DSIGN(SQRT(scale),wfn(n-2)) 2848 wfn(1:n)=wfn(1:n)*scale 2849 !WRITE(STD_OUT,*) 'exiting boundsep with energy ', l,energy 2850 DEALLOCATE(a,b,y,u,p1) 2851 END SUBROUTINE boundsep 2852 2853 !********************************************************* 2854 ! program to to transform smooth wavefunction to all-electron 2855 ! wavefunction within Blochl's paw formalism 2856 ! otp == projector function 2857 ! odphi == ophi - otphi 2858 !********************************************************* 2859 SUBROUTINE PStoAE(Grid,PAW,nr,l,tpsi,psi) 2860 TYPE(GridInfo), INTENT(IN) :: Grid 2861 TYPE(PseudoInfo), INTENT(IN) :: PAW 2862 INTEGER, INTENT(IN) :: nr,l 2863 REAL(8), INTENT(IN) :: tpsi(:) 2864 REAL(8), INTENT(INOUT) :: psi(:) 2865 2866 INTEGER :: n,ib,nbase,irc 2867 REAL(8) :: h,scale 2868 2869 n=Grid%n; h=Grid%h; nbase=PAW%nbase; irc=PAW%irc 2870 2871 IF (nr<irc) THEN 2872 WRITE(STD_OUT,*) 'error in PStoAE -- nr < irc', nr,irc 2873 STOP 2874 ENDIF 2875 2876 psi(1:nr)=tpsi(1:nr) 2877 DO ib=1,nbase 2878 IF (l==PAW%l(ib)) psi(1:nr)=psi(1:nr)+& 2879& (PAW%ophi(1:nr,ib)-PAW%otphi(1:nr,ib))*& 2880& overlap(Grid,PAW%otp(:,ib),tpsi,1,irc) 2881 ENDDO 2882 END SUBROUTINE PStoAE 2883 2884 SUBROUTINE Set_PAW_MatrixElements(Grid,PAW,ifen) 2885 ! Calculate PAW matrix elements from reference configuration before 2886 ! unscreening/rescreening 2887 TYPE(GridInfo), INTENT(IN) :: Grid 2888 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 2889 INTEGER, INTENT(IN) :: ifen 2890 2891 INTEGER :: nbase,l,ib,ic,io 2892 REAL(8) :: x 2893 TYPE(OrbitInfo), POINTER :: PSO 2894 REAL(8), allocatable :: wij(:,:) 2895 2896 PAW%oij=0 2897 PAW%dij=0 2898 PAW%wij=0 2899 nbase=PAW%nbase 2900 PSO=>PAW%TOCCWFN 2901 ALLOCATE (wij(nbase,nbase)) 2902 2903 do ib=1,nbase 2904 do ic=1,nbase 2905 IF (PAW%l(ib)==PAW%l(ic)) THEN 2906 CALL dqij(Grid,PAW,ib,ic,PAW%oij(ib,ic)) 2907 CALL altdtij(Grid,PAW,ib,ic,PAW%dij(ib,ic)) 2908 CALL avij(Grid,PAW,ib,ic,x) 2909 PAW%dij(ib,ic)=PAW%dij(ib,ic)+x 2910 write(std_out,'(2i5, 1p,3e17.8)') ib,ic,PAW%oij(ib,ic),PAW%dij(ib,ic),x 2911 Endif 2912 enddo 2913 enddo 2914 2915 wij=0 2916 Do io=1,PSO%norbit 2917 l=PSO%l(io) 2918 if (PSO%occ(io)>1.d-8.and..NOT.PSO%iscore(io)) then 2919 CALL calcwij(Grid,PAW,l,PSO%occ(io),PSO%wfn(:,io),wij) 2920 endif 2921 enddo 2922 2923 do ib=1,nbase 2924 do ic=1,nbase 2925 PAW%wij(ib,ic)=wij(ib,ic) 2926 enddo 2927 enddo 2928 2929 DEALLOCATE(wij) 2930 2931 if (.not.Check_overlap_of_projectors(Grid,PAW,ifen)) then 2932 write(std_out,*) "The overlap operator has at leat one negative eigenvalue!" 2933 write(std_out,*) "It might be not positive definite..." 2934 write(std_out,*) "Program is stopping." 2935 write(std_out,*) "This probably means that your projectors are too similar" 2936 write(std_out,*) "or that your PAW basis is incomplete." 2937 write(std_out,*) "Advice: try to change your input parameters (f.i. : reference energies)." 2938 stop 2939 end if 2940 2941 END SUBROUTINE Set_PAW_MatrixElements 2942 2943 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2944 ! FUNCTION Check_overlap_of_projectors(Grid,PAW,ifen) 2945 ! function written by Francois Jollet and Marc Torrent 2946 ! Approximate assessment of "completeness" of PAW basis set 2947 ! measured by the eigenvalues of the 1+O operator int 2948 ! the "basis" of the projectors May 2019 2949 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2950 FUNCTION Check_overlap_of_projectors(Grid,PAW,ifen) 2951 ! Check that the overlap operator is positive definite 2952 ! in the "basis" of projectors 2953 LOGICAL :: Check_overlap_of_projectors 2954 TYPE(GridInfo), INTENT(IN) :: Grid 2955 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 2956 INTEGER, INTENT(IN) :: ifen 2957 2958 REAL(8),PARAMETER :: tol=1.d-10 2959 2960 INTEGER :: n,nbase,i,j,k,m,info 2961 REAL(8) :: ovlp_ik,ovlp_jm 2962 REAL(8),allocatable :: ovlp(:,:),dum(:) 2963! REAL(8),allocatable :: wr(:),wi(:),work(:),vl(:,:),vr(:,:) 2964 REAL(8), allocatable :: eig(:),work(:) 2965 2966 Check_overlap_of_projectors=.true. 2967 2968 n=PAW%irc ; nbase=PAW%nbase 2969 allocate(ovlp(nbase,nbase)) 2970 ovlp(:,:)=0.d0 2971 2972 allocate(dum(n)) 2973 do k=1,nbase 2974 do m=1,nbase 2975 if (PAW%l(k)==PAW%l(m)) then 2976 dum(1:n)=PAW%otp(1:n,k)*PAW%otp(1:n,m) 2977 ovlp(m,k)=integrator(Grid,dum(1:n),1,n) 2978 2979 do i=1,nbase 2980 if (PAW%l(i)==PAW%l(k)) then 2981 dum(1:n)=PAW%otp(1:n,i)*PAW%otp(1:n,k) 2982 ovlp_ik=integrator(Grid,dum(1:n),1,n) 2983 2984 do j=1,nbase 2985 if (PAW%l(j)==PAW%l(m)) then 2986 dum(1:n)=PAW%otp(1:n,j)*PAW%otp(1:n,m) 2987 ovlp_jm=integrator(Grid,dum(1:n),1,n) 2988 2989 ovlp(m,k)=ovlp(m,k)+ovlp_ik*PAW%oij(i,j)*ovlp_jm 2990 2991 end if ! l_j=l_m 2992 end do ! Loop j 2993 end if ! l_i=l_k 2994 end do ! Loop i 2995 end if ! l_k=l_m 2996 end do ! Loop m 2997 end do ! Loop k 2998 deallocate(dum) 2999 3000! allocate(wr(nbase),wi(nbase),work(4*nbase)) 3001! allocate(vl(nbase,nbase),vr(nbase,nbase)) 3002! call dgeev('N','N',nbase,ovlp,nbase,wr,wi,vl,nbase,vr,nbase,work,4*nbase,info) 3003 3004 allocate(eig(nbase),work(4*nbase)) 3005 call DSYEV('N','U',nbase,ovlp,nbase,eig,work,4*nbase,info) 3006 write(std_out,'(" Completed diagonalization of ovlp with info = ", i8)')info 3007 write(ifen,'(" Completed diagonalization of ovlp with info = ", i8)')info 3008 if (info==0) then 3009 write(std_out,*) " " 3010 write(std_out,*) "Eigenvalues of overlap operator (in the basis of projectors):" 3011 write(ifen,*) " " 3012 write(ifen,*) "Eigenvalues of overlap operator (in the basis of projectors):" 3013 do i=1,nbase 3014 write(std_out,'( i5,5x,1p,1e17.8)') i,eig(i) 3015 write(ifen,'( i5,5x,1p,1e17.8)') i,eig(i) 3016 if (eig(i)<tol) Check_overlap_of_projectors=.false. 3017 enddo 3018 else 3019 write(std_out,*) 'Stopping due to failure of ovlp diagonalization' 3020 write(ifen,*) 'Stopping due to failure of ovlp diagonalization' 3021 stop 3022 endif 3023 write(std_out,*) " " 3024 write(ifen,*) " " 3025 3026! do i=1,nbase 3027! write(std_out,*) i,wr(i) 3028! if (wr(i)<tol) Check_overlap_of_projectors=.false. 3029! end do 3030 3031 deallocate(eig,work) 3032 deallocate(ovlp) 3033 3034 END FUNCTION Check_overlap_of_projectors 3035 3036 !************************************************************************ 3037 ! program to calculate logerivatives of paw wavefunctions 3038 ! and to compare them with all electron wavefunctions 3039 ! optionally, wavefunctions are written to a file 3040 ! Assumes prior call to SUBROUTINE Set_PAW_MatrixElements(Grid,PAW) 3041 ! 7/2018 program modified by Casey Brock to output atan(logderiv) 3042 ! values, with pi jumps taken into account 3043 ! 3044 !************************************************************************ 3045 SUBROUTINE logderiv(Grid,Pot,PAW) 3046 TYPE(GridInfo), INTENT(IN) :: Grid 3047 TYPE(PotentialInfo), INTENT(IN) :: Pot 3048 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 3049 3050 TYPE(PotentialInfo) :: PS 3051 INTEGER :: n,l,ie,ke,nbase,ib,ic,nr,i,nodes,mbase,irc,lng 3052 REAL(8), PARAMETER :: e0=-5.d0 ! Where to print WFn in file 3053 REAL(8) :: de,h,x,dwdr,dcwdr,scale,energy 3054 REAL(8), ALLOCATABLE :: psi(:),tpsi(:),ttpsi(:) 3055 CHARACTER(4) :: flnm 3056 REAL(8) :: y_ae, x_ae, y_ps, x_ps 3057 REAL(8) :: cumshift_ae, cumshift_ps 3058 REAL(8) :: atan_ae, atan_ps, atan_prev_ae, atan_prev_ps 3059 REAL(8) :: slope_ae, slope_ps, slope_prev_ae, slope_prev_ps 3060 REAL(8) :: atan_shifted_ae, atan_shifted_ps 3061 INTEGER :: cnt_ae, cnt_ps 3062 3063 3064 n=Grid%n; h=Grid%h; nbase=PAW%nbase; irc=PAW%irc; nr=irc+10 3065 3066 ALLOCATE(psi(nr),tpsi(nr),ttpsi(nr),PS%rv(n),stat=ie) 3067 IF (ie/=0) THEN 3068 WRITE(STD_OUT,*) 'Error in logderiv allocation',n,ie 3069 STOP 3070 ENDIF 3071 3072 ! load PS 3073 PS%rv=PAW%rveff ; PS%nz=0.d0 3074 call zeropot(Grid,PS%rv,PS%v0,PS%v0p) 3075 ! 3076 ! calculate logderivatives at irc 3077 ! 3078 WRITE(STD_OUT,*) 'calculating log derivatives at irc',Grid%r(irc) 3079 ! 3080 de=(maxlogderiv-minlogderiv)/dble(nlogderiv-1) 3081 ke=1+anint((e0-minlogderiv)/de) 3082 3083 mbase=nbase 3084 DO l=0,PAW%lmax+1 3085 CALL mkname(l,flnm) 3086 OPEN(56,file='logderiv.'//TRIM(flnm),form='formatted') 3087 3088! initialize variables for atan calculation 3089 cnt_ae = 1 3090 slope_ae = 0.d0 3091 cumshift_ae = 0.d0 3092 cnt_ps = 1 3093 slope_ps = 0.d0 3094 cumshift_ps = 0.d0 3095 3096 DO ie=1,nlogderiv 3097 energy=minlogderiv+de*(ie-1) 3098 psi=0;tpsi=0;ttpsi=0 3099 if (scalarrelativistic) then 3100 CALL unboundsr(Grid,Pot,nr,l,energy,psi,nodes) 3101 elseif (TRIM(PAW%exctype)=='HF') then 3102 CALL HFunocc(Grid,PAW%OCCWFN,l,energy,Pot%rv,Pot%v0,Pot%v0p,& 3103& psi,lng) 3104 else 3105 CALL unboundsch(Grid,Pot%rv,Pot%v0,Pot%v0p,nr,l,energy,psi,nodes) 3106 endif 3107 CALL unboundsep(Grid,PS,PAW,nr,l,energy,tpsi,nodes) 3108 CALL PStoAE(Grid,PAW,nr,l,tpsi,ttpsi) 3109 ! 3110! dwdr=Gfirstderiv(Grid,irc,psi)/psi(irc) 3111! dcwdr=Gfirstderiv(Grid,irc,ttpsi)/ttpsi(irc) 3112 !!! CALCULATE ATANS 3113 y_ae = Gfirstderiv(Grid,irc,psi) 3114 x_ae = psi(irc) 3115 y_ps = Gfirstderiv(Grid,irc,ttpsi) 3116 x_ps = ttpsi(irc) 3117 3118 dwdr = y_ae/x_ae 3119 dcwdr = y_ps/x_ps 3120 3121! calculate atan and 3122! phase shift by n*pi if necessary so the atan curve is continous 3123 atan_ae = atan2(y_ae, x_ae) 3124 IF (ie>1) THEN 3125 slope_ae = atan_ae - atan_prev_ae 3126 CALL phase_unwrap(ie, y_ae, x_ae, atan_ae, atan_prev_ae, slope_ae, slope_prev_ae, cumshift_ae, cnt_ae) 3127 ENDIF 3128 atan_shifted_ae = atan_ae + cumshift_ae 3129 slope_prev_ae = slope_ae 3130 atan_prev_ae = atan_ae 3131 atan_ps = atan2(y_ps, x_ps) 3132 IF (ie>1) THEN 3133 slope_ps = atan_ps - atan_prev_ps 3134 CALL phase_unwrap(ie, y_ps, x_ps, atan_ps, atan_prev_ps, slope_ps, slope_prev_ps, cumshift_ps, cnt_ps) 3135 ENDIF 3136 atan_shifted_ps = atan_ps + cumshift_ps 3137 slope_prev_ps = slope_ps 3138 atan_prev_ps = atan_ps 3139 3140 WRITE(56,'(1p,5e12.4)') energy,dwdr,dcwdr,atan_shifted_ae,& 3141& atan_shifted_ps 3142 IF (ie.EQ.ke) THEN 3143 mbase=mbase+1 3144 CALL mkname(mbase,flnm) 3145 OPEN(57,file='wfn'//TRIM(flnm),form='formatted') 3146 WRITE(57,*) '# l=',l,'energy=',energy 3147 ! 3148 ! form converted wavefunction and rescale exact wavefunction 3149 ! 3150 scale=ttpsi(irc)/psi(irc) 3151 DO i=1,nr 3152 WRITE(57, '(1p,5e12.4)') Grid%r(i),tpsi(i),ttpsi(i),psi(i)*scale 3153 ENDDO 3154 CLOSE(57) 3155 ENDIF 3156 ENDDO !ie 3157 CLOSE(56) 3158 3159 ENDDO !l 3160 3161 DEALLOCATE(psi,tpsi,PS%rv) 3162 END SUBROUTINE logderiv 3163 3164 3165 ! Assumes prior call to SUBROUTINE Set_PAW_MatrixElements(Grid,PAW) 3166 SUBROUTINE FindVlocfromVeff(Grid,Orbit,PAW) 3167 TYPE(GridInfo), INTENT(INOUT) :: Grid 3168 TYPE(OrbitInfo), INTENT(IN) :: Orbit 3169 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 3170 3171 REAL(8), POINTER :: r(:) 3172 REAL(8) :: h,qeff,tq,rat,q00,ecoul,etxc,eexc,occ,fac 3173 INTEGER :: n,i,irc,io,nbase,ib,ic 3174 REAL(8), allocatable :: d(:),dv(:),dvx(:),v(:),vv(:),dpdr(:),pbr(:) 3175 3176 CALL FillHat(Grid,PAW) 3177 3178 if (Vlocalindex == SETVLOC) then 3179 write(std_out,*) 'Vloc == VlocCoef*shapefunc ' 3180 return 3181 endif 3182 3183 n=Grid%n;irc=PAW%irc 3184 nbase=PAW%nbase 3185 h=Grid%h ; r=>Grid%r 3186 irc=max(PAW%irc,PAW%irc_shap,PAW%irc_vloc,PAW%irc_core) 3187 3188 PAW%den=0.d0;PAW%tden=0.d0 3189 PAW%valetau=0.d0;PAW%tvaletau=0.d0 3190 allocate(pbr(n),dpdr(n),d(n)) 3191 do io=1,PAW%OCCWFN%norbit 3192 if (.not.PAW%OCCwfn%iscore(io)) then 3193 occ=PAW%OCCwfn%occ(io) 3194 fac=PAW%OCCwfn%l(io)*(PAW%OCCwfn%l(io)+1) 3195 PAW%den=PAW%den+occ*PAW%OCCwfn%wfn(:,io)**2 3196 PAW%tden=PAW%tden+occ*PAW%TOCCwfn%wfn(:,io)**2 3197 dpdr=0.d0;pbr=0.d0 3198 pbr(2:n)=PAW%OCCwfn%wfn(2:n,io)/Grid%r(2:n) 3199 CALL derivative(Grid,pbr,dpdr,2,Grid%n) 3200 CALL extrapolate(Grid,pbr) 3201 CALL extrapolate(Grid,dpdr) 3202 d(:)=((Grid%r(:)*dpdr(:)))**2+fac*(pbr(:))**2 3203 PAW%valetau=PAW%valetau+occ*d 3204 dpdr=0.d0;pbr=0.d0 3205 pbr(2:n)=PAW%TOCCwfn%wfn(2:n,io)/Grid%r(2:n) 3206 CALL derivative(Grid,pbr,dpdr,2,Grid%n) 3207 CALL extrapolate(Grid,pbr) 3208 CALL extrapolate(Grid,dpdr) 3209 d(:)=((Grid%r(:)*dpdr(:)))**2+fac*(pbr(:))**2 3210 PAW%tvaletau=PAW%tvaletau+occ*d 3211 endif 3212 enddo 3213 deallocate(pbr,dpdr,d) 3214 3215 allocate(d(n),dv(n),dvx(n),STAT=i) 3216 if (i/=0) stop 'Error (1) in allocating arrays in findvlocfromveff' 3217 3218 d=PAW%tden+PAW%tcore 3219 call poisson(Grid,tq,d,dv,ecoul) 3220 3221 d=PAW%core-PAW%tcore 3222 q00=0.5d0*PAW%AErefrv(1)+integrator(Grid,d) 3223 write(std_out,*) 'nucleus and core q00 ', q00 3224 do ib=1,nbase 3225 do ic=1,nbase 3226 q00=q00+PAW%wij(ib,ic)*PAW%oij(ib,ic) 3227 enddo 3228 enddo 3229 write(std_out,*) 'Complete q00 ', q00 3230 3231 dv=dv+q00*PAW%hatpot 3232 If (TRIM(PAW%exctype)=='EXX') THEN 3233 CALL EXX_pseudoVx(Grid,Orbit,PAW,dvx) 3234 PAW%trvx=dvx 3235 dv=dv+dvx 3236 ELSEIf (TRIM(PAW%exctype)=='EXXKLI') THEN 3237 write(std_out,*) 'before EXX_pseudo'; call flush_unit(std_out) 3238 CALL EXXKLI_pseudoVx(Grid,PAW,dvx) 3239 PAW%trvx=dvx 3240 dv=dv+dvx 3241 ELSEIf (TRIM(PAW%exctype)=='HF') THEN 3242 dvx=0.d0 3243 PAW%trvx=0.d0 3244 ELSE 3245 d=PAW%tden+PAW%tcore 3246 CALL exch(Grid,d,dvx,etxc,eexc) 3247 PAW%trvx=dvx 3248 dv=dv+dvx 3249 endif 3250 3251 rat=0.d0 3252 DO i=2,n 3253 PAW%vloc(i)=(PAW%rveff(i)-dv(i))/r(i) 3254 IF (i>=irc) rat=rat+ABS(PAW%vloc(i)) 3255 ENDDO 3256 WRITE(STD_OUT,*) 'Error in vloc -- ',rat 3257 call extrapolate(Grid,PAW%vloc) 3258 PAW%vloc(irc:n)=0.d0 3259 3260!!!!!!!!!!!This part does not work for HF!!!! 3261! Construct ionic local potential for abinit from screened pseudopotential 3262! in addition to ionic unscreening include hathat density 3263! in exchange-correlation functional 3264! also generate version without hathat density in vxc 3265 3266 PAW%abinitvloc=0.d0 3267 PAW%abinitnohat=0.d0 3268 allocate(v(n),vv(n),STAT=i) 3269 if (i/=0) stop 'Error (2) in allocating arrays in findvlocfromveff' 3270 3271 d=PAW%den-PAW%tden 3272 tq=integrator(Grid,d,1,irc) 3273 write(std_out,*) ' abinit tq = ', tq 3274 3275! Compute VH(tDEN+hatDEN) 3276 d=PAW%tden+tq*PAW%hatden 3277 CALL poisson(Grid,q00,d,v,rat) 3278 3279! Compute Vxc(tcore+tDEN) 3280 d=PAW%tcore+PAW%tden 3281 CALL exch(Grid,d,v,etxc,eexc) 3282 3283! Compute Vxc(tcore+tDEN+hatDEN) 3284 d=PAW%tcore+PAW%tden+tq*PAW%hatden 3285 CALL exch(Grid,d,vv,etxc,eexc) 3286 3287! Store Vxc(tcore+tDEN)-Vxc(tcore+tDEN+hatDEN) 3288 do i=2,n 3289 PAW%abinitvloc(i)=(v(i)-vv(i))/Grid%r(i) 3290 enddo 3291 call extrapolate(Grid,PAW%abinitvloc) 3292 3293! Compute VH(tcore+(Nc-tNc-Z)*hatDEN) 3294! (For consistency with ABINIT, need to integrate Poisson 3295! equation up to coretailpoints (not the whole mesh)) 3296 d=PAW%core-PAW%tcore 3297 qeff=0.5d0*PAW%AErefrv(1)+integrator(Grid,d,1,PAW%irc_core) 3298 d=PAW%tcore+qeff*PAW%hatden 3299 Grid%n=PAW%coretailpoints 3300 CALL poisson_marc(Grid,q00,d(1:PAW%coretailpoints),vv(1:PAW%coretailpoints),rat) 3301 Grid%n=n 3302 3303! Compute ionic potential following Bloechl formalism 3304 do i=2,PAW%coretailpoints 3305 PAW%abinitnohat(i)=PAW%vloc(i)+vv(i)/r(i) ! in Rydberg units 3306 enddo 3307 call extrapolate(Grid,PAW%abinitnohat) 3308 PAW%abinitnohat(PAW%coretailpoints+1:n)=PAW%vloc(PAW%coretailpoints+1:n) & 3309& +vv(PAW%coretailpoints)/Grid%r(PAW%coretailpoints+1:n) 3310 3311! Compute ionic potential following Kresse formalism 3312 do i=1,n 3313 PAW%abinitvloc(i)=PAW%abinitvloc(i)+PAW%abinitnohat(i) ! in Rydberg units 3314 enddo 3315 3316! check if PAW%tcore+tq*PAW%hatden is positive 3317 PAW%poscorenhat=.true. 3318 PAW%nhatv=tq*PAW%hatden 3319 do i=1,irc 3320 if ((PAW%tcore(i)+PAW%nhatv(i))<-1.d-8) then 3321 PAW%poscorenhat=.false. 3322 exit 3323 endif 3324 enddo 3325 3326 open(123,file='compare.abinit', form='formatted') 3327 do i=2,n 3328 write(123,'(1p,5e16.7)') r(i),PAW%rveff(i)/r(i),& 3329& PAW%vloc(i),PAW%abinitvloc(i),PAW%abinitnohat(i) 3330 enddo 3331 close(123) 3332 3333 deallocate(v,vv) 3334 deallocate(d,dv,dvx) 3335 3336 END SUBROUTINE FindVlocfromVeff 3337 3338!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3339! SCFPAW 3340! It is assumed that the new configuration involves only occupation 3341! changes 3342!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3343 SUBROUTINE SCFPAW(Grid,PAW) 3344 Type(GridInfo), INTENT(IN) :: Grid 3345 Type(PseudoInfo), INTENT(INOUT) :: PAW 3346 3347 INTEGER :: i,j,k,l,n,nvorbit,nbase,io,ib,jo,jb,ip,nfix 3348 INTEGER :: loop 3349 INTEGER :: mxloop=1000 3350 REAL(8) :: err0=1.d-7,mix=0.5d0 3351 TYPE(OrbitInfo), POINTER :: AEO,PSO 3352 INTEGER :: firsttime=0,norbit_mod,np(5) 3353 REAL(8) :: xocc,err 3354 LOGICAL :: success 3355 CHARACTER(4) :: stuff 3356 INTEGER, ALLOCATABLE :: tmap(:) 3357 INTEGER, ALLOCATABLE :: orbit_mod_l(:),orbit_mod_n(:),orbit_mod_k(:) 3358 REAL(8), ALLOCATABLE :: orbit_mod_occ(:) 3359 3360 AEO=>PAW%OCCWFN; PSO=>PAW%TOCCWFN 3361 WRITE(STD_OUT,*) 'Current occupancies:' 3362 IF (.NOT.diracrelativistic) THEN 3363 WRITE(STD_OUT,*) ' n l occupancy energy ' 3364 ELSE 3365 WRITE(STD_OUT,*) ' n l kap occupancy energy ' 3366 ENDIF 3367 DO io=1,PSO%norbit 3368 IF (.NOT.PSO%iscore(io)) THEN 3369 IF (.NOT.diracrelativistic) THEN 3370 WRITE(STD_OUT,'(i2,1x,i2,4x,1p,2e15.7)') PSO%np(io),& 3371& PSO%l(io),PSO%occ(io), PSO%eig(io) 3372 ELSE 3373 WRITE(STD_OUT,'(i2,1x,i2,2x,i2,4x,1p,2e15.7)') PSO%np(io),& 3374& PSO%l(io),PSO%kappa(io),PSO%occ(io), PSO%eig(io) 3375 ENDIF 3376 if (firsttime==0) then 3377 ib=PAW%valencemap(io) 3378 write(std_out,*) 'Setting pseudo orbital ', io,ib 3379 PSO%wfn(:,io)=PAW%tphi(:,ib) 3380 endif 3381 ENDIF 3382 ENDDO 3383 firsttime=firsttime+1 3384 3385 !Read modified occupations from standard input 3386 np(1)=AEO%nps;np(2)=AEO%npp;np(3)=AEO%npd;np(4)=AEO%npf;np(5)=AEO%npg 3387 CALL input_dataset_read_occ(norbit_mod,orbit_mod_l,orbit_mod_n,orbit_mod_k,& 3388& orbit_mod_occ,np,diracrelativistic) 3389 3390 DO jo=1,norbit_mod 3391 nfix=-100 3392 DO io=1,PSO%norbit 3393 IF (orbit_mod_n(jo)==PSO%np(io).AND.orbit_mod_l(jo)==PSO%l(io).AND.& 3394& ((.NOT.diracrelativistic).OR.orbit_mod_k(jo)==PSO%kappa(io))) THEN 3395 IF (PSO%iscore(io)) THEN 3396 write(std_out,*) 'Core orbitals cannot be changed',& 3397 & orbit_mod_n(jo),orbit_mod_l(jo),orbit_mod_occ(jo) 3398 stop 3399 endif 3400 nfix=io 3401 EXIT 3402 ENDIF 3403 ENDDO 3404 IF (nfix.LE.0) THEN 3405 WRITE(STD_OUT,*) 'error in occupations -- ip,l,xocc', & 3406& orbit_mod_n(jo),orbit_mod_l(jo),orbit_mod_occ(jo),nfix 3407 STOP 3408 ENDIF 3409 PSO%occ(nfix)=orbit_mod_occ(jo) 3410 AEO%occ(nfix)=orbit_mod_occ(jo) 3411 ENDDO 3412 DEALLOCATE(orbit_mod_l) 3413 DEALLOCATE(orbit_mod_n) 3414 DEALLOCATE(orbit_mod_k) 3415 DEALLOCATE(orbit_mod_occ) 3416 3417 WRITE(STD_OUT,*) 'New configuration:' 3418 DO io=1,PSO%norbit 3419 If (.NOT.PSO%iscore(io)) then 3420 IF (.NOT.diracrelativistic) THEN 3421 WRITE(STD_OUT,'(i2,1x,i2,4x,1p,2e15.7)') PSO%np(io),& 3422& PSO%l(io),PSO%occ(io) 3423 ELSE 3424 WRITE(STD_OUT,'(i2,1x,i2,2x,i2,4x,1p,2e15.7)') PSO%np(io),& 3425& PSO%l(io),PSO%kappa(io),PSO%occ(io) 3426 END IF 3427 Endif 3428 Enddo 3429 3430 success=.FALSE. 3431 do loop=1,mxloop 3432 if (TRIM(PSO%exctype)=='HF') then 3433 call PAWIter_HF(Grid,PAW,mix*0.01,err0,err,success) 3434 else 3435 call PAWIter_LDA(Grid,PAW,mix,err0,err,success) 3436 endif 3437 write(std_out,*) '--Results for Iter -- ', loop 3438 IF (.NOT.diracrelativistic) THEN 3439 write(std_out,*) ' n l occupancy energy ' 3440 do io=1,PSO%norbit 3441 If (.NOT.PSO%iscore(io)) then 3442 WRITE(STD_OUT,'(i2,1x,i2,4x,1p,2e15.7)') PSO%np(io),& 3443& PSO%l(io),PSO%occ(io),PSO%eig(io) 3444 endif 3445 enddo 3446 ELSE 3447 write(std_out,*) ' n l kap occupancy energy ' 3448 do io=1,PSO%norbit 3449 If (.NOT.PSO%iscore(io)) then 3450 WRITE(STD_OUT,'(i2,1x,i2,2x,i2,4x,1p,2e15.7)') PSO%np(io),& 3451& PSO%l(io),PSO%kappa(io),PSO%occ(io),PSO%eig(io) 3452 endif 3453 enddo 3454 END IF 3455 IF (success) then 3456 WRITE(STD_OUT,*) ' PS wfn iteration converged ', loop 3457 write(std_out,*) '--Results for Iter -- ', loop 3458 IF (.NOT.diracrelativistic) THEN 3459 write(std_out,*) ' n l occupancy energy ' 3460 do io=1,PSO%norbit 3461 If (.NOT.PSO%iscore(io)) then 3462 WRITE(STD_OUT,'(i2,1x,i2,4x,1p,2e15.7)') PSO%np(io),& 3463& PSO%l(io),PSO%occ(io),PSO%eig(io) 3464 endif 3465 enddo 3466 ELSE 3467 write(std_out,*) ' n l kap occupancy energy ' 3468 do io=1,PSO%norbit 3469 If (.NOT.PSO%iscore(io)) then 3470 WRITE(STD_OUT,'(i2,1x,i2,2x,i2,4x,1p,2e15.7)') PSO%np(io),& 3471& PSO%l(io),PSO%kappa(io),PSO%occ(io),PSO%eig(io) 3472 endif 3473 enddo 3474 ENDIF 3475 exit 3476 ENDIF 3477 enddo 3478 3479 Call mkname(firsttime,stuff) 3480 allocate(tmap(PSO%norbit)) 3481 ip=0; tmap=0 3482 do io=1,PSO%norbit 3483 if (.not.PSO%iscore(io)) then 3484 ip=ip+1 3485 tmap(ip)=io 3486 call PStoAE(Grid,PAW,Grid%n,PSO%l(io),PSO%wfn(:,io),& 3487& PAW%OCCWFN%wfn(:,io)) 3488 endif 3489 enddo 3490 OPEN(unit=1001,file='PAWwfn.'//TRIM(stuff),form='formatted') 3491 do i=1,Grid%n 3492 write(1001,'(1p,60e15.7)') Grid%r(i),(PSO%wfn(i,tmap(k)),& 3493& PAW%OCCWFN%wfn(i,tmap(k)),k=1,ip) 3494 enddo 3495 CLOSE(1001) 3496 3497 deallocate(tmap) 3498 3499 END SUBROUTINE SCFPAW 3500 3501!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3502! PAWIter_LDA(Grid,PAW,err0,err,success) 3503! On input PAW%TOCCWFN contains initial guess of smooth wfn's 3504! On output the guess is updated 3505!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3506 SUBROUTINE PAWIter_LDA(Grid,PAW,mix,err0,err,success) 3507 Type(GridInfo), INTENT(IN) :: Grid 3508 Type(PseudoInfo), INTENT(INOUT) :: PAW 3509 REAL(8) , INTENT(IN) :: mix,err0 3510 REAL(8) , INTENT(OUT) :: err 3511 LOGICAL , INTENT(OUT) :: success 3512 3513 INTEGER :: i,j,k,l,n,io,jo,ib,jb,kb,lb,irc,nbase,nocc 3514 INTEGER, allocatable :: tmap(:) 3515 REAL(8) , ALLOCATABLE :: arg(:),rhs(:),rv(:),aden(:),v1(:),v2(:),o(:,:) 3516 Type(OrbitInfo), POINTER :: PSO 3517 Type(OrbitInfo) :: tmpOrbit 3518 REAL(8) :: occ,x,y,q,v0term,en,val,mix1 3519 INTEGER :: fcount=0 3520 3521 success=.false. 3522 n=Grid%n; nbase=PAW%nbase; irc=PAW%irc 3523 ALLOCATE(arg(n),rhs(n),rv(n),aden(n),v1(n),v2(n),& 3524& tmap(PAW%OCCWFN%norbit),o(PAW%OCCWFN%norbit,nbase)) 3525 3526 PSO=>PAW%TOCCWFN 3527 ! orthonormalize 3528 nocc=0 3529 do io=1,PSO%norbit 3530 if (.not.PSO%iscore(io).and.io>1) then 3531 do jo=1,io-1 3532 if (PSO%l(jo)==PSO%l(io).and..not.PSO%iscore(jo)) then 3533 call genOrthog(Grid,PAW,n,PSO%l(io),& 3534& PSO%wfn(:,io),PSO%wfn(:,jo)) 3535 write(std_out,*) 'orthog', io,jo 3536 endif 3537 enddo 3538 x=genoverlap(Grid,PAW,n,PSO%l(io),PSO%wfn(:,io),PSO%wfn(:,io)) 3539 PSO%wfn(:,io)=PSO%wfn(:,io)/sqrt(x) 3540 CALL ADJUSTSIGN(PSO%wfn(:,io),3) 3541 write(std_out,*) 'normalize ', io, x 3542 nocc=nocc+1 3543 tmap(nocc)=io 3544 endif 3545 enddo 3546 3547 ! prepare overlap terms 3548 o=0 3549 do k=1,nocc 3550 io=tmap(k); l=PSO%l(io) 3551 do ib=1,PAW%nbase 3552 if (l==PAW%l(ib)) then 3553 o(io,ib)=overlap(Grid,PSO%wfn(:,io),PAW%otp(:,ib),1,irc) 3554 write(std_out,'("<p|psi> ", 2i5,1p,e15.7)') io,ib,o(io,ib) 3555 endif 3556 enddo 3557 enddo 3558 3559 Call CopyOrbit(PSO,tmpOrbit) 3560 3561 rv=PAW%rtVf; PAW%tkin=0; PAW%wij=0 3562 PAW%tden=0;rhs=0;aden=0;x=0 3563 do k=1,nocc 3564 io=tmap(k) ; l=PSO%l(io) 3565 occ=PSO%occ(io) 3566 PAW%tden=PAW%tden+occ*(PSO%wfn(:,io))**2 3567 call kinetic_ij(Grid,PSO%wfn(:,io),PSO%wfn(:,io),l,y) 3568 PAW%tkin= PAW%tkin + occ*y 3569 PSO%eig(io)=y 3570 write(std_out,*) 'Kinetic ',io,y 3571 do ib=1,PAW%nbase 3572 do jb=1,PAW%nbase 3573 if (PAW%l(ib)==l.and.PAW%l(jb)==l) then 3574 x=x+occ*o(io,ib)*o(io,jb)*PAW%mLij(ib,jb,1) 3575 PAW%wij(ib,jb)=PAW%wij(ib,jb)+occ*o(io,ib)*o(io,jb) 3576 endif 3577 enddo 3578 enddo 3579 enddo 3580 3581 aden=PAW%tden+x*PAW%g(:,1) 3582 arg=0; arg(2:n)=(aden(2:n)*PAW%hatpot(2:n))/Grid%r(2:n) 3583 v0term=integrator(Grid,arg) 3584 write(std_out,*) 'v0term, aden', v0term,integrator(Grid,aden) 3585 3586 call poisson(Grid,q,aden,rhs,x) 3587 write(std_out,*) 'PAW poisson ', q,x 3588 PAW%tvale=x 3589 arg=0; arg(2:n)=PAW%rtVf(2:n)/Grid%r(2:n) 3590 PAW%tion=overlap(Grid,arg,PAW%tden) 3591 rv=rv+rhs ! vion + valence-Hartree 3592 arg=PAW%tden+PAW%tcore 3593 call exch(Grid,arg,rhs,x,y) 3594 PAW%txc=y 3595 rv=rv+rhs ! + vxc 3596 3597 do k=1,nocc 3598 io=tmap(k);l=PSO%l(io) 3599 arg=rv*(PSO%wfn(:,io)**2); 3600 arg(1)=0; arg(2:n)=arg(2:n)/Grid%r(2:n) 3601 x=integrator(Grid,arg) 3602 write(std_out,*) ' potential term ',io, x 3603 PSO%eig(io)=PSO%eig(io)+x 3604 enddo 3605 3606 PAW%dij=0; PAW%Ea=0 3607 do ib=1,PAW%nbase 3608 do jb=1,PAW%nbase 3609 if (PAW%l(ib)==PAW%l(jb)) then 3610 PAW%dij(ib,jb)=PAW%dij(ib,jb) + PAW%Kij(ib,jb) + & 3611& PAW%VFij(ib,jb)+PAW%mLij(ib,jb,1)*v0term 3612 PAW%Ea=PAW%Ea + PAW%wij(ib,jb)*(PAW%Kij(ib,jb) + & 3613& PAW%VFij(ib,jb)) 3614 x=0; ! accumulate Hartree term 3615 do kb=1,PAW%nbase 3616 do lb=1,PAW%nbase 3617 if (PAW%l(kb)==PAW%l(lb)) then 3618 x=x+PAW%wij(kb,lb)*PAW%DR(ib,jb,kb,lb,1) 3619 endif 3620 enddo 3621 enddo 3622 PAW%dij(ib,jb)=PAW%dij(ib,jb)+x 3623 PAW%Ea=PAW%Ea + 0.5d0*PAW%wij(ib,jb)*x 3624 endif 3625 enddo 3626 enddo 3627 3628 ! exchange-correlation part 3629 arg=PAW%core; rhs=PAW%tcore 3630 do ib=1,PAW%nbase 3631 do jb=1,PAW%nbase 3632 if (PAW%l(ib)==PAW%l(jb)) then 3633 arg=arg+PAW%wij(ib,jb)*PAW%ophi(:,ib)*PAW%ophi(:,jb) 3634 rhs=rhs+PAW%wij(ib,jb)*PAW%otphi(:,ib)*PAW%otphi(:,jb) 3635 endif 3636 enddo 3637 enddo 3638 3639 Write(STD_OUT,*) 'Before EXC ', PAW%Ea 3640 irc=PAW%irc 3641 call exch(Grid,arg,v1,x,y,fin=irc) 3642 PAW%Ea=PAW%Ea+y ; write(std_out,*) 'AE EXC ' ,y 3643 call exch(Grid,rhs,v2,x,y,fin=irc) 3644 PAW%Ea=PAW%Ea-y ; write(std_out,*) 'PS EXC ' ,y 3645 3646 PAW%Etotal=PAW%tkin+PAW%tion+PAW%tvale+PAW%txc+PAW%Ea 3647 Write(STD_OUT,*) '*******Total energy*********', PAW%Etotal 3648 Write(STD_OUT,*) 'PAW%tkin ',PAW%tkin 3649 Write(STD_OUT,*) 'PAW%tion ',PAW%tion 3650 Write(STD_OUT,*) 'PAW%tvale ',PAW%tvale 3651 Write(STD_OUT,*) 'PAW%txc ',PAW%txc 3652 Write(STD_OUT,*) 'PAW%Ea ',PAW%Ea 3653 3654 do ib=1,PAW%nbase 3655 do jb=1,PAW%nbase 3656 if (PAW%l(ib)==PAW%l(jb)) then 3657 arg=PAW%ophi(:,ib)*PAW%ophi(:,jb)*v1(:) & 3658& - PAW%otphi(:,ib)*PAW%otphi(:,jb)*v2(:) 3659 arg(2:n)=arg(2:n)/Grid%r(2:n) 3660 call extrapolate(Grid,arg) 3661 PAW%dij(ib,jb)=PAW%dij(ib,jb)+integrator(Grid,arg,1,irc) 3662 endif 3663 enddo 3664 enddo 3665 3666 3667 ! complete estimate of PSO%eig(io) 3668 do k=1,nocc 3669 io=tmap(k); l=PSO%l(io) 3670 write(std_out,*) 'Eig before dij ', PSO%eig(io) 3671 do ib=1,PAW%nbase 3672 do jb=1,PAW%nbase 3673 if (PAW%l(ib)==l.and.PAW%l(jb)==l) then 3674 PSO%eig(io)=PSO%eig(io)+PAW%dij(ib,jb)*o(io,ib)*o(io,jb) 3675 endif 3676 enddo 3677 enddo 3678 write(std_out,*) 'Eig after dij ', PSO%eig(io) 3679 enddo 3680 3681 ! Solve inhomogeneous diffeq. and store result in tmpOrbit 3682 err=0; 3683 do k=1,nocc 3684 io=tmap(k); l=PSO%l(io) 3685 en=PSO%eig(io) 3686 rhs=0 3687 do ib=1,PAW%nbase 3688 do jb=1,PAW%nbase 3689 if (PAW%l(ib)==l.and.PAW%l(jb)==l) then 3690 rhs=rhs-PAW%otp(:,ib)*(en*PAW%oij(ib,jb)-& 3691& PAW%dij(ib,jb))*o(io,jb) 3692 endif 3693 enddo 3694 enddo 3695 tmpOrbit%wfn(:,io)=0 3696 ! note rhs is -(what we expect) 3697 CALL inhomo_bound_numerov(Grid,l,n,en,rv,rhs,tmpOrbit%wfn(:,io)) 3698 arg=(PSO%wfn(:,io)-tmpOrbit%wfn(:,io))**2 3699 err=err+PSO%occ(io)*Integrator(Grid,arg) 3700 !do i=1,Grid%n 3701 ! write(800+k,'(1p,20e15.7)') Grid%r(i),PSO%wfn(i,io),& 3702 !& tmpOrbit%wfn(i,io),rhs(i),rv(i) 3703 !enddo 3704 enddo 3705 3706 write(std_out,*) 'PAWIter ', fcount,err 3707 ! update wfn if tolerance not satisfied 3708 IF (err>err0) THEN 3709 mix1=MAX(mix,mix/err) 3710 val=(1.d0-mix1) 3711 WRITE(STD_OUT,*) 'mixing wfns ', val 3712 DO k=1,nocc 3713 io=tmap(k) 3714 PSO%wfn(:,io)=val*PSO%wfn(:,io)+mix1*tmpOrbit%wfn(:,io) 3715 ENDDO 3716 ELSE 3717 success=.TRUE. 3718 ENDIF 3719 fcount=fcount+1 3720 DEALLOCATE(arg,rhs,rv,aden,v1,v2,tmap,o) 3721 Call DestroyOrbit(tmpOrbit) 3722 END SUBROUTINE PAWIter_LDA 3723 3724 3725!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3726! exploreparms reads ndata sets of PAW parameters and calculates 3727! the corresponding logderivs It is intended that this routine 3728! would be used to choose the parameters with the best logderivs. 3729! It should be called after the program has first run a "normal" 3730! calculation. exploreparms recalculates the wfni (i=1,2,..) and 3731! logderivl (l=0,1,2,..) values but does not generate the full 3732! PAW dataset. Upon analyzing the logderiv results, the atompaw 3733! program should be run in the "normal" mode to generate the best 3734! dataset. exploreparams should only be called once. 3735! 3736!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3737 3738 3739 SUBROUTINE exploreparms(Grid,Pot,FC,Orbit,PAW) 3740 Type(GridInfo), INTENT(IN) :: Grid 3741 Type(PotentialInfo), INTENT(IN) :: Pot 3742 Type(FCInfo), INTENT(IN) :: FC 3743 Type (OrbitInfo), INTENT(IN) :: Orbit 3744 Type(PseudoInfo), INTENT(INOUT) :: PAW 3745 3746 INTEGER :: ndata,i,j,k,l,nrcs,ircs,ictrc,iprev 3747 INTEGER :: ifen=9 3748 CHARACTER(4) ::fdata 3749 CHARACTER(132) :: inputline,keyword 3750 REAL(8), allocatable :: logderiverror(:,:) 3751 REAL(8) :: thisrc 3752 REAL(8) :: EBEGIN=-10, EEND=10 ! default logderiv range (Ry) 3753 LOGICAL :: success 3754 TYPE(input_dataset_t) :: backup_dataset 3755 3756 Type rcresults 3757 INTEGER :: beginindex, endindex 3758 REAL(8) :: rc 3759 End type rcresults 3760 3761 Type (rcresults), POINTER :: rcsummary(:) 3762 3763 success=.true. 3764 write(std_out,*) 'Input the number of PAW parameter sets in this run' 3765 READ(STD_IN,'(a)') inputline 3766 CALL eliminate_comment(inputline) 3767 read(inputline,*) ndata,nrcs 3768 if (ndata>9999) then 3769 write(std_out,*) 'Error : ndata must be <= 9999', ndata 3770 stop 3771 endif 3772 if (nrcs<1.or.nrcs>9999) then 3773 write(std_out,*) 'Error : nrcs must be <= 9999 and >=1', nrcs 3774 stop 3775 endif 3776 3777 call UpperCase(inputline) 3778 i=0; i=INDEX(inputline,'LOGDERIVERANGE') 3779 if (i>0) then 3780 read(unit=inputline(i+14:80),fmt=*,err=111,end=111,iostat=i) & 3781& EBEGIN , EEND 3782 111 continue 3783 endif 3784 3785 allocate(rcsummary(nrcs)) 3786 allocate(logderiverror(6,ndata)) ! 6 represents the max l+1 3787 logderiverror=9.d20 3788 3789 write(std_out,*) 'Begin explore runs' 3790 OPEN(20,file='EXPLORERESULTS',form='formatted') 3791 OPEN(21,file='EXPLORESUMMARY',form='formatted') 3792 write(20,'("#dataset"," rc ",6(2x,i3,10x))') (l,l=0,PAW%lmax+1) 3793 write(21,'(" Logderiv errors based on energy range", 2f12.2)') & 3794& EBEGIN, EEND 3795 3796! Save input dataset 3797 CALL input_dataset_copy(input_dataset,backup_dataset) 3798 3799 ircs=0; thisrc=-1 3800 do i=1,ndata 3801 write(std_out,*) '===================== #',i,'==================' 3802 call mkname(i,fdata) 3803 OPEN(ifen, file='EXPLOREout.'//TRIM(fdata), form='formatted') 3804 3805 3806 !Read new basis parameters 3807 CALL input_dataset_read(echofile='EXPLOREIN.'//TRIM(fdata),& 3808& read_global_data=.false.,read_elec_data=.false.,& 3809& read_coreval_data=.false.,read_basis_data=.true.) 3810 3811 CALL SetPAWOptions1(ifen,Grid) 3812 Call InitPAW(PAW,Grid,FCOrbit) 3813 CALL setbasis(Grid,FCPot,FCOrbit) 3814 Call setcoretail(Grid,FC%coreden) 3815 Call setttau(Grid,FC%coretau) 3816 If (TRIM(FCorbit%exctype)=='HF'.or.TRIM(FCorbit%exctype)=='EXXKLI') PAW%tcore=0 3817 If (TRIM(FCorbit%exctype)=='EXXKLI') Call fixtcorewfn(Grid,PAW) 3818 Call SetPAWOptions2(ifen,Grid,FCOrbit,FCPot,success) 3819 Call Report_PseudobasisRP(Grid,PAW,ifen,fdata) 3820 if (success) then 3821 Call Set_PAW_MatrixElements(Grid,PAW,ifen) 3822 CALL EXPLORElogderiv(Grid,FCPot,PAW,fdata,EBEGIN,EEND,& 3823& logderiverror(:,i)) 3824 endif 3825 write(20,'(i5,2x,f12.5,1p,6e15.7)')& 3826& i,PAW%rc,(logderiverror(l+1,i),l=0,PAW%lmax+1) 3827 close(ifen) 3828 Call DestroyPAW(PAW) 3829 if (abs(PAW%rc-thisrc)>1.d-10 ) then 3830 ircs=ircs+1 3831 rcsummary(ircs)%beginindex=i 3832 rcsummary(ircs)%endindex=i 3833 thisrc=PAW%rc 3834 rcsummary(ircs)%rc=thisrc 3835 else 3836 rcsummary(ircs)%rc=thisrc 3837 rcsummary(ircs)%endindex=i 3838 endif 3839 write(std_out,*) 'test ', i,ircs, rcsummary(ircs)%beginindex,& 3840& rcsummary(ircs)%endindex 3841 3842 enddo 3843 3844! Restore input dataset 3845 CALL input_dataset_copy(backup_dataset,input_dataset) 3846 3847 Write(STD_OUT,*) 'Results for minimum logderiverror' 3848 Write(21,*) 'Results for minimum logderiverror' 3849 Do ircs=1,nrcs 3850 write(std_out,'( "=== Rc = ", f20.6, " ====")') rcsummary(ircs)%rc 3851 write(21,'( "=== Rc = ", f20.6, " ====")') rcsummary(ircs)%rc 3852 j=rcsummary(ircs)%beginindex 3853 k=rcsummary(ircs)%endindex 3854 Do l=0,PAW%lmax+1 3855 i=MINLOC(logderiverror(l+1,j:k),1)+j-1 3856 Write(STD_OUT,'(" l =", i5,2x, i6, 1pe15.7)') l,i,logderiverror(l+1,i) 3857 Write(21,'(" l =", i5,2x, i6, 1pe15.7)') l,i,logderiverror(l+1,i) 3858 enddo 3859 enddo 3860 3861 CLOSE(20); CLOSE(21) 3862 Deallocate(logderiverror,rcsummary) 3863 3864 END SUBROUTINE exploreparms 3865 3866 3867 !************************************************************************ 3868 ! program to calculate logerivatives of paw wavefunctions 3869 ! and to compare them with all electron wavefunctions 3870 ! optionally, wavefunctions are written to a file 3871 ! Assumes prior call to SUBROUTINE Set_PAW_MatrixElements(Grid,PAW) 3872 ! Version for exploreparms 3873 !************************************************************************ 3874 SUBROUTINE EXPLORElogderiv(Grid,Pot,PAW,label,EBEGIN,EEND,lderror) 3875 TYPE(GridInfo), INTENT(IN) :: Grid 3876 TYPE(PotentialInfo), INTENT(IN) :: Pot 3877 TYPE(PseudoInfo), INTENT(INOUT) :: PAW 3878 CHARACTER(4), INTENT(IN) :: label 3879 REAL(8), INTENT(IN) :: EBEGIN,EEND 3880 REAL(8), INTENT(INOUT) :: lderror(:) 3881 3882 TYPE(PotentialInfo) :: PS 3883 INTEGER :: n,l,ie,nbase,ib,ic,nr,i,nodes,mbase,irc,lng,ne 3884 REAL(8), PARAMETER :: e0=-5.d0,de=0.05d0 3885 REAL(8) :: h,x,dwdr,dcwdr,scale,energy 3886 REAL(8), ALLOCATABLE :: psi(:),tpsi(:),ttpsi(:) 3887 CHARACTER(4) :: flnm 3888 LOGICAL :: OK 3889 3890 n=Grid%n; h=Grid%h; nbase=PAW%nbase; irc=PAW%irc; nr=irc+10 3891 ne=(EEND-EBEGIN)/de + 1 3892 3893 ALLOCATE(psi(nr),tpsi(nr),ttpsi(nr),PS%rv(n),stat=ie) 3894 IF (ie/=0) THEN 3895 WRITE(STD_OUT,*) 'Error in logderiv allocation',n,ie 3896 STOP 3897 ENDIF 3898 3899 ! load PS 3900 PS%rv=PAW%rveff ; PS%nz=0.d0 3901 call zeropot(Grid,PS%rv,PS%v0,PS%v0p) 3902 ! 3903 ! calculate logderivatives at irc 3904 ! 3905 WRITE(STD_OUT,*) 'calculating log derivatives at irc',Grid%r(irc) 3906 ! 3907 mbase=nbase; irc=PAW%irc 3908 write(std_out,*) 'Nodes counted to radius ', Grid%r(irc) 3909 DO l=0,PAW%lmax+1 3910 OK=.true. 3911 If (l<=PAW%lmax) then 3912 Basislist: do ib=1,nbase 3913 if (PAW%l(ib)==l) then 3914 nodes=countnodes(2,irc,PAW%ophi(:,ib)) 3915 if (nodes/=PAW%nodes(ib).and.PAW%eig(ib)>0) then 3916 lderror(l+1)=9.d20 3917 write(std_out,*) 'Warning node problems for case ', l,ib,& 3918& nodes,PAW%nodes(ib) 3919 OK=.false. 3920 exit Basislist 3921 endif 3922 endif 3923 enddo Basislist 3924 endif 3925 if (OK) then 3926 CALL mkname(l,flnm) 3927 OPEN(56,file=TRIM(label)//'.logderiv.'//TRIM(flnm),form='formatted') 3928 lderror(l+1)=0 3929 3930 DO ie=1,ne 3931 energy=EBEGIN+de*(ie-1) 3932 psi=0;tpsi=0;ttpsi=0 3933 if (scalarrelativistic) then 3934 CALL unboundsr(Grid,Pot,nr,l,energy,psi,nodes) 3935 elseif (TRIM(PAW%exctype)=='HF') then 3936 CALL HFunocc(Grid,PAW%OCCWFN,l,energy,Pot%rv,Pot%v0,Pot%v0p,& 3937 psi,lng) 3938 else 3939 CALL unboundsch(Grid,Pot%rv,Pot%v0,Pot%v0p,nr,l,energy,psi,nodes) 3940 endif 3941 CALL unboundsep(Grid,PS,PAW,nr,l,energy,tpsi,nodes) 3942 CALL PStoAE(Grid,PAW,nr,l,tpsi,ttpsi) 3943 ! 3944 dwdr=Gfirstderiv(Grid,irc,psi)/psi(irc) 3945 dcwdr=Gfirstderiv(Grid,irc,ttpsi)/ttpsi(irc) 3946 3947 WRITE(56,'(1p,5e12.4)') energy,dwdr,dcwdr ; 3948 lderror(l+1)=lderror(l+1)+abs(atan(dwdr)-atan(dcwdr)) 3949 ENDDO !ie 3950 CLOSE(56) 3951 ENDIF !OK 3952 3953 ENDDO !l 3954 3955 DEALLOCATE(psi,tpsi,ttpsi,PS%rv) 3956 END SUBROUTINE EXPLORElogderiv 3957 3958 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3959 ! Repeated version of output for used with EXPLORElogderiv 3960 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3961 3962 SUBROUTINE Report_PseudobasisRP(Grid,PAW,ifen,fdata) 3963 Type(GridInfo), INTENT(IN) :: Grid 3964 Type(PseudoInfo), INTENT(IN) :: PAW 3965 INTEGER, INTENT(IN) :: ifen 3966 CHARACTER(4) :: fdata 3967 3968 INTEGER, parameter :: ifout=15 3969 INTEGER :: i,j,io,nbase,irc,icount,n 3970 INTEGER, ALLOCATABLE :: mapp(:) 3971 CHARACTER (len=4) :: flnm 3972 3973 nbase=PAW%nbase;irc=PAW%irc;n=Grid%n 3974 WRITE(ifen,'(/"Number of basis functions ",i5)') nbase 3975 WRITE(ifen,*)'No. n l Energy Cp coeff Occ' 3976 3977 DO io=1,nbase 3978 WRITE(ifen,'(3i5,1p,3e15.7)') io,PAW%np(io),PAW%l(io),PAW%eig(io),& 3979& PAW%ck(io),PAW%occ(io) 3980 CALL mkname(io,flnm) 3981 OPEN(ifout,file=TRIM(fdata)//'.wfn'//TRIM(flnm),form='formatted') 3982 WRITE(ifout,*) '# l=',PAW%l(io),'basis function with energy ',& 3983& PAW%eig(io) 3984 DO i=1,irc+50 3985 WRITE(ifout,'(1p,5e12.4)') Grid%r(i),PAW%ophi(i,io),& 3986& PAW%otphi(i,io),PAW%otp(i,io) 3987 ENDDO 3988 CLOSE(ifout) 3989 ENDDO 3990 3991 ! also write "raw" wavefunctions 3992 DO io=1,nbase 3993 CALL mkname(io,flnm) 3994 OPEN(ifout,file=TRIM(fdata)//'.wfn00'//TRIM(flnm),form='formatted') 3995 WRITE(ifout,*) '# l=',PAW%l(io),'basis function with energy ',& 3996& PAW%eig(io) 3997 DO i=1,irc+50 3998 WRITE(ifout,'(1p,5e12.4)') Grid%r(i),PAW%phi(i,io),& 3999& PAW%tphi(i,io),PAW%tp(i,io) 4000 ENDDO 4001 CLOSE(ifout) 4002 ENDDO 4003 4004 allocate(mapp(PAW%OCCWFN%norbit)) 4005 mapp=0 4006 icount=0 4007 do io=1,PAW%OCCWFN%norbit 4008 if(PAW%OCCWFN%iscore(io)) then 4009 else 4010 icount=icount+1 4011 mapp(icount)=io 4012 endif 4013 enddo 4014 OPEN(ifout,file='OCCWFN',form='formatted') 4015 WRITE(ifout,'("# ",50i30)') (mapp(j),j=1,icount) 4016 do i=1,n 4017 write(ifout,'(1p,51e15.7)') Grid%r(i),(PAW%OCCWFN%wfn(i,mapp(j)),& 4018& PAW%TOCCWFN%wfn(i,mapp(j)),j=1,icount) 4019 enddo 4020 close(ifout) 4021 4022 deallocate(mapp) 4023 4024 END SUBROUTINE Report_PseudobasisRP 4025 4026 4027!************************************************************************ 4028! 7/2018 Program written by Casey Brock from Vanderbilt U. 4029! Unwraps phase so it is continuous (no jumps of pi) 4030! The algorithm was designed by Alan Tackett. 4031! 4032! x, y: psi and psi' used to calculate atan 4033! atan_curr, atan_prev: value of atan, and values from prev. iteration 4034! These should always be unshifted. 4035! s2, s1: current (s2) and previous (s1) values of atan slope 4036! cumshift: an integer multiple of pi, the cumulative shift applied 4037! to atan 4038! cnt: the number of previous points that were shifted 4039!************************************************************************ 4040SUBROUTINE phase_unwrap(ie, y, x, atan_curr, atan_prev, s2, s1, cumshift, cnt) 4041 INTEGER, INTENT(IN) :: ie 4042 REAL(8), INTENT(IN) :: x, y, s1, s2 4043 REAL(8), INTENT(IN) :: atan_curr, atan_prev 4044 REAL(8), INTENT(INOUT) :: cumshift 4045 INTEGER, INTENT(INOUT) :: cnt 4046 4047 cnt = cnt + 1 4048 4049! if (x<0) AND (slope changes sign) 4050 IF ((x<0) .AND. (atan_curr*atan_prev < 0)) THEN 4051 cnt = 0 4052 cumshift = cumshift - sign(pi, y) 4053! if (slope changes sign) AND (last two points weren't shifted) 4054! AND (absolute change in slope > 0.01) 4055 ELSE IF ((s1*s2 < 0) .AND. (cnt>2) .AND. (ABS(s2-s1)>0.01)) THEN 4056 cnt = 0 4057 cumshift = cumshift + sign(pi, s1) 4058 ENDIF 4059END SUBROUTINE phase_unwrap 4060 4061END MODULE pseudo 4062