1!! ----------------- Show the list of all supported real space functions 2subroutine funclist 3use defvar 4write(*,*) " ----------- Avaliable real space functions -----------" 5if (allocated(b)) then 6 write(*,*) "1 Electron density 2 Gradient norm of electron density" 7 write(*,*) "3 Laplacian of electron density 4 Value of orbital wavefunction" 8 if (ipolarpara==0) write(*,*) "5 Electron spin density" 9 if (ipolarpara==1) write(*,*) "5 Spin polarization parameter function" 10 write(*,*) "6 Hamiltonian kinetic energy density K(r)" 11 write(*,*) "7 Lagrangian kinetic energy density G(r)" 12 if (ifiletype==4) then 13 write(*,*) "8 Electrostatic potential from atomic charges" 14 else 15 write(*,*) "8 Electrostatic potential from nuclear charges" 16 end if 17 if (ELFLOL_type==0) write(*,*) "9 Electron Localization Function (ELF)" 18 if (ELFLOL_type==1) write(*,*) "9 Electron Localization Function (ELF) defined by Tsirelson" 19 if (ELFLOL_type==2) write(*,*) "9 Electron Localization Function (ELF) defined by Lu, Tian" 20 if (ELFLOL_type==0) write(*,*) "10 Localized orbital locator (LOL)" 21 if (ELFLOL_type==1) write(*,*) "10 Localized orbital locator (LOL) defined by Tsirelson" 22 if (ELFLOL_type==2) write(*,*) "10 Localized orbital locator (LOL) defined by Lu, Tian" 23 write(*,*) "11 Local information entropy" 24 write(*,*) "12 Total electrostatic potential (ESP)" 25 write(*,*) "13 Reduced density gradient (RDG) 14 RDG with promolecular approximation" 26 write(*,*) "15 Sign(lambda2)*rho 16 Sign(lambda2)*rho with promolecular approximation" 27 !Fermi hole function only available to single-determinant wavefunction 28 if (pairfunctype==1) write(*,"(a,3f10.5)") " 17 Correlation hole for alpha, ref. point:",refx,refy,refz 29 if (pairfunctype==2) write(*,"(a,3f10.5)") " 17 Correlation hole for beta, ref. point:",refx,refy,refz 30 if (pairfunctype==4) write(*,"(a,3f10.5)") " 17 Correlation factor for alpha, ref. point:",refx,refy,refz 31 if (pairfunctype==5) write(*,"(a,3f10.5)") " 17 Correlation factor for beta, ref. point:",refx,refy,refz 32 if (pairfunctype==7) write(*,"(a,3f10.5)") " 17 Exc.-corr. density for alpha, ref. point:",refx,refy,refz 33 if (pairfunctype==8) write(*,"(a,3f10.5)") " 17 Exc.-corr. density for beta, ref. point:",refx,refy,refz 34 if (pairfunctype==10) write(*,"(a,3f10.5)") " 17 Pair density for alpha, ref. point:",refx,refy,refz 35 if (pairfunctype==11) write(*,"(a,3f10.5)") " 17 Pair density for beta, ref. point:",refx,refy,refz 36 if (pairfunctype==12) write(*,"(a,3f10.5)") " 17 Pair density for all electrons, ref. point:",refx,refy,refz 37 write(*,*) "18 Average local ionization energy" 38 write(*,"(a,i2,a,3f10.5)") " 19 Source function, mode:",srcfuncmode,", ref. point:",refx,refy,refz 39 write(*,*) "20 Electron delocalization range function EDR(r;d)" 40 write(*,*) "21 Orbital overlap distance function D(r)" 41 write(*,"(a,i5)") " 100 User-defined real space function, iuserfunc=",iuserfunc 42else !No wavefunction information is available 43 if (ifiletype==4) then 44 write(*,*) "8 Electrostatic potential from atomic charges" 45 else 46 write(*,*) "8 Electrostatic potential from nuclear charges" 47 end if 48 write(*,*) "14 Reduced density gradient(RDG) with promolecular approximation" 49 write(*,*) "16 Sign(lambda2)*rho with promolecular approximation" 50 write(*,"(a,i3)") " 100 User-defined real space function, iuserfunc=",iuserfunc 51end if 52end subroutine 53 54 55!---- Standard interface for selecting real space function 56!Note that iorbsel is a global variable 57subroutine selfunc_interface(ifunc) 58use defvar 59integer ifunc,edrmaxpara,wrtnumedr 60real*8 wrtstart 61call funclist 62read(*,*) ifunc 63if (ifunc==4) then 64 write(*,"(a,i10)") " Input orbital index, between 1 and",nmo 65 read(*,*) iorbsel 66else if (ifunc==20) then !Read length scale to evaluate EDR(r;d) 67 write(*,*) "The EDR(r;d) computing code was contributed by Arshad Mehmood" 68 write(*,"(a,/)") " References: J. Chem. Phys., 141, 144104 (2014); J. Chem. Theory Comput., 12, 79 (2016); Angew. Chem. Int. Ed., 56, 6878 (2017)" 69 write(*,*) "Input length scale d (Bohr) e.g. 0.85" 70 read(*,*) dedr 71else if (ifunc==21) then 72 write(*,*) "The D(r) computing code was contributed by Arshad Mehmood" 73 write(*,"(a,/)") " References: J. Chem. Theory Comput., 12, 3185 (2016); Phys. Chem. Chem. Phys., 17, 18305 (2015)" 74 write(*,*) "1 Manually input total number, start and increment in EDR exponents" 75 write(*,*) "2 Use default values i.e. 20,2.50,1.50" 76 read(*,*) edrmaxpara 77 if (edrmaxpara==1) then 78 write(*,*) "Please input in order: exponents start increment e.g. 20 2.5 1.5" 79 write(*,*) "Note: Max. allowed exponents are 50 and min. allowed increment is 1.01" 80 read(*,*) nedr,edrastart,edrainc 81 if (nedr<1) then 82 write(*,*) "Error: Bad Number of EDR exponents. Should be between 1 to 50" 83 write(*,*) "Press ENTER to exit" 84 read(*,*) 85 stop 86 else if (nedr>50) then 87 write(*,*) "Error: Bad Number of EDR exponents. Should be between 1 to 50" 88 write(*,*) "Press ENTER to exit" 89 read(*,*) 90 stop 91 end if 92 if (edrainc<1.01d0) then 93 write(*,*) "Error: Bad increment in EDR exponents. Should not be less than 1.01" 94 write(*,*) "Press ENTER to exit" 95 read(*,*) 96 stop 97 end if 98 else if (edrmaxpara==2) then 99 nedr=20 100 edrastart=2.5d0 101 edrainc=1.5d0 102 end if 103 write(*,*) "The following EDR exponents will be used in calculation:" 104 wrtstart=edrastart 105 do wrtnumedr=1,nedr 106 wrtexpo(wrtnumedr)=wrtstart 107 wrtstart=wrtstart/edrainc 108 write(*,"(E13.5)") wrtexpo(wrtnumedr) 109 end do 110 write(*,*) 111end if 112end subroutine 113 114 115 116 117 118 119! ====================================== 120! =========== Module function ========== 121! ====================================== 122 123module function 124use defvar 125implicit real*8 (a-h,o-z) 126 127contains 128 129 130!-------- Calculate any supported real space function at a given point 131real*8 function calcfuncall(ifunc,x,y,z) 132integer ifunc 133real*8 x,y,z 134if (ifunc==1) then 135calcfuncall=fdens(x,y,z) 136else if (ifunc==2) then 137calcfuncall=fgrad(x,y,z,'t') 138else if (ifunc==3) then 139calcfuncall=flapl(x,y,z,'t') 140else if (ifunc==4) then 141calcfuncall=fmo(x,y,z,iorbsel) 142else if (ifunc==5) then 143calcfuncall=fspindens(x,y,z,'s') 144else if (ifunc==6) then 145calcfuncall=Hamkin(x,y,z,0) 146else if (ifunc==7) then 147calcfuncall=lagkin(x,y,z,0) 148else if (ifunc==8) then 149calcfuncall=nucesp(x,y,z) 150else if (ifunc==9) then 151calcfuncall=ELF_LOL(x,y,z,"ELF") 152else if (ifunc==10) then 153calcfuncall=ELF_LOL(x,y,z,"LOL") 154else if (ifunc==11) then 155calcfuncall=infoentro(1,x,y,z) 156else if (ifunc==12) then 157calcfuncall=totesp(x,y,z) 158else if (ifunc==13) then 159calcfuncall=fgrad(x,y,z,'r') 160else if (ifunc==14) then 161calcfuncall=RDGprodens(x,y,z) 162else if (ifunc==15) then 163calcfuncall=signlambda2rho(x,y,z) 164else if (ifunc==16) then 165calcfuncall=signlambda2rho_prodens(x,y,z) 166else if (ifunc==17) then 167calcfuncall=pairfunc(refx,refy,refz,x,y,z) 168else if (ifunc==18) then 169calcfuncall=avglocion(x,y,z) 170else if (ifunc==19) then 171calcfuncall=srcfunc(x,y,z,srcfuncmode) 172else if (ifunc==20) then 173calcfuncall=edr(x,y,z) 174else if (ifunc==21) then 175calcfuncall=edrdmax(x,y,z) 176else if (ifunc==100) then 177calcfuncall=userfunc(x,y,z) 178end if 179end function 180 181!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 182!! Calculate wavefunction value of a range of orbitals and their derivatives at a given point, up to third-order 183!! istart and iend is the range of the orbitals will be calculated, to calculate all orbitals, use 1,nmo 184!! runtype=1: value =2: value+dx/y/z =3: value+dxx/yy/zz(diagonal of hess) =4: value+dx/y/z+Hessian 185!! =5: value+dx/y/z+hess+3-order derivative tensor 186subroutine orbderv(runtype,istart,iend,x,y,z,wfnval,grad,hess,tens3) 187real*8 x,y,z,wfnval(nmo) 188real*8,optional :: grad(3,nmo),hess(3,3,nmo),tens3(3,3,3,nmo) 189integer runtype,istart,iend 190 191wfnval=0D0 192if (present(grad)) grad=0D0 193if (present(hess)) hess=0D0 194if (present(tens3)) tens3=0D0 195lastcen=-1 !arbitrary value 196 197! if the center/exp of current GTF is the same as previous, then needn't recalculate them 198do j=1,nprims 199 ix=type2ix(b(j)%functype) 200 iy=type2iy(b(j)%functype) 201 iz=type2iz(b(j)%functype) 202 ep=b(j)%exp 203 204 if (b(j)%center/=lastcen) then 205 sftx=x-a(b(j)%center)%x 206 sfty=y-a(b(j)%center)%y 207 sftz=z-a(b(j)%center)%z 208 sftx2=sftx*sftx 209 sfty2=sfty*sfty 210 sftz2=sftz*sftz 211 rr=sftx2+sfty2+sftz2 212 end if 213 if (expcutoff>0.or.-ep*rr>expcutoff) then 214 expterm=exp(-ep*rr) 215 else 216 expterm=0D0 217 end if 218 lastcen=b(j)%center 219! expterm=exp(-ep*dsqrt(rr)) 220 if (expterm==0D0) cycle 221 222 !Calculate value for current GTF 223 if (b(j)%functype==1) then !Some functype use manually optimized formula for cutting down computational time 224 GTFval=expterm 225 else if (b(j)%functype==2) then 226 GTFval=sftx*expterm 227 else if (b(j)%functype==3) then 228 GTFval=sfty*expterm 229 else if (b(j)%functype==4) then 230 GTFval=sftz*expterm 231 else if (b(j)%functype==5) then 232 GTFval=sftx2*expterm 233 else if (b(j)%functype==6) then 234 GTFval=sfty2*expterm 235 else if (b(j)%functype==7) then 236 GTFval=sftz2*expterm 237 else if (b(j)%functype==8) then 238 GTFval=sftx*sfty*expterm 239 else if (b(j)%functype==9) then 240 GTFval=sftx*sftz*expterm 241 else if (b(j)%functype==10) then 242 GTFval=sfty*sftz*expterm 243 else !If above condition is not satisfied(Angular moment higher than f), the function will calculated explicitly 244 GTFval=sftx**ix *sfty**iy *sftz**iz *expterm 245 end if 246 !Calculate orbital wavefunction value 247 do imo=istart,iend 248 wfnval(imo)=wfnval(imo)+co(imo,j)*GTFval 249 end do 250 251 if (runtype>=2) then 252 !Calculate 1-order derivative for current GTF 253 tx=0.0D0 254 ty=0.0D0 255 tz=0.0D0 256 if (ix/=0) tx=ix*sftx**(ix-1) 257 if (iy/=0) ty=iy*sfty**(iy-1) 258 if (iz/=0) tz=iz*sftz**(iz-1) 259 GTFdx=sfty**iy *sftz**iz *expterm*(tx-2*ep*sftx**(ix+1)) 260 GTFdy=sftx**ix *sftz**iz *expterm*(ty-2*ep*sfty**(iy+1)) 261 GTFdz=sftx**ix *sfty**iy *expterm*(tz-2*ep*sftz**(iz+1)) 262 !Calculate 1-order derivative for orbitals 263 do imo=istart,iend 264 grad(1,imo)=grad(1,imo)+co(imo,j)*GTFdx 265 grad(2,imo)=grad(2,imo)+co(imo,j)*GTFdy 266 grad(3,imo)=grad(3,imo)+co(imo,j)*GTFdz 267 end do 268 269 if (runtype>=3) then 270 !Calculate 2-order derivative for current GTF 271 txx=0.0D0 272 tyy=0.0D0 273 tzz=0.0D0 274 if (ix>=2) txx=ix*(ix-1)*sftx**(ix-2) 275 if (iy>=2) tyy=iy*(iy-1)*sfty**(iy-2) 276 if (iz>=2) tzz=iz*(iz-1)*sftz**(iz-2) 277 GTFdxx=sfty**iy *sftz**iz *expterm*( txx + 2*ep*sftx**ix*(-2*ix+2*ep*sftx2-1) ) 278 GTFdyy=sftx**ix *sftz**iz *expterm*( tyy + 2*ep*sfty**iy*(-2*iy+2*ep*sfty2-1) ) 279 GTFdzz=sftx**ix *sfty**iy *expterm*( tzz + 2*ep*sftz**iz*(-2*iz+2*ep*sftz2-1) ) 280 ttx=tx-2*ep*sftx**(ix+1) 281 tty=ty-2*ep*sfty**(iy+1) 282 ttz=tz-2*ep*sftz**(iz+1) 283 GTFdxy=sftz**iz *expterm*ttx*tty 284 GTFdyz=sftx**ix *expterm*tty*ttz 285 GTFdxz=sfty**iy *expterm*ttx*ttz 286 !Calculate diagonal Hessian elements for orbitals 287 do imo=istart,iend 288 hess(1,1,imo)=hess(1,1,imo)+co(imo,j)*GTFdxx !dxx 289 hess(2,2,imo)=hess(2,2,imo)+co(imo,j)*GTFdyy !dyy 290 hess(3,3,imo)=hess(3,3,imo)+co(imo,j)*GTFdzz !dzz 291 end do 292 if (runtype>=4) then !Also process nondiagonal elements 293 do imo=istart,iend 294 hess(1,2,imo)=hess(1,2,imo)+co(imo,j)*GTFdxy !dxy 295 hess(2,3,imo)=hess(2,3,imo)+co(imo,j)*GTFdyz !dyz 296 hess(1,3,imo)=hess(1,3,imo)+co(imo,j)*GTFdxz !dxz 297 end do 298 hess(2,1,:)=hess(1,2,:) 299 hess(3,2,:)=hess(2,3,:) 300 hess(3,1,:)=hess(1,3,:) 301 end if 302 303 if (runtype>=5) then 304 !Calculate 3-order derivative for current GTF 305 ep2=ep*2D0 306 ep4=ep*4D0 307 epep4=ep2*ep2 308 epep8=epep4*2D0 309 !dxyz 310 a1=0D0 311 b1=0D0 312 c1=0D0 313 if (ix>=1) a1=ix*sftx**(ix-1) 314 if (iy>=1) b1=iy*sfty**(iy-1) 315 if (iz>=1) c1=iz*sftz**(iz-1) 316 a2=-ep2*sftx**(ix+1) 317 b2=-ep2*sfty**(iy+1) 318 c2=-ep2*sftz**(iz+1) 319 GTFdxyz=(a1+a2)*(b1+b2)*(c1+c2)*expterm 320 !dxyy,dxxy,dxxz,dxzz,dyzz,dyyz 321 atmp=0D0 322 btmp=0D0 323 ctmp=0D0 324 if (ix>=2) atmp=ix*(ix-1)*sftx**(ix-2) 325 if (iy>=2) btmp=iy*(iy-1)*sfty**(iy-2) 326 if (iz>=2) ctmp=iz*(iz-1)*sftz**(iz-2) 327 GTFdxyy=(a1+a2)*sftz**iz *expterm*(-ep4*iy*sfty**iy+epep4*sfty**(iy+2)+btmp-ep2*sfty**iy) !ok 328 GTFdxxy=(b1+b2)*sftz**iz *expterm*(-ep4*ix*sftx**ix+epep4*sftx**(ix+2)+atmp-ep2*sftx**ix) !=dyxx 329 GTFdxxz=(c1+c2)*sfty**iy *expterm*(-ep4*ix*sftx**ix+epep4*sftx**(ix+2)+atmp-ep2*sftx**ix) !=dzxx 330 GTFdxzz=(a1+a2)*sfty**iy *expterm*(-ep4*iz*sftz**iz+epep4*sftz**(iz+2)+ctmp-ep2*sftz**iz) 331 GTFdyzz=(b1+b2)*sftx**ix *expterm*(-ep4*iz*sftz**iz+epep4*sftz**(iz+2)+ctmp-ep2*sftz**iz) 332 GTFdyyz=(c1+c2)*sftx**ix *expterm*(-ep4*iy*sfty**iy+epep4*sfty**(iy+2)+btmp-ep2*sfty**iy) !=dzyy,ok 333 !dxxx,dyyy,dzzz 334 aatmp1=0D0 335 bbtmp1=0D0 336 cctmp1=0D0 337 if (ix>=1) aatmp1=ep2*ix*sftx**(ix-1) 338 if (iy>=1) bbtmp1=ep2*iy*sfty**(iy-1) 339 if (iz>=1) cctmp1=ep2*iz*sftz**(iz-1) 340 aatmp2=0D0 341 bbtmp2=0D0 342 cctmp2=0D0 343 if (ix>=2) aatmp2=ep2*ix*(ix-1)*sftx**(ix-1) 344 if (iy>=2) bbtmp2=ep2*iy*(iy-1)*sfty**(iy-1) 345 if (iz>=2) cctmp2=ep2*iz*(iz-1)*sftz**(iz-1) 346 aatmp3=0D0 347 bbtmp3=0D0 348 cctmp3=0D0 349 if (ix>=3) aatmp3=ix*(ix-1)*(ix-2)*sftx**(ix-3) 350 if (iy>=3) bbtmp3=iy*(iy-1)*(iy-2)*sfty**(iy-3) 351 if (iz>=3) cctmp3=iz*(iz-1)*(iz-2)*sftz**(iz-3) 352 GTFdxxx=sfty**iy*sftz**iz*expterm*( (-2*ix+ep2*sftx2-1)*(-epep4*sftx**(ix+1) + aatmp1) - aatmp2 + epep8*sftx**(ix+1) + aatmp3 ) 353 GTFdyyy=sftx**ix*sftz**iz*expterm*( (-2*iy+ep2*sfty2-1)*(-epep4*sfty**(iy+1) + bbtmp1) - bbtmp2 + epep8*sfty**(iy+1) + bbtmp3 ) 354 GTFdzzz=sfty**iy*sftx**ix*expterm*( (-2*iz+ep2*sftz2-1)*(-epep4*sftz**(iz+1) + cctmp1) - cctmp2 + epep8*sftz**(iz+1) + cctmp3 ) 355 356 !Calculate 3-order derivative tensor for orbital wavefunction 357 do imo=istart,iend 358 tens3(1,1,1,imo)=tens3(1,1,1,imo)+co(imo,j)*GTFdxxx !dxxx 359 tens3(2,2,2,imo)=tens3(2,2,2,imo)+co(imo,j)*GTFdyyy !dyyy 360 tens3(3,3,3,imo)=tens3(3,3,3,imo)+co(imo,j)*GTFdzzz !dzzz 361 tens3(1,2,2,imo)=tens3(1,2,2,imo)+co(imo,j)*GTFdxyy !dxyy* 362 tens3(1,1,2,imo)=tens3(1,1,2,imo)+co(imo,j)*GTFdxxy !dxxy* 363 tens3(1,1,3,imo)=tens3(1,1,3,imo)+co(imo,j)*GTFdxxz !dxxz* 364 tens3(1,3,3,imo)=tens3(1,3,3,imo)+co(imo,j)*GTFdxzz !dxzz* 365 tens3(2,3,3,imo)=tens3(2,3,3,imo)+co(imo,j)*GTFdyzz !dyzz* 366 tens3(2,2,3,imo)=tens3(2,2,3,imo)+co(imo,j)*GTFdyyz !dyyz* 367 tens3(1,2,3,imo)=tens3(1,2,3,imo)+co(imo,j)*GTFdxyz !dxyz 368 end do 369 tens3(1,2,1,:)=tens3(1,1,2,:) !dxyx=dxxy 370 tens3(1,3,1,:)=tens3(1,1,3,:) !dxzx=dxxz 371 tens3(1,3,2,:)=tens3(1,2,3,:) !dxzy=dxyz 372 tens3(2,1,1,:)=tens3(1,1,2,:) !dyxx=dxxy 373 tens3(2,1,2,:)=tens3(1,2,2,:) !dyxy=dxyy 374 tens3(2,1,3,:)=tens3(1,2,3,:) !dyxz=dxyz 375 tens3(2,2,1,:)=tens3(1,2,2,:) !dyyx=dxyy 376 tens3(2,3,1,:)=tens3(1,2,3,:) !dyzx=dxyz 377 tens3(2,3,2,:)=tens3(2,2,3,:) !dyzy=dyyz 378 tens3(3,1,1,:)=tens3(1,1,3,:) !dzxx=dxxz 379 tens3(3,1,2,:)=tens3(1,2,3,:) !dzxy=dxyz 380 tens3(3,1,3,:)=tens3(1,3,3,:) !dzxz=dxzz 381 tens3(3,2,1,:)=tens3(1,2,3,:) !dzyx=dxyz 382 tens3(3,2,2,:)=tens3(2,2,3,:) !dzyy=dyyz 383 tens3(3,2,3,:)=tens3(2,3,3,:) !dzyz=dyzz 384 tens3(3,3,1,:)=tens3(1,3,3,:) !dzzx=dxzz 385 tens3(3,3,2,:)=tens3(2,3,3,:) !dzzy=dyzz 386 end if !end runtype>=5 387 388 end if !end runtype>=3 389 end if !end runtype>=2 390end do 391end subroutine 392 393 394!!!----------- Calculate contribution from EDFs (recorded in wfx file) to density and corresponding derivatives (up to third-order) 395!Only S-type GTFs are supported 396! In wfx files, GTFs are used to expand core density 397! runtype=1: Only calculate rho, =2: rho+dx/dy/dz =3: rho+dx/dy/dz+dxx/dyy/dzz 398! =4: rho+dx/dy/dz+full Hessian =5: rho+dx/dy/dz+full Hessian+tens3 399subroutine EDFrho(runtype,x,y,z,value,grad,hess,tens3) 400integer runtype 401real*8 x,y,z,value 402real*8,optional :: grad(3),hess(3,3),tens3(3,3,3) 403value=0D0 404if (present(grad)) grad=0D0 405if (present(hess)) hess=0D0 406if (present(tens3)) tens3=0D0 407 408do i=1,nEDFprims 409 sftx=x-a(b_EDF(i)%center)%x 410 sfty=y-a(b_EDF(i)%center)%y 411 sftz=z-a(b_EDF(i)%center)%z 412 sftx2=sftx*sftx 413 sfty2=sfty*sfty 414 sftz2=sftz*sftz 415 rr=sftx2+sfty2+sftz2 416 ep=b_EDF(i)%exp 417 expterm=exp(-ep*rr) 418 value=value+CO_EDF(i)*expterm 419! write(11,"(i5,3D20.10)") i,value,CO_EDF(i),expterm 420 if (runtype>=2) then 421 tmp=2*CO_EDF(i)*expterm*ep 422 grad(1)=grad(1)-tmp*sftx 423 grad(2)=grad(2)-tmp*sfty 424 grad(3)=grad(3)-tmp*sftz 425 if (runtype>=3) then 426 hess(1,1)=hess(1,1)+tmp*(2*ep*sftx2-1) 427 hess(2,2)=hess(2,2)+tmp*(2*ep*sfty2-1) 428 hess(3,3)=hess(3,3)+tmp*(2*ep*sftz2-1) 429 if (runtype>=4) then 430 epep4=ep*ep*4 431 tmp2=CO_EDF(i)*epep4*expterm 432 hess(1,2)=hess(1,2)+tmp2*sftx*sfty 433 hess(1,3)=hess(1,3)+tmp2*sftx*sftz 434 hess(2,3)=hess(2,3)+tmp2*sfty*sftz 435 hess(2,1)=hess(1,2) 436 hess(3,1)=hess(1,3) 437 hess(3,2)=hess(2,3) 438 if (runtype>=5) then 439 tmp3=CO_EDF(i)*epep4*expterm 440 tens3(1,1,1)=tens3(1,1,1)+tmp3*sftx*(3-2*ep*sftx2) 441 tens3(2,2,2)=tens3(2,2,2)+tmp3*sfty*(3-2*ep*sfty2) 442 tens3(3,3,3)=tens3(3,3,3)+tmp3*sftz*(3-2*ep*sftz2) 443 tens3(1,2,2)=tens3(1,2,2)+tmp3*sftx*(1-2*ep*sfty2) 444 tens3(1,1,2)=tens3(1,1,2)+tmp3*sfty*(1-2*ep*sftx2) 445 tens3(1,1,3)=tens3(1,1,3)+tmp3*sftz*(1-2*ep*sftx2) 446 tens3(1,3,3)=tens3(1,3,3)+tmp3*sftx*(1-2*ep*sftz2) 447 tens3(2,3,3)=tens3(2,3,3)+tmp3*sfty*(1-2*ep*sftz2) 448 tens3(2,2,3)=tens3(2,2,3)+tmp3*sftz*(1-2*ep*sfty2) 449 tens3(1,2,3)=tens3(1,2,3)-CO_EDF(i)*8*ep**3*sftx*sfty*sftz*expterm 450 tens3(1,2,1)=tens3(1,1,2) !dxyx=dxxy 451 tens3(1,3,1)=tens3(1,1,3) !dxzx=dxxz 452 tens3(1,3,2)=tens3(1,2,3) !dxzy=dxyz 453 tens3(2,1,1)=tens3(1,1,2) !dyxx=dxxy 454 tens3(2,1,2)=tens3(1,2,2) !dyxy=dxyy 455 tens3(2,1,3)=tens3(1,2,3) !dyxz=dxyz 456 tens3(2,2,1)=tens3(1,2,2) !dyyx=dxyy 457 tens3(2,3,1)=tens3(1,2,3) !dyzx=dxyz 458 tens3(2,3,2)=tens3(2,2,3) !dyzy=dyyz 459 tens3(3,1,1)=tens3(1,1,3) !dzxx=dxxz 460 tens3(3,1,2)=tens3(1,2,3) !dzxy=dxyz 461 tens3(3,1,3)=tens3(1,3,3) !dzxz=dxzz 462 tens3(3,2,1)=tens3(1,2,3) !dzyx=dxyz 463 tens3(3,2,2)=tens3(2,2,3) !dzyy=dyyz 464 tens3(3,2,3)=tens3(2,3,3) !dzyz=dyzz 465 tens3(3,3,1)=tens3(1,3,3) !dzzx=dxzz 466 tens3(3,3,2)=tens3(2,3,3) !dzzy=dyzz 467 end if 468 end if 469 end if 470 end if 471end do 472end subroutine 473 474 475!!!------ A general routine used to calculate value, gradient and Hessian matrix at a given point for some real space functions 476! itype=1 Only calculate value and grad 477! itype=2 Calculate value, gradient and Hessian 478subroutine gencalchessmat(itype,ifunc,x,y,z,value,grad,hess) 479integer ifunc,itype 480real*8 x,y,z,value,grad(3),hess(3,3) 481real*8 gradaddx(3),gradminx(3),gradaddy(3),gradminy(3),gradaddz(3),gradminz(3) 482character selELFLOL*3 483diff=8D-4 484denom=2D0*diff 485if (ifunc==9) selELFLOL="ELF" 486if (ifunc==10) selELFLOL="LOL" 487 488!For a few functions whose both analytic gradient and Hessian are available, evaluate them and then return 489if (ifunc==1) then 490 call calchessmat_dens(itype,x,y,z,value,grad,hess) 491 return 492else if (ifunc==4) then 493 call calchessmat_mo(itype,iorbsel,x,y,z,value,grad,hess) 494 return 495end if 496 497!Calculate gradient 498!For other functions aside from above ones, analytical Hessian or even gradient hasn't been realized 499if (ifunc==3) then 500 call calchessmat_lapl(1,x,y,z,value,grad,hess) 501else if (ifunc==9.or.ifunc==10) then 502 if (ELFLOL_type==0) then !Analytical gradient for Becke's ELF/LOL 503 call calchessmat_ELF_LOL(1,x,y,z,value,grad,hess,selELFLOL) 504 else !Numerical gradient for other definition of ELF/LOL 505 value=ELF_LOL(x,y,z,selELFLOL) 506 xadd=ELF_LOL(x+diff,y,z,selELFLOL) 507 xmin=ELF_LOL(x-diff,y,z,selELFLOL) 508 yadd=ELF_LOL(x,y+diff,z,selELFLOL) 509 ymin=ELF_LOL(x,y-diff,z,selELFLOL) 510 zadd=ELF_LOL(x,y,z+diff,selELFLOL) 511 zmin=ELF_LOL(x,y,z-diff,selELFLOL) 512 grad(1)=(xadd-xmin)/denom 513 grad(2)=(yadd-ymin)/denom 514 grad(3)=(zadd-zmin)/denom 515 end if 516else if (ifunc==12) then 517 value=totesp(x,y,z) 518 xadd=totesp(x+diff,y,z) 519 xmin=totesp(x-diff,y,z) 520 yadd=totesp(x,y+diff,z) 521 ymin=totesp(x,y-diff,z) 522 zadd=totesp(x,y,z+diff) 523 zmin=totesp(x,y,z-diff) 524 grad(1)=(xadd-xmin)/denom 525 grad(2)=(yadd-ymin)/denom 526 grad(3)=(zadd-zmin)/denom 527else if (ifunc==100) then 528 value=userfunc(x,y,z) 529 xadd=userfunc(x+diff,y,z) 530 xmin=userfunc(x-diff,y,z) 531 yadd=userfunc(x,y+diff,z) 532 ymin=userfunc(x,y-diff,z) 533 zadd=userfunc(x,y,z+diff) 534 zmin=userfunc(x,y,z-diff) 535 grad(1)=(xadd-xmin)/denom 536 grad(2)=(yadd-ymin)/denom 537 grad(3)=(zadd-zmin)/denom 538end if 539 540!Calculate Hessian semi (namely based on analyical gradient) or pure (based on function value) numerically 541if (itype==2) then 542 if (ifunc==3) then !Use semi-analytical for Hessian of laplacian 543 call calchessmat_lapl(1,x+diff,y,z,tmpval,gradaddx,hess) 544 call calchessmat_lapl(1,x-diff,y,z,tmpval,gradminx,hess) 545 call calchessmat_lapl(1,x,y+diff,z,tmpval,gradaddy,hess) 546 call calchessmat_lapl(1,x,y-diff,z,tmpval,gradminy,hess) 547 call calchessmat_lapl(1,x,y,z+diff,tmpval,gradaddz,hess) 548 call calchessmat_lapl(1,x,y,z-diff,tmpval,gradminz,hess) 549 hess(1,1)=(gradaddx(1)-gradminx(1))/denom 550 hess(2,2)=(gradaddy(2)-gradminy(2))/denom 551 hess(3,3)=(gradaddz(3)-gradminz(3))/denom 552 hess(1,2)=(gradaddy(1)-gradminy(1))/denom 553 hess(2,3)=(gradaddz(2)-gradminz(2))/denom 554 hess(1,3)=(gradaddz(1)-gradminz(1))/denom 555 hess(2,1)=hess(1,2) 556 hess(3,2)=hess(2,3) 557 hess(3,1)=hess(1,3) 558 return !Don't do below procedures for generating pure numerical Hessian 559 else if (ifunc==9.or.ifunc==10) then 560 if (ELFLOL_type==0) then !Use semi-analytical for Hessian of Becke's ELF/LOL 561 call calchessmat_ELF_LOL(1,x+diff,y,z,tmpval,gradaddx,hess,selELFLOL) 562 call calchessmat_ELF_LOL(1,x-diff,y,z,tmpval,gradminx,hess,selELFLOL) 563 call calchessmat_ELF_LOL(1,x,y+diff,z,tmpval,gradaddy,hess,selELFLOL) 564 call calchessmat_ELF_LOL(1,x,y-diff,z,tmpval,gradminy,hess,selELFLOL) 565 call calchessmat_ELF_LOL(1,x,y,z+diff,tmpval,gradaddz,hess,selELFLOL) 566 call calchessmat_ELF_LOL(1,x,y,z-diff,tmpval,gradminz,hess,selELFLOL) 567 hess(1,1)=(gradaddx(1)-gradminx(1))/denom 568 hess(2,2)=(gradaddy(2)-gradminy(2))/denom 569 hess(3,3)=(gradaddz(3)-gradminz(3))/denom 570 hess(1,2)=(gradaddy(1)-gradminy(1))/denom 571 hess(2,3)=(gradaddz(2)-gradminz(2))/denom 572 hess(1,3)=(gradaddz(1)-gradminz(1))/denom 573 hess(2,1)=hess(1,2) 574 hess(3,2)=hess(2,3) 575 hess(3,1)=hess(1,3) 576 return 577 else !for other definition of ELF/LOL, use pure numerical Hessian 578 xaddxadd=ELF_LOL(x+2D0*diff,y,z,selELFLOL) 579 xminxmin=ELF_LOL(x-2D0*diff,y,z,selELFLOL) 580 yaddyadd=ELF_LOL(x,y+2D0*diff,z,selELFLOL) 581 yminymin=ELF_LOL(x,y-2D0*diff,z,selELFLOL) 582 zaddzadd=ELF_LOL(x,y,z+2D0*diff,selELFLOL) 583 zminzmin=ELF_LOL(x,y,z-2D0*diff,selELFLOL) 584 xaddyadd=ELF_LOL(x+diff,y+diff,z,selELFLOL) 585 xminyadd=ELF_LOL(x-diff,y+diff,z,selELFLOL) 586 xaddymin=ELF_LOL(x+diff,y-diff,z,selELFLOL) 587 xminymin=ELF_LOL(x-diff,y-diff,z,selELFLOL) 588 yaddzadd=ELF_LOL(x,y+diff,z+diff,selELFLOL) 589 yminzadd=ELF_LOL(x,y-diff,z+diff,selELFLOL) 590 yaddzmin=ELF_LOL(x,y+diff,z-diff,selELFLOL) 591 yminzmin=ELF_LOL(x,y-diff,z-diff,selELFLOL) 592 xaddzadd=ELF_LOL(x+diff,y,z+diff,selELFLOL) 593 xminzadd=ELF_LOL(x-diff,y,z+diff,selELFLOL) 594 xaddzmin=ELF_LOL(x+diff,y,z-diff,selELFLOL) 595 xminzmin=ELF_LOL(x-diff,y,z-diff,selELFLOL) 596 end if 597 else if (ifunc==12) then !pure numerical Hessian for total ESP 598 xaddxadd=totesp(x+2D0*diff,y,z) 599 xminxmin=totesp(x-2D0*diff,y,z) 600 yaddyadd=totesp(x,y+2D0*diff,z) 601 yminymin=totesp(x,y-2D0*diff,z) 602 zaddzadd=totesp(x,y,z+2D0*diff) 603 zminzmin=totesp(x,y,z-2D0*diff) 604 xaddyadd=totesp(x+diff,y+diff,z) 605 xminyadd=totesp(x-diff,y+diff,z) 606 xaddymin=totesp(x+diff,y-diff,z) 607 xminymin=totesp(x-diff,y-diff,z) 608 yaddzadd=totesp(x,y+diff,z+diff) 609 yminzadd=totesp(x,y-diff,z+diff) 610 yaddzmin=totesp(x,y+diff,z-diff) 611 yminzmin=totesp(x,y-diff,z-diff) 612 xaddzadd=totesp(x+diff,y,z+diff) 613 xminzadd=totesp(x-diff,y,z+diff) 614 xaddzmin=totesp(x+diff,y,z-diff) 615 xminzmin=totesp(x-diff,y,z-diff) 616 else if (ifunc==100) then !pure numerical Hessian for user-defined function 617 xaddxadd=userfunc(x+2D0*diff,y,z) 618 xminxmin=userfunc(x-2D0*diff,y,z) 619 yaddyadd=userfunc(x,y+2D0*diff,z) 620 yminymin=userfunc(x,y-2D0*diff,z) 621 zaddzadd=userfunc(x,y,z+2D0*diff) 622 zminzmin=userfunc(x,y,z-2D0*diff) 623 xaddyadd=userfunc(x+diff,y+diff,z) 624 xminyadd=userfunc(x-diff,y+diff,z) 625 xaddymin=userfunc(x+diff,y-diff,z) 626 xminymin=userfunc(x-diff,y-diff,z) 627 yaddzadd=userfunc(x,y+diff,z+diff) 628 yminzadd=userfunc(x,y-diff,z+diff) 629 yaddzmin=userfunc(x,y+diff,z-diff) 630 yminzmin=userfunc(x,y-diff,z-diff) 631 xaddzadd=userfunc(x+diff,y,z+diff) 632 xminzadd=userfunc(x-diff,y,z+diff) 633 xaddzmin=userfunc(x+diff,y,z-diff) 634 xminzmin=userfunc(x-diff,y,z-diff) 635 end if 636 !Collect above temporary data to evaluate pure numerical Hessian 637 gradx_yadd=(xaddyadd-xminyadd)/denom 638 gradx_ymin=(xaddymin-xminymin)/denom 639 grady_zadd=(yaddzadd-yminzadd)/denom 640 grady_zmin=(yaddzmin-yminzmin)/denom 641 gradx_zadd=(xaddzadd-xminzadd)/denom 642 gradx_zmin=(xaddzmin-xminzmin)/denom 643 hess(1,1)=(xaddxadd-2*value+xminxmin)/denom**2 644 hess(2,2)=(yaddyadd-2*value+yminymin)/denom**2 645 hess(3,3)=(zaddzadd-2*value+zminzmin)/denom**2 646 hess(1,2)=(gradx_yadd-gradx_ymin)/denom 647 hess(2,3)=(grady_zadd-grady_zmin)/denom 648 hess(1,3)=(gradx_zadd-gradx_zmin)/denom 649 hess(2,1)=hess(1,2) 650 hess(3,2)=hess(2,3) 651 hess(3,1)=hess(1,3) 652end if 653end subroutine 654 655 656 657!============================================================ 658!=================== Real space functions =================== 659!============================================================ 660 661 662!!!--------- User-defined function, the content is needed to be filled by users or selected by iuserfunc 663!The units should be given in a.u. 664real*8 function userfunc(x,y,z) 665real*8 x,y,z 666userfunc=1D0 !Default value. Note: default "iuserfunc" is 0 667!Below functions can be selected by "iuserfunc" parameter in settings.ini 668if (iuserfunc==-2) userfunc=calcprodens(x,y,z,0) !Promolecular density 669if (iuserfunc==-1) userfunc=linintp3d(x,y,z,1) !The function value evaluated by trilinear interpolation from cubmat 670if (iuserfunc==1) userfunc=fspindens(x,y,z,'a') !Alpha density 671if (iuserfunc==2) userfunc=fspindens(x,y,z,'b') !Beta density 672if (iuserfunc==3) userfunc=(x*x+y*y+z*z)*fdens(x,y,z) !Integrand of electronic spatial extent <R**2> 673if (iuserfunc==4) userfunc=weizpot(x,y,z) !Weizsacker potential 674if (iuserfunc==5) userfunc=weizsacker(x,y,z) !Integrand of weizsacker functional 675if (iuserfunc==6) userfunc=4*pi*fdens(x,y,z)*(x*x+y*y+z*z) !Radial distribution function (assume that density is sphericalized) 676if (iuserfunc==7) userfunc=2D0/3D0*lagkin(x,y,z,0)/fdens(x,y,z) !Local Temperature(Kelvin), PNAS,81,8028 677if (iuserfunc==8) userfunc=totesp(x,y,z)/fdens(x,y,z) !Average local electrostatic potential, useful to exhibit atomic shell structure, see Chapter 8 of Theoretical Aspects of Chemical Reactivity 678if (iuserfunc==9) userfunc=fdens(x,y,z)/nelec !Shape function 679if (iuserfunc==10) userfunc=-Hamkin(x,y,z,0)-lagkin(x,y,z,0) !Potential energy density, also known as virial field 680if (iuserfunc==11) userfunc=-Hamkin(x,y,z,0) !Energy density 681if (iuserfunc==12) userfunc=-nucesp(x,y,z)*fdens(x,y,z) !Local nuclear attraction potential energy 682if (iuserfunc==13) userfunc=lagkin(x,y,z,0)/fdens(x,y,z) !This quantity at bond critical point is useful to discriminate covalent bonding and close-shell interaction 683if (iuserfunc==14) userfunc=eleesp(x,y,z) !Electrostatic potential from electrons 684if (iuserfunc==15) userfunc=fdens(x,y,z)/flapl(x,y,z,'t') !Bond metallicity 685if (iuserfunc==16) userfunc=36*(3*pi*pi)**(2D0/3D0)/5D0*fdens(x,y,z)**(5D0/3D0)/flapl(x,y,z,'t') !Dimensionless bond metallicity 686if (iuserfunc==17) userfunc=-Hamkin(x,y,z,0)/fdens(x,y,z) !Energy density per electron 687if (iuserfunc==18) then !Region of Slow Electrons (RoSE), defined in Chem. Phys. Lett., 582, 144 (2013) 688 rho=fdens(x,y,z) 689 if (wfntype==0.or.wfntype==3) then !close shell cases 690 Dh=2.871234000D0*rho**(5.0D0/3.0D0) 691 else if (wfntype==1.or.wfntype==2.or.wfntype==4) then !Open shell cases 692 rhospin=fspindens(x,y,z,'s') !rhospin=rhoa-rhob, rho=rhoa+rhob 693 rhoa=(rhospin+rho)/2D0 694 rhob=(rho-rhospin)/2D0 695 Dh=4.557799872D0*(rhoa**(5.0D0/3.0D0)+rhob**(5.0D0/3.0D0)) !kinetic energy of HEG 696 end if 697 G=Lagkin(x,y,z,0) 698 userfunc=(Dh-G)/(Dh+G) 699end if 700if (iuserfunc==19) userfunc=SEDD(x,y,z) !SEDD 701if (iuserfunc==20) userfunc=DORI(x,y,z) !DORI 702if (iuserfunc==21) userfunc=-x*fdens(x,y,z) !Integrand of X component of electric dipole moment 703if (iuserfunc==22) userfunc=-y*fdens(x,y,z) !Integrand of Y component of electric dipole moment 704if (iuserfunc==23) userfunc=-z*fdens(x,y,z) !Integrand of Z component of electric dipole moment 705if (iuserfunc==24) userfunc=linrespkernel(x,y,z) !Approximate form of DFT linear response kernel for close-shell 706if (iuserfunc==25) userfunc=fgrad(x,y,z,'t')/fdens(x,y,z)/2D0 !Magnitude of fluctuation of the electronic momentum 707if (iuserfunc==26) userfunc=2.871234D0*rho**(5D0/3D0) !Thomas-Fermi kinetic energy density 708if (iuserfunc==27) userfunc=loceleaff(x,y,z) !Local electron affinity 709if (iuserfunc==28) userfunc=(avglocion(x,y,z)+loceleaff(x,y,z))/2 !Local Mulliken electronegativity 710if (iuserfunc==29) userfunc=(avglocion(x,y,z)-loceleaff(x,y,z))/2 !Local hardness 711if (iuserfunc==30) userfunc=densellip(x,y,z,1) !Ellipticity of electron density 712if (iuserfunc==31) userfunc=densellip(x,y,z,2) !eta index, Angew. Chem. Int. Ed., 53, 2766-2770 (2014) 713if (iuserfunc==32) userfunc=densellip(x,y,z,2)-1 !Modified eta index 714if (iuserfunc==33) userfunc=PAEM(x,y,z,1) !PAEM, potential acting on one electron in a molecule, defined by Zhongzhi Yang 715if (iuserfunc==34) userfunc=PAEM(x,y,z,2) !The same as 33, but using DFT XC potential directly rather than evaluating the XC potential based on pair density 716if (iuserfunc==35) then !|V(r)|/G(r) 717 tmpval=lagkin(x,y,z,0) 718 userfunc=abs(-Hamkin(x,y,z,0)-tmpval)/tmpval 719end if 720if (iuserfunc==36) then !On-top pair density, i.e. r1=r2 case of pair density. paircorrtype affects result 721 pairfunctypeold=pairfunctype 722 pairfunctype=12 723 userfunc=pairfunc(x,y,z,x,y,z) 724 pairfunctype=pairfunctypeold 725end if 726if (iuserfunc==38) userfunc=Ang_rhoeigvec_ple(x,y,z,2) !The angle between the second eigenvector of rho and the plane defined by option 4 of main function 1000 727if (iuserfunc==39) userfunc=totespskip(x,y,z,iskipnuc) !ESP without contribution of nuclues "iskipnuc" 728if (iuserfunc==40) userfunc=weizsacker(x,y,z) !Steric energy 729if (iuserfunc==41) userfunc=stericpot(x,y,z) !Steric potential 730if (iuserfunc==42) userfunc=stericcharge(x,y,z) !Steric charge 731if (iuserfunc==43) userfunc=stericforce(x,y,z) !The magnitude of steric force 732if (iuserfunc==44) userfunc=stericpot_damp(x,y,z) !Steric potential with damping function to a given constant value 733if (iuserfunc==45) userfunc=stericforce_damp(x,y,z) !Steric force based on damped potential 734if (iuserfunc==46) userfunc=stericforce_directdamp(x,y,z) !Steric force directly damped to zero 735if (iuserfunc==50) userfunc=infoentro(2,x,y,z) !Shannon entropy density, see JCP,126,191107 for example 736if (iuserfunc==51) userfunc=Fisherinfo(1,x,y,z) !Fisher information density, see JCP,126,191107 for example 737if (iuserfunc==52) userfunc=Fisherinfo(2,x,y,z) !Second Fisher information density, see JCP,126,191107 for derivation 738if (iuserfunc==53) userfunc=Ghoshentro(x,y,z,1) !Ghosh entropy density with G(r) as kinetic energy density, PNAS,81,8028 739if (iuserfunc==54) userfunc=Ghoshentro(x,y,z,2) !Ghosh entropy density with G(r)-der2rho/8 as kinetic energy density, exactly corresponds to Eq.22 in PNAS,81,8028 740if (iuserfunc==55) userfunc=fdens(x,y,z)**2 !Integrand of quadratic form of Renyi entropy 741if (iuserfunc==56) userfunc=fdens(x,y,z)**3 !Integrand of cubic form of Renyi entropy 742if (iuserfunc==60) userfunc=paulipot(x,y,z) !Pauli potential, Comp. Theor. Chem., 1006, 92-99 743if (iuserfunc==61) userfunc=pauliforce(x,y,z) !The magnitude of Pauli force 744if (iuserfunc==62) userfunc=paulicharge(x,y,z) !Pauli charge 745if (iuserfunc==63) userfunc=quantumpot(x,y,z) !Quantum potential 746if (iuserfunc==64) userfunc=quantumforce(x,y,z) !The magnitude of quantum force 747if (iuserfunc==65) userfunc=quantumcharge(x,y,z) !Quantum charge 748if (iuserfunc==66) userfunc=elestatforce(x,y,z) !The magnitude of electrostatic force 749if (iuserfunc==67) userfunc=elestatcharge(x,y,z) !Electrostatic charge 750if (iuserfunc==70) userfunc=4.5D0*fdens(x,y,z)**2/lagkin(x,y,z,0) !Phase-space-defined Fisher information density 751if (iuserfunc==71) userfunc=elemomdens(x,y,z,1) !X component of electron linear momentum density in 3D representation 752if (iuserfunc==72) userfunc=elemomdens(x,y,z,2) !Y component of electron linear momentum density in 3D representation 753if (iuserfunc==73) userfunc=elemomdens(x,y,z,3) !Z component of electron linear momentum density in 3D representation 754if (iuserfunc==74) userfunc=elemomdens(x,y,z,0) !Magnitude of electron linear momentum density in 3D representation 755if (iuserfunc==75) userfunc=magmomdens(x,y,z,1) !X component of magnetic dipole moment density 756if (iuserfunc==76) userfunc=magmomdens(x,y,z,2) !Y component of magnetic dipole moment density 757if (iuserfunc==77) userfunc=magmomdens(x,y,z,3) !Z component of magnetic dipole moment density 758if (iuserfunc==78) userfunc=magmomdens(x,y,z,0) !Magnitude of magnetic dipole moment density 759if (iuserfunc==79) userfunc=energydens_grdn(x,y,z) !Gradient norm of energy density 760if (iuserfunc==80) userfunc=energydens_lapl(x,y,z) !Laplacian of energy density 761if (iuserfunc==81) userfunc=hamkin(x,y,z,1) !X component of Hamiltonian kinetic energy density 762if (iuserfunc==82) userfunc=hamkin(x,y,z,2) !Y component of Hamiltonian kinetic energy density 763if (iuserfunc==83) userfunc=hamkin(x,y,z,3) !Z component of Hamiltonian kinetic energy density 764if (iuserfunc==84) userfunc=Lagkin(x,y,z,1) !X component of Lagrangian kinetic energy density 765if (iuserfunc==85) userfunc=Lagkin(x,y,z,2) !Y component of Lagrangian kinetic energy density 766if (iuserfunc==86) userfunc=Lagkin(x,y,z,3) !Z component of Lagrangian kinetic energy density 767if (iuserfunc==87) userfunc=localcorr(x,y,z,1) !Local total electron correlation function 768if (iuserfunc==88) userfunc=localcorr(x,y,z,2) !Local dynamic electron correlation function 769if (iuserfunc==89) userfunc=localcorr(x,y,z,3) !Local nondynamic electron correlation function 770if (iuserfunc==90) then 771 tmpELF=ELF_LOL(x,y,z,"ELF") 772 userfunc=tmpELF*tmpELF*(x*x+y*y+z*z) 773else if (iuserfunc==91) then 774 userfunc=tmpELF*tmpELF 775end if 776if (iuserfunc==100) userfunc=fdens(x,y,z)**2 !Disequilibrium (also known as semi-similarity), DOI: 10.1002/qua.24510 777if (iuserfunc==101) then !Positive part of ESP 778 userfunc=totesp(x,y,z) 779 if (userfunc<0D0) userfunc=0D0 780else if (iuserfunc==102) then !Negative part of ESP 781 userfunc=totesp(x,y,z) 782 if (userfunc>0D0) userfunc=0D0 783end if 784if (iuserfunc>=802.and.iuserfunc<=807) userfunc=funcvalLSB(x,y,z,iuserfunc-800) 785if (iuserfunc>=812.and.iuserfunc<=817) userfunc=1/funcvalLSB(x,y,z,iuserfunc-810) 786if (iuserfunc==900) userfunc=x !X coordinate 787if (iuserfunc==901) userfunc=y !Y coordinate 788if (iuserfunc==902) userfunc=z !Z coordinate 789if (iuserfunc==1000) userfunc=DFTxcfunc(x,y,z) !Various kinds of DFT exchange-correlation functions 790if (iuserfunc==1100) userfunc=DFTxcpot(x,y,z) !Various kinds of DFT exchange-correlation potentials 791!Below are other examples 792! userfunc=hamkin(x,y,z,3)-0.5D0*(hamkin(x,y,z,1)+hamkin(x,y,z,2)) !Anisotropy of Hamiltonian kinetic energy in Z, namely K_Z-0.5*(K_X+K_Y) 793! userfunc=-x*y*fdens(x,y,z) !Integrand of XY component of electric quadrupole moment 794! userfunc=-x*y*z*fdens(x,y,z)*au2debye*b2a**2 !Integrand of XYZ component of electric octapole moment in Debye-Ang**2 795end function 796 797 798!!!----------- Output orbital wavefunction value at a given point (fmo=function for outputting MO) 799real*8 function fmo(x,y,z,id) 800real*8 x,y,z,orbval(nmo) 801integer id 802call orbderv(1,id,id,x,y,z,orbval) 803fmo=orbval(id) 804end function 805 806!!!----- Calculate orbital wavefunction Hessian matrix at x,y,z, store to hess, output value and gradient vector at same time 807!itype=1 Only calculate value and gradient, not Hessian 808!itype=2 Calculate value, gradient and Hessian 809subroutine calchessmat_mo(itype,id,x,y,z,value,grad,hess) 810integer id,itype 811real*8 x,y,z,grad(3),hess(3,3),value,wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo) 812if (itype==1) call orbderv(2,id,id,x,y,z,wfnval,wfnderv,wfnhess) 813if (itype==2) call orbderv(4,id,id,x,y,z,wfnval,wfnderv,wfnhess) 814value=wfnval(id) 815grad=wfnderv(:,id) 816hess=wfnhess(:,:,id) 817end subroutine 818 819 820!--------Output electron density at a point 821real*8 function fdens(x,y,z) 822real*8 x,y,z,wfnval(nmo) 823call orbderv(1,1,nmo,x,y,z,wfnval) 824fdens=0D0 825do i=1,nmo 826 fdens=fdens+MOocc(i)*wfnval(i)**2 827end do 828! Add in contribution of Electron density function 829if (allocated(b_EDF)) then 830 call EDFrho(1,x,y,z,EDFdens) 831 fdens=fdens+EDFdens 832end if 833! if (fdens>0.5D0) fdens=0 834end function 835 836 837!!!------------------------- Output spin or Alpha or Beta electron density at a point 838!itype='s' output spin density, ='a' output alpha density, ='b' output beta density 839real*8 function fspindens(x,y,z,itype) 840real*8 :: x,y,z,wfnval(nmo) 841character itype 842call orbderv(1,1,nmo,x,y,z,wfnval) 843adens=0.0D0 844bdens=0.0D0 845do i=1,nmo 846 if (MOtype(i)==1) then 847 adens=adens+MOocc(i)*wfnval(i)**2 848 else if (MOtype(i)==2) then 849 bdens=bdens+MOocc(i)*wfnval(i)**2 850 else if (MOtype(i)==0) then 851 adens=adens+MOocc(i)/2D0*wfnval(i)**2 852 bdens=bdens+MOocc(i)/2D0*wfnval(i)**2 853 end if 854end do 855if (itype=='s') then 856 fspindens=adens-bdens 857 if (ipolarpara==1) fspindens=fspindens/(adens+bdens) 858else if (itype=='a') then 859 fspindens=adens 860else if (itype=='b') then 861 fspindens=bdens 862end if 863end function 864 865 866!!!------------------------- Output gradient of rho and RDG(reduced density gradient) at a point 867!label=x/y/z output 1-order derivation of x/y/z, =t get norm, =r get RDG, =s get |der_rho|/rho^(4/3) 868real*8 function fgrad(x,y,z,label) 869real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo),gradrho(3),EDFgrad(3),sumgrad2 870character label 871call orbderv(2,1,nmo,x,y,z,wfnval,wfnderv) 872rho=0D0 873gradrho=0D0 874do i=1,nmo 875 rho=rho+MOocc(i)*wfnval(i)**2 876 gradrho(:)=gradrho(:)+MOocc(i)*wfnval(i)*wfnderv(:,i) 877end do 878gradrho=2*gradrho 879! Add in contribution of Electron density function 880if (allocated(b_EDF)) then 881 call EDFrho(2,x,y,z,EDFdens,EDFgrad) 882 rho=rho+EDFdens 883 gradrho=gradrho+EDFgrad 884end if 885if (label=='x') then 886 fgrad=gradrho(1) 887else if (label=='y') then 888 fgrad=gradrho(2) 889else if (label=='z') then 890 fgrad=gradrho(3) 891else if (label=='t') then 892 fgrad=dsqrt( sum(gradrho(:)**2) ) 893else if (label=='r') then 894 sumgrad2=sum(gradrho(:)**2) 895 if (RDG_maxrho/=0.0D0.and.rho>=RDG_maxrho) then 896 fgrad=100D0 897!This occurs at distant region when exponent cutoff is used, the actual value should be very large. In order to avoid denominator become zero, we set it artifically to a big value 898 else if (sumgrad2==0D0.or.rho==0D0) then 899 RDG=999D0 900 else 901 fgrad=0.161620459673995D0*dsqrt(sumgrad2)/rho**(4.0D0/3.0D0) !0.161620459673995D0=1/(2*(3*pi**2)**(1/3)) 902 end if 903else if (label=='s') then 904 if (rho==0D0) then 905 fgrad=0 906 else 907 fgrad=dsqrt(sum(gradrho(:)**2))/rho**(4.0D0/3.0D0) 908 end if 909end if 910end function 911 912 913!!--- Simultaneously generate electron density, gradient norm for alpha and beta electrons, as well as dot product between grada and gradb 914!---- Mainly used to evalute DFT functional. EDF is not taken into account 915!adens/bdens/tdens means the density of alpha/beta/total density, similar for *grad 916subroutine gendensgradab(x,y,z,adens,bdens,tdens,agrad,bgrad,tgrad,abgrad) 917real*8 x,y,z,adens,bdens,agrad,bgrad,abgrad,wfnval(nmo),wfnderv(3,nmo),gradrhoa(3),gradrhob(3),gradrhot(3),tmparr(3) 918call orbderv(2,1,nmo,x,y,z,wfnval,wfnderv) 919adens=0D0 920bdens=0D0 921gradrhoa=0D0 922gradrhob=0D0 923do i=1,nmo 924 if (MOtype(i)==1) then 925 adens=adens+MOocc(i)*wfnval(i)**2 926 gradrhoa(:)=gradrhoa(:)+MOocc(i)*wfnval(i)*wfnderv(:,i) 927 else if (MOtype(i)==2) then 928 bdens=bdens+MOocc(i)*wfnval(i)**2 929 gradrhob(:)=gradrhob(:)+MOocc(i)*wfnval(i)*wfnderv(:,i) 930 else if (MOtype(i)==0) then 931 tmpval=MOocc(i)/2D0*wfnval(i)**2 932 adens=adens+tmpval 933 bdens=bdens+tmpval 934 tmparr(:)=MOocc(i)/2D0*wfnval(i)*wfnderv(:,i) 935 gradrhoa(:)=gradrhoa(:)+tmparr(:) 936 gradrhob(:)=gradrhob(:)+tmparr(:) 937 end if 938end do 939tdens=adens+bdens 940gradrhoa=gradrhoa*2 941gradrhob=gradrhob*2 942gradrhot=gradrhoa+gradrhob 943agrad=dsqrt(sum(gradrhoa**2)) 944bgrad=dsqrt(sum(gradrhob**2)) 945tgrad=dsqrt(sum(gradrhot**2)) 946abgrad=sum(gradrhoa*gradrhob) 947end subroutine 948 949 950!!!------------------------- Output Laplacian of electron density at a point 951!label=x/y/z output 2-order derivative of electron density respect to xx/yy/zz; & 952!=t get their summing; =s get der2rho/rho^(5/3), which is used LSB's project 953real*8 function flapl(x,y,z,label) 954real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo),laplx,laply,laplz,EDFgrad(3),EDFhess(3,3) 955character label 956call orbderv(3,1,nmo,x,y,z,wfnval,wfnderv,wfnhess) 957laplx=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)**2 + wfnval(1:nmo)*wfnhess(1,1,1:nmo) ) ) 958laply=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)**2 + wfnval(1:nmo)*wfnhess(2,2,1:nmo) ) ) 959laplz=2*sum( MOocc(1:nmo)*( wfnderv(3,1:nmo)**2 + wfnval(1:nmo)*wfnhess(3,3,1:nmo) ) ) 960!Add in contribution of electron density function, assume EDFs are S type 961if (allocated(b_EDF)) then 962 call EDFrho(3,x,y,z,EDFdens,EDFgrad,EDFhess) 963 laplx=laplx+EDFhess(1,1) 964 laply=laply+EDFhess(2,2) 965 laplz=laplz+EDFhess(3,3) 966end if 967if (label=='t') then 968 flapl=laplx+laply+laplz 969 flapl=flapl*laplfac !laplfac is an external variable 970else if (label=='x') then 971 flapl=laplx 972else if (label=='y') then 973 flapl=laply 974else if (label=='z') then 975 flapl=laplz 976else if (label=='s') then 977 dens=sum(MOocc(1:nmo)*wfnval(1:nmo)**2) 978 if (allocated(b_EDF)) dens=dens+EDFdens 979 flapl=(laplx+laply+laplz)/dens**(5D0/3D0) 980end if 981end function 982 983 984!!!----- Calculate electron density, its gradient and Hessian matrix 985!itype=1 Only calculate value and grad, not Hessian 986!itype=2 calculate value, gradient and Hessian 987subroutine calchessmat_dens(itype,x,y,z,elerho,elegrad,elehess) 988real*8 x,y,z,elerho,elegrad(3),elehess(3,3),wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo),EDFgrad(3),EDFhess(3,3) 989integer itype 990call orbderv(4,1,nmo,x,y,z,wfnval,wfnderv,wfnhess) 991elerho=sum( MOocc(1:nmo)*wfnval(1:nmo)*wfnval(1:nmo) ) 992do itmp=1,3 993 elegrad(itmp)=2*sum( MOocc(1:nmo)*wfnval(1:nmo)*wfnderv(itmp,1:nmo) ) 994end do 995if (itype==2) then 996 elehess(1,1)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)**2 + wfnval(1:nmo)*wfnhess(1,1,1:nmo) ) ) 997 elehess(2,2)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)**2 + wfnval(1:nmo)*wfnhess(2,2,1:nmo) ) ) 998 elehess(3,3)=2*sum( MOocc(1:nmo)*( wfnderv(3,1:nmo)**2 + wfnval(1:nmo)*wfnhess(3,3,1:nmo) ) ) 999 elehess(1,2)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(2,1:nmo)+wfnhess(1,2,1:nmo)*wfnval(1:nmo) ) ) 1000 elehess(2,3)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)*wfnderv(3,1:nmo)+wfnhess(2,3,1:nmo)*wfnval(1:nmo) ) ) 1001 elehess(1,3)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(3,1:nmo)+wfnhess(1,3,1:nmo)*wfnval(1:nmo) ) ) 1002 elehess(2,1)=elehess(1,2) 1003 elehess(3,2)=elehess(2,3) 1004 elehess(3,1)=elehess(1,3) 1005end if 1006 1007!Add in contribution of electron density function, assume EDFs are S type 1008if (allocated(b_EDF)) then 1009 call EDFrho(4,x,y,z,EDFdens,EDFgrad,EDFhess) 1010 elerho=elerho+EDFdens 1011 elegrad=elegrad+EDFgrad 1012 elehess=elehess+EDFhess 1013end if 1014end subroutine 1015 1016 1017!!------------- Calculate Laplacian of electron density, its gradient and Hessian matrix 1018!itype=1 calculate value, gradient 1019!itype=2 calculate value, gradient and Hessian (Not available) 1020subroutine calchessmat_lapl(itype,x,y,z,value,grad,hess) 1021use util 1022real*8 x,y,z,value,grad(3),hess(3,3) 1023real*8 wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo),wfntens3(3,3,3,nmo),rhotens3(3,3,3) 1024real*8 EDFgrad(3),EDFhess(3,3),EDFtens3(3,3,3) 1025integer itype 1026!Numerically verify 3-order derivative of orbital wavefunction 1027! diff=1D-5 1028! call orbderv(5,1,nmo,x+diff,y,z,wfnval,wfnderv,wfnhess,wfntens3) 1029! t1=wfnhess(3,3,1) 1030! call orbderv(5,1,nmo,x-diff,y,z,wfnval,wfnderv,wfnhess,wfntens3) 1031! t2=wfnhess(3,3,1) 1032! write(*,*) (t1-t2)/(2*diff) 1033 1034call orbderv(5,1,nmo,x,y,z,wfnval,wfnderv,wfnhess,wfntens3) 1035rhotens3=0D0 1036dxx=0D0 1037dyy=0D0 1038dzz=0D0 1039do i=1,nmo 1040 dxx=dxx+MOocc(i)*(wfnderv(1,i)**2+wfnval(i)*wfnhess(1,1,i)) 1041 dyy=dyy+MOocc(i)*(wfnderv(2,i)**2+wfnval(i)*wfnhess(2,2,i)) 1042 dzz=dzz+MOocc(i)*(wfnderv(3,i)**2+wfnval(i)*wfnhess(3,3,i)) 1043 rhotens3(1,1,1)=rhotens3(1,1,1)+MOocc(i)*( 3*wfnderv(1,i)*wfnhess(1,1,i)+wfnval(i)*wfntens3(1,1,1,i) ) 1044 rhotens3(2,2,2)=rhotens3(2,2,2)+MOocc(i)*( 3*wfnderv(2,i)*wfnhess(2,2,i)+wfnval(i)*wfntens3(2,2,2,i) ) 1045 rhotens3(3,3,3)=rhotens3(3,3,3)+MOocc(i)*( 3*wfnderv(3,i)*wfnhess(3,3,i)+wfnval(i)*wfntens3(3,3,3,i) ) 1046 rhotens3(1,1,2)=rhotens3(1,1,2)+MOocc(i)*( 2*wfnderv(1,i)*wfnhess(1,2,i)+wfnderv(2,i)*wfnhess(1,1,i)+wfnval(i)*wfntens3(1,1,2,i) ) 1047 rhotens3(1,1,3)=rhotens3(1,1,3)+MOocc(i)*( 2*wfnderv(1,i)*wfnhess(1,3,i)+wfnderv(3,i)*wfnhess(1,1,i)+wfnval(i)*wfntens3(1,1,3,i) ) 1048 rhotens3(2,2,3)=rhotens3(2,2,3)+MOocc(i)*( 2*wfnderv(2,i)*wfnhess(2,3,i)+wfnderv(3,i)*wfnhess(2,2,i)+wfnval(i)*wfntens3(2,2,3,i) ) 1049 rhotens3(1,2,2)=rhotens3(1,2,2)+MOocc(i)*( 2*wfnderv(2,i)*wfnhess(2,1,i)+wfnderv(1,i)*wfnhess(2,2,i)+wfnval(i)*wfntens3(2,2,1,i) ) !=2,2,1 exchange 1<->2 from (1,1,2) to derive this 1050 rhotens3(1,3,3)=rhotens3(1,3,3)+MOocc(i)*( 2*wfnderv(3,i)*wfnhess(3,1,i)+wfnderv(1,i)*wfnhess(3,3,i)+wfnval(i)*wfntens3(3,3,1,i) ) !2.758568947939382D-002 1051 rhotens3(2,3,3)=rhotens3(2,3,3)+MOocc(i)*( 2*wfnderv(3,i)*wfnhess(3,2,i)+wfnderv(2,i)*wfnhess(3,3,i)+wfnval(i)*wfntens3(3,3,2,i) ) 1052! write(*,*) "A",wfnhess(3,1,i),wfntens3(1,3,3,i) 1053end do 1054dxx=2D0*dxx 1055dyy=2D0*dyy 1056dzz=2D0*dzz 1057value=laplfac*(dxx+dyy+dzz) 1058rhotens3=rhotens3*2D0*laplfac 1059grad(1)=rhotens3(1,1,1)+rhotens3(1,2,2)+rhotens3(1,3,3) 1060grad(2)=rhotens3(1,1,2)+rhotens3(2,2,2)+rhotens3(2,3,3) 1061grad(3)=rhotens3(1,1,3)+rhotens3(2,2,3)+rhotens3(3,3,3) 1062 1063! diff=1D-5 1064!Check of flapldx,dy,dz is correct! 1065! write(*,*) grad(:) 1066! difflapldx=(flapl(x+diff,y,z,'t')-flapl(x-diff,y,z,'t'))/(2*diff) 1067! difflapldy=(flapl(x,y+diff,z,'t')-flapl(x,y-diff,z,'t'))/(2*diff) 1068! difflapldz=(flapl(x,y,z+diff,'t')-flapl(x,y,z-diff,'t'))/(2*diff) 1069! write(*,*) difflapldx,difflapldy,difflapldz 1070 1071!Check deviation between analytic and numerical solution 1072! write(*,*) rhotens3(1,1,1),rhotens3(1,2,2),rhotens3(1,3,3) 1073! diffrhodxxx=(flapl(x+diff,y,z,'x')-flapl(x-diff,y,z,'x'))/(2*diff) 1074! diffrhodxyy=(flapl(x+diff,y,z,'y')-flapl(x-diff,y,z,'y'))/(2*diff) 1075! diffrhodxzz=(flapl(x+diff,y,z,'z')-flapl(x-diff,y,z,'z'))/(2*diff) 1076! write(*,*) diffrhodxxx,diffrhodxyy,diffrhodxzz 1077! write(*,*) 1078 1079!Check diagonal term with finite difference 1080! write(*,*) rhotens3(1,1,1),rhotens3(2,2,2),rhotens3(3,3,3) 1081! diffrhodxxx=(flapl(x+diff,y,z,'x')-flapl(x-diff,y,z,'x'))/(2*diff) 1082! diffrhodyyy=(flapl(x,y+diff,z,'y')-flapl(x,y-diff,z,'y'))/(2*diff) 1083! diffrhodzzz=(flapl(x,y,z+diff,'z')-flapl(x,y,z-diff,'z'))/(2*diff) 1084! write(*,*) diffrhodxxx,diffrhodyyy,diffrhodzzz 1085 1086if (allocated(b_EDF)) then 1087 call EDFrho(5,x,y,z,EDFdens,EDFgrad,EDFhess,EDFtens3) 1088 grad(1)=grad(1)+EDFtens3(1,1,1)+EDFtens3(1,2,2)+EDFtens3(1,3,3) 1089 grad(2)=grad(2)+EDFtens3(1,1,2)+EDFtens3(2,2,2)+EDFtens3(2,3,3) 1090 grad(3)=grad(3)+EDFtens3(1,1,3)+EDFtens3(2,2,3)+EDFtens3(3,3,3) 1091end if 1092 1093!Don't consider Laplacian currently 1094if (itype==2) then !Calculate Hessian of laplacian 1095end if 1096end subroutine 1097 1098 1099!!!------------------------- Output Lagrangian kinetic G(r) at a point. idir=0/1/2/3 means total/x/y/z kinetic energy 1100real*8 function Lagkin(x,y,z,idir) 1101real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo) 1102integer idir 1103call orbderv(2,1,nmo,x,y,z,wfnval,wfnderv) 1104lagkin=0D0 1105if (idir==0) then 1106 do imo=1,nmo 1107 lagkin=lagkin+MOocc(imo)*sum(wfnderv(:,imo)**2) 1108 end do 1109else 1110 do imo=1,nmo 1111 lagkin=lagkin+MOocc(imo)*wfnderv(idir,imo)**2 1112 end do 1113end if 1114lagkin=lagkin/2D0 1115end function 1116 1117 1118!!------------- Output Hamiltonian kinetic K(r) at a point. idir=0/1/2/3 means total/X/Y/Z kinetic energy 1119real*8 function Hamkin(x,y,z,idir) 1120real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo) 1121integer idir 1122call orbderv(3,1,nmo,x,y,z,wfnval,wfnderv,wfnhess) 1123if (idir==0) then 1124 hamx=sum( MOocc(1:nmo)*wfnhess(1,1,1:nmo)*wfnval(1:nmo) ) 1125 hamy=sum( MOocc(1:nmo)*wfnhess(2,2,1:nmo)*wfnval(1:nmo) ) 1126 hamz=sum( MOocc(1:nmo)*wfnhess(3,3,1:nmo)*wfnval(1:nmo) ) 1127 Hamkin=hamx+hamy+hamz 1128else if (idir==1) then 1129 Hamkin=sum( MOocc(1:nmo)*wfnhess(1,1,1:nmo)*wfnval(1:nmo) ) 1130else if (idir==2) then 1131 Hamkin=sum( MOocc(1:nmo)*wfnhess(2,2,1:nmo)*wfnval(1:nmo) ) 1132else if (idir==3) then 1133 Hamkin=sum( MOocc(1:nmo)*wfnhess(3,3,1:nmo)*wfnval(1:nmo) ) 1134end if 1135Hamkin=-Hamkin/2D0 1136end function 1137 1138 1139!!------- Output analytic gradient vector of energy density (i.e. negative of K(r)) 1140subroutine energydens_grad(x,y,z,grad) 1141real*8 x,y,z,grad(3),wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo),wfntens(3,3,3,nmo) 1142call orbderv(5,1,nmo,x,y,z,wfnval,wfnderv,wfnhess,wfntens) 1143grad=0 1144do imo=1,nmo 1145 wfnlapl=wfnhess(1,1,imo)+wfnhess(2,2,imo)+wfnhess(3,3,imo) 1146 gx=wfnderv(1,imo)*wfnlapl+wfnval(imo)*(wfntens(1,1,1,imo)+wfntens(1,2,2,imo)+wfntens(1,3,3,imo)) 1147 gy=wfnderv(2,imo)*wfnlapl+wfnval(imo)*(wfntens(1,1,2,imo)+wfntens(2,2,2,imo)+wfntens(2,3,3,imo)) 1148 gz=wfnderv(3,imo)*wfnlapl+wfnval(imo)*(wfntens(1,1,3,imo)+wfntens(2,2,3,imo)+wfntens(3,3,3,imo)) 1149 grad(1)=grad(1)+gx*MOocc(imo) 1150 grad(2)=grad(2)+gy*MOocc(imo) 1151 grad(3)=grad(3)+gz*MOocc(imo) 1152end do 1153grad=grad/2 1154end subroutine 1155!!------- Output gradient norm of energy density (i.e. negative of K(r)) 1156real*8 function energydens_grdn(x,y,z) 1157real*8 x,y,z,grad(3) 1158call energydens_grad(x,y,z,grad) 1159energydens_grdn=dsqrt(sum(grad**2)) 1160end function 1161!!------- Output Laplacian of energy density (i.e. negative of K(r)) 1162real*8 function energydens_lapl(x,y,z) 1163real*8 x,y,z,gradaddx(3),gradminx(3),gradaddy(3),gradminy(3),gradaddz(3),gradminz(3) 1164diff=1D-5 1165denom=2*diff 1166call energydens_grad(x+diff,y,z,gradaddx) 1167call energydens_grad(x-diff,y,z,gradminx) 1168call energydens_grad(x,y+diff,z,gradaddy) 1169call energydens_grad(x,y-diff,z,gradminy) 1170call energydens_grad(x,y,z+diff,gradaddz) 1171call energydens_grad(x,y,z-diff,gradminz) 1172xlapl=(gradaddx(1)-gradminx(1))/denom 1173ylapl=(gradaddy(2)-gradminy(2))/denom 1174zlapl=(gradaddz(3)-gradminz(3))/denom 1175energydens_lapl=xlapl+ylapl+zlapl 1176end function 1177 1178 1179 1180!!!--------- Output electron linear momentum density in 3D representation at a point. idir=0/1/2/3 means magnitude/x/y/z component 1181real*8 function elemomdens(x,y,z,idir) 1182real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo),comp(0:3) 1183integer idir 1184call orbderv(2,1,nmo,x,y,z,wfnval,wfnderv) 1185comp=0 1186do imo=1,nmo 1187 comp(1:3)=comp(1:3)-MOocc(imo)*wfnval(imo)*wfnderv(1:3,imo) 1188end do 1189comp(0)=sum(comp(1:3)**2) 1190elemomdens=comp(idir) 1191end function 1192 1193 1194!!!--------- Output magnetic dipole moment density at a point. idir=0/1/2/3 means magnitude/x/y/z component 1195real*8 function magmomdens(x,y,z,idir) 1196real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo),comp(0:3) 1197integer idir 1198call orbderv(2,1,nmo,x,y,z,wfnval,wfnderv) 1199comp=0 1200do imo=1,nmo 1201 tmpx=-wfnval(imo)*(y*wfnderv(3,imo)-z*wfnderv(2,imo)) 1202 tmpy=-wfnval(imo)*(z*wfnderv(1,imo)-x*wfnderv(3,imo)) 1203 tmpz=-wfnval(imo)*(x*wfnderv(2,imo)-y*wfnderv(1,imo)) 1204 comp(1)=comp(1)+MOocc(imo)*tmpx 1205 comp(2)=comp(2)+MOocc(imo)*tmpy 1206 comp(3)=comp(3)+MOocc(imo)*tmpz 1207end do 1208comp(0)=sum(comp(1:3)**2) 1209magmomdens=comp(idir) 1210end function 1211 1212 1213!!!----- Calculate Sign(lambda2(r))*rho(r) function, this is a warpper used to convert subroutine to function form 1214real*8 function signlambda2rho(x,y,z) 1215real*8 x,y,z,sl2r,RDG,rho 1216call signlambda2rho_RDG(x,y,z,sl2r,RDG,rho) 1217signlambda2rho=sl2r 1218end function 1219!!!------ Calculate signlambda2rho and RDG at the same time 1220subroutine signlambda2rho_RDG(x,y,z,sl2r,RDG,elerho) 1221use util 1222real*8 x,y,z,elerho,sl2r,RDG 1223real*8 eigvecmat(3,3),eigval(3),elehess(3,3),elegrad(3) !Hessian of electron density 1224call calchessmat_dens(2,x,y,z,elerho,elegrad,elehess) 1225call diagmat(elehess,eigvecmat,eigval,100,1D-10) 1226call sort(eigval) 1227if (eigval(2)==0D0) then !When eigval(2)==0.0D0, eigval(2)/abs(eigval(2)) can't be calculated, elerho generally will be zero, so sign is not important 1228 sl2r=elerho 1229else 1230 sl2r=elerho*eigval(2)/abs(eigval(2)) 1231end if 1232sumgrad2=sum(elegrad(:)**2) 1233if (RDG_maxrho/=0D0.and.elerho>=RDG_maxrho) then 1234 RDG=100D0 1235!This occurs at distant region when exponent cutoff is used, the actual value should be very large. In order to avoid denominator become zero, we set it artifically to a big value 1236else if (sumgrad2==0D0.or.elerho==0D0) then 1237 RDG=999D0 1238else 1239 RDG=0.161620459673995D0*dsqrt(sumgrad2)/elerho**(4D0/3D0) !0.161620459673995D0=1/(2*(3*pi**2)**(1/3)) 1240end if 1241end subroutine 1242 1243 1244 1245!!!-----Output ELF(Electron Localization Function) or LOL(Localized orbital locator) and similar function value at a point 1246!label="ELF" or "LOL" 1247real*8 function ELF_LOL(x,y,z,label) 1248real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo) 1249real*8 D,Dh,gradrho(3),gradrhoa(3),gradrhob(3),rho,rhoa,rhob,rhospin,MOoccnow 1250real*8 :: Fc=2.871234000D0 ! Fermi constant = (3/10)*(3*Pi^2)**(2/3) = 2.871234, 1/2.871234=0.34828 1251real*8 :: Fc_pol=4.557799872D0 ! Fermi constant for spin polarized = (3/10)*(6*Pi^2)**(2/3) = 4.5578, 1/4.5578=0.2194 1252character*3 label 1253 1254!The ELF defined by Tsirelson, CPL,351,142. 1255!The LOL defined by Tsirelson, Acta Cryst. (2002). B58, 780-785. 1256!These functions are not important, so I don't write a special code 1257!for it, since rho, nebla-rho, nebla^2-rho support EDF, this function also support EDF 1258if (ELFLOL_type==1) then 1259 rho=fdens(x,y,z) 1260 if (wfntype==0.or.wfntype==3) then !close shell cases 1261 Dh=Fc*rho**(5.0D0/3.0D0) 1262 else if (wfntype==1.or.wfntype==2.or.wfntype==4) then !Open shell cases 1263 rhospin=fspindens(x,y,z,'s') !rhospin=rhoa-rhob, rho=rhoa+rhob 1264 rhoa=(rhospin+rho)/2D0 1265 rhob=(rho-rhospin)/2D0 1266 Dh=Fc_pol*(rhoa**(5.0D0/3.0D0)+rhob**(5.0D0/3.0D0)) !kinetic energy of HEG 1267 end if 1268 if (label=="ELF") then 1269 !Restrictly speaking, the kinetic energy expansion should be replace by polarized form for open-shell 1270 D=Dh-(1/9D0)*fgrad(x,y,z,'t')**2/rho+(1/6D0)*flapl(x,y,z,'t') 1271 ELF_LOL=1/(1+(D/Dh)**2) 1272 else if (label=="LOL") then 1273 D=Dh+(1/72D0)*fgrad(x,y,z,'t')**2/rho+(1/6D0)*flapl(x,y,z,'t') 1274 t=Dh/D 1275 ELF_LOL=1D0/(1D0/t+1) 1276 end if 1277 if (ELF_LOL<ELFLOL_cut) ELF_LOL=0 1278 return 1279end if 1280 1281!For Becke's and Lu Tian's definition below 1282call orbderv(2,1,nmo,x,y,z,wfnval,wfnderv) 1283 1284D=0.0D0 1285rho=0.0D0 1286rhoa=0D0 1287rhob=0D0 1288gradrho=0D0 1289gradrhoa=0D0 1290gradrhob=0D0 1291if (label=="ELF") then 1292 if (wfntype==0.or.wfntype==3) then !spin-unpolarized case 1293 do i=1,nmo 1294 rho=rho+MOocc(i)*wfnval(i)**2 1295 gradrho(:)=gradrho(:)+2.0D0*MOocc(i)*wfnval(i)*wfnderv(:,i) 1296 D=D+MOocc(i)*(sum(wfnderv(:,i)**2)) !Calculate actual kinetic term 1297 end do 1298 Dh=Fc*rho**(5.0D0/3.0D0) !Thomas-Fermi uniform electron gas kinetic energy 1299 D=D/2D0 1300 if (rho/=0D0) D=D-(sum(gradrho(:)**2))/rho/8D0 1301 else if (wfntype==1.or.wfntype==2.or.wfntype==4) then !spin-polarized case 1302 do i=1,nmo 1303 MOoccnow=MOocc(i) 1304 if (MOtype(i)==0) MOoccnow=MOocc(i)/2D0 !double occupied, present when wfntype==2 (ROHF), alpha and beta get half part 1305 if (MOtype(i)==1.or.MOtype(i)==0) then 1306 rhoa=rhoa+MOoccnow*wfnval(i)**2 1307 gradrhoa(:)=gradrhoa(:)+2.0D0*MOoccnow*wfnval(i)*wfnderv(:,i) 1308 end if 1309 if (MOtype(i)==2.or.MOtype(i)==0) then 1310 rhob=rhob+MOoccnow*wfnval(i)**2 1311 gradrhob(:)=gradrhob(:)+2.0D0*MOoccnow*wfnval(i)*wfnderv(:,i) 1312 end if 1313 D=D+MOocc(i)*(sum(wfnderv(:,i)**2)) !Calculate actual kinetic term 1314 end do 1315 Dh=Fc_pol*(rhoa**(5.0D0/3.0D0)+rhob**(5.0D0/3.0D0)) 1316 D=D/2D0 1317 if (rhoa/=0D0) D=D-(sum(gradrhoa(:)**2))/rhoa/8 1318 if (rhob/=0D0) D=D-(sum(gradrhob(:)**2))/rhob/8 1319 end if 1320 if (ELFLOL_type==0.or.ELFLOL_type==2) then 1321 if (ELF_addminimal==1) D=D+1D-5 !add 1D-5 to avoid D become zero, leads to unity in infinite 1322 if (ELFLOL_type==0) ELF_LOL=1/(1+(D/Dh)**2) 1323 if (ELFLOL_type==2) ELF_LOL=1/(1+(D/Dh)) !New formalism defined by Tian Lu 1324 end if 1325 if (ELFLOL_type==3) ELF_LOL=D/Dh 1326 if (ELFLOL_type==995) ELF_LOL=Dh !Thomas-Fermi kinetic energy density 1327 if (ELFLOL_type==996) ELF_LOL=D/Dh !For test 1328 if (ELFLOL_type==997) ELF_LOL=D !For test 1329 if (ELFLOL_type==998) ELF_LOL=1/(1+D) !For test 1330 if (ELFLOL_type==999) ELF_LOL=1/(1+D**2) !For test 1331!---------------------------- 1332else if (label=="LOL") then 1333 t=0.0D0 1334 if (wfntype==0.or.wfntype==3) then !Spin-unpolarized case 1335 do i=1,nmo !Store actual kinetic energy to t first 1336 rho=rho+MOocc(i)*wfnval(i)**2 1337 t=t+MOocc(i)*(sum(wfnderv(:,i)**2)) 1338 end do 1339 t=t/2D0 1340 Dh=Fc*rho**(5.0D0/3.0D0) 1341 else if (wfntype==1.or.wfntype==2.or.wfntype==4) then !spin-polarized case 1342 do i=1,nmo !Store actual kinetic energy to t first 1343 MOoccnow=MOocc(i) 1344 if (MOtype(i)==0) MOoccnow=MOocc(i)/2D0 !double occupied, present when wfntype==2 (ROHF), alpha and beta get half part 1345 if (MOtype(i)==1.or.MOtype(i)==0) rhoa=rhoa+MOoccnow*wfnval(i)**2 1346 if (MOtype(i)==2.or.MOtype(i)==0) rhob=rhob+MOoccnow*wfnval(i)**2 1347 t=t+MOocc(i)*(sum(wfnderv(:,i)**2)) 1348 end do 1349 t=t/2D0 1350 Dh=Fc_pol*(rhoa**(5.0D0/3.0D0)+rhob**(5.0D0/3.0D0)) 1351 end if 1352!--------- A new definition of LOL, however the value range is not as good as LOL 1353! ELF_LOL=t-Dh 1354! if (ELF_LOL>0) then 1355! ELF_LOL=1D0/(1D0+1D0/ELF_LOL) 1356! else if (ELF_LOL<0) then 1357! ELF_LOL=-1D0/(1D0+1D0/abs(ELF_LOL)) 1358! end if 1359!------------- 1360!If there is very long distance between molecule and current point, t (above) is zero, 1361!and t=Dh/t is also zero (because rho converges faster), but can't be calculate directly, so simply skip 1362 if (t/=0.0D0) t=Dh/t 1363 if (ELFLOL_type==0) ELF_LOL=1D0/(1D0/t+1) !namely t/(1+t). This is default case 1364 if (ELFLOL_type==2) ELF_LOL=1D0/((1D0/t)**2+1) !New form defined by Tian Lu 1365end if 1366if (ELF_LOL<ELFLOL_cut) ELF_LOL=0 1367end function 1368 1369 1370!!!----- Calculate ELF/LOL, its gradient and Hessian matrix at x,y,z, store to hess 1371!!!!!!!!!!!! currently can not calculate Hessian 1372!funsel="ELF" or "LOL" 1373!itype=1 Only calculate value and gradient, not Hessian 1374!itype=2 Calculate value, gradient and Hessian (Not available) 1375subroutine calchessmat_ELF_LOL(itype,x,y,z,value,grad,hess,funsel) 1376use util 1377integer itype 1378real*8 x,y,z,value,grad(3),hess(3,3),MOoccnow 1379real*8 wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo),wfntens3(3,3,3,nmo) 1380real*8 rho,gradrho(3),hessrho(3,3),Dh,gradDh(3),hessDh(3,3),Ts,gradTs(3),hessTs(3,3),Wei,gradWei(3),hessWei(3,3) 1381real*8 rhoa,rhob,gradrhoa(3),gradrhob(3),hessrhoa(3,3),hessrhob(3,3),Dha,Dhb,gradDha(3),gradDhb(3),Weia,Weib,gradWeia(3),gradWeib(3) 1382! real*8 hessDha(3,3),hessDhb(3,3),hessWeia(3,3),hessWeib(3,3) 1383real*8 :: Fc=2.871234000D0 ! Fermi constant = (3/10)*(3*Pi^2)**(2/3) = 2.871234, 1/2.871234=0.34828 1384real*8 :: Fc_pol=4.557799872D0 ! Fermi constant for spin polarized = (3/10)*(6*Pi^2)**(2/3) = 4.5578, 1/4.5578=0.2194 1385real*8 :: corrELF=1D-5 1386character funsel*3 1387 1388if (itype==1) call orbderv(4,1,nmo,x,y,z,wfnval,wfnderv,wfnhess) !Get Hessian of GTF, needn't 3-order tensor 1389if (itype==2) call orbderv(5,1,nmo,x,y,z,wfnval,wfnderv,wfnhess,wfntens3) 1390 1391!spin-unpolarized case 1392if (wfntype==0.or.wfntype==3) then 1393 rho=0D0 1394 gradrho=0D0 1395 Ts=0D0 1396 gradTs=0D0 1397 do i=1,nmo 1398 rho=rho+MOocc(i)*wfnval(i)**2 1399 gradrho(:)=gradrho(:)+MOocc(i)*wfnval(i)*wfnderv(:,i) 1400 Ts=Ts+MOocc(i)*(sum(wfnderv(:,i)**2)) 1401 end do 1402 gradrho=2*gradrho 1403 Ts=Ts/2D0 1404 Dh=Fc*rho**(5D0/3D0) 1405 gradDh(:)=5D0/3D0*Fc*rho**(2D0/3D0)*gradrho(:) 1406 do i=1,nmo 1407 gradTs(1)=gradTs(1)+MOocc(i)*(wfnderv(1,i)*wfnhess(1,1,i)+wfnderv(2,i)*wfnhess(1,2,i)+wfnderv(3,i)*wfnhess(1,3,i)) 1408 gradTs(2)=gradTs(2)+MOocc(i)*(wfnderv(2,i)*wfnhess(2,2,i)+wfnderv(1,i)*wfnhess(2,1,i)+wfnderv(3,i)*wfnhess(2,3,i)) 1409 gradTs(3)=gradTs(3)+MOocc(i)*(wfnderv(3,i)*wfnhess(3,3,i)+wfnderv(1,i)*wfnhess(3,1,i)+wfnderv(2,i)*wfnhess(3,2,i)) 1410 end do 1411 if (funsel=="ELF") then 1412 !Calculate Hessian for rho 1413 hessrho(1,1)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)**2 + wfnval(1:nmo)*wfnhess(1,1,1:nmo) ) ) 1414 hessrho(2,2)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)**2 + wfnval(1:nmo)*wfnhess(2,2,1:nmo) ) ) 1415 hessrho(3,3)=2*sum( MOocc(1:nmo)*( wfnderv(3,1:nmo)**2 + wfnval(1:nmo)*wfnhess(3,3,1:nmo) ) ) 1416 hessrho(1,2)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(2,1:nmo)+wfnhess(1,2,1:nmo)*wfnval(1:nmo) ) ) 1417 hessrho(2,3)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)*wfnderv(3,1:nmo)+wfnhess(2,3,1:nmo)*wfnval(1:nmo) ) ) 1418 hessrho(1,3)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(3,1:nmo)+wfnhess(1,3,1:nmo)*wfnval(1:nmo) ) ) 1419 hessrho(2,1)=hessrho(1,2) 1420 hessrho(3,1)=hessrho(1,3) 1421 hessrho(3,2)=hessrho(2,3) 1422 !Calculate Weizsacker functional and its derivatives 1423 Wei=sum(gradrho(:)**2)/8D0/rho 1424 D=Ts-Wei+corrELF 1425 chi=D/Dh 1426 value=1D0/(1D0+chi**2) 1427 gradWei(1)= 0.25D0/rho*( gradrho(1)*hessrho(1,1)+gradrho(2)*hessrho(1,2)+gradrho(3)*hessrho(1,3) ) - wei/rho*gradrho(1) 1428 gradWei(2)= 0.25D0/rho*( gradrho(2)*hessrho(2,2)+gradrho(1)*hessrho(2,1)+gradrho(3)*hessrho(2,3) ) - wei/rho*gradrho(2) 1429 gradWei(3)= 0.25D0/rho*( gradrho(3)*hessrho(3,3)+gradrho(2)*hessrho(3,2)+gradrho(1)*hessrho(3,1) ) - wei/rho*gradrho(3) 1430 chidx=(gradTs(1)-gradWei(1))/Dh - gradDh(1)/Dh**2 *D 1431 chidy=(gradTs(2)-gradWei(2))/Dh - gradDh(2)/Dh**2 *D 1432 chidz=(gradTs(3)-gradWei(3))/Dh - gradDh(3)/Dh**2 *D 1433 grad(1)=-2D0*chi/(1+chi**2)**2 * chidx 1434 grad(2)=-2D0*chi/(1+chi**2)**2 * chidy 1435 grad(3)=-2D0*chi/(1+chi**2)**2 * chidz 1436 else if (funsel=="LOL") then 1437 value=1D0/(1D0+Ts/Dh) 1438 tmp=-1D0/Dh/(1D0+Ts/Dh)**2 1439 grad(1)=tmp*(gradTs(1)-Ts/Dh*gradDh(1)) 1440 grad(2)=tmp*(gradTs(2)-Ts/Dh*gradDh(2)) 1441 grad(3)=tmp*(gradTs(3)-Ts/Dh*gradDh(3)) 1442 end if 1443!spin-polarized case 1444else if (wfntype==1.or.wfntype==2.or.wfntype==4) then 1445 rhoa=0D0 1446 rhob=0D0 1447 gradrhoa=0D0 1448 gradrhob=0D0 1449 Ts=0D0 1450 gradTs=0D0 1451 do i=1,nmo 1452 MOoccnow=MOocc(i) 1453 if (MOtype(i)==0) MOoccnow=MOocc(i)/2D0 1454 if (MOtype(i)==1.or.MOtype(i)==0) then 1455 rhoa=rhoa+MOoccnow*wfnval(i)**2 1456 gradrhoa(:)=gradrhoa(:)+MOoccnow*wfnval(i)*wfnderv(:,i) 1457 end if 1458 if (MOtype(i)==2.or.MOtype(i)==0) then 1459 rhob=rhob+MOoccnow*wfnval(i)**2 1460 gradrhob(:)=gradrhob(:)+MOoccnow*wfnval(i)*wfnderv(:,i) 1461 end if 1462 Ts=Ts+MOocc(i)*(sum(wfnderv(:,i)**2)) 1463 end do 1464 gradrhoa=2*gradrhoa 1465 gradrhob=2*gradrhob 1466 Ts=Ts/2D0 1467 Dha=Fc_pol*rhoa**(5D0/3D0) 1468 Dhb=Fc_pol*rhob**(5D0/3D0) 1469 Dh=Dha+Dhb 1470 gradDha(:)=5D0/3D0*Fc_pol*rhoa**(2D0/3D0)*gradrhoa(:) 1471 gradDhb(:)=5D0/3D0*Fc_pol*rhob**(2D0/3D0)*gradrhob(:) 1472 gradDh=gradDha+gradDhb 1473 do i=1,nmo 1474 gradTs(1)=gradTs(1)+MOocc(i)*(wfnderv(1,i)*wfnhess(1,1,i)+wfnderv(2,i)*wfnhess(1,2,i)+wfnderv(3,i)*wfnhess(1,3,i)) 1475 gradTs(2)=gradTs(2)+MOocc(i)*(wfnderv(2,i)*wfnhess(2,2,i)+wfnderv(1,i)*wfnhess(2,1,i)+wfnderv(3,i)*wfnhess(2,3,i)) 1476 gradTs(3)=gradTs(3)+MOocc(i)*(wfnderv(3,i)*wfnhess(3,3,i)+wfnderv(1,i)*wfnhess(3,1,i)+wfnderv(2,i)*wfnhess(3,2,i)) 1477 end do 1478 1479 if (funsel=="ELF") then 1480! !Calculate Hessian for rho 1481 hessrhoa=0D0 1482 hessrhob=0D0 1483 do i=1,nmo 1484 MOoccnow=MOocc(i) 1485 if (MOtype(i)==0) MOoccnow=MOocc(i)/2D0 !double occupied, present when wfntype==2 (ROHF), alpha and beta get half part 1486 if (MOtype(i)==1.or.MOtype(i)==0) then 1487 hessrhoa(1,1)=hessrhoa(1,1)+MOoccnow*( wfnderv(1,i)**2 + wfnval(i)*wfnhess(1,1,i) ) 1488 hessrhoa(2,2)=hessrhoa(2,2)+MOoccnow*( wfnderv(2,i)**2 + wfnval(i)*wfnhess(2,2,i) ) 1489 hessrhoa(3,3)=hessrhoa(3,3)+MOoccnow*( wfnderv(3,i)**2 + wfnval(i)*wfnhess(3,3,i) ) 1490 hessrhoa(1,2)=hessrhoa(1,2)+MOoccnow*( wfnderv(1,i)*wfnderv(2,i)+wfnhess(1,2,i)*wfnval(i) ) 1491 hessrhoa(2,3)=hessrhoa(2,3)+MOoccnow*( wfnderv(2,i)*wfnderv(3,i)+wfnhess(2,3,i)*wfnval(i) ) 1492 hessrhoa(1,3)=hessrhoa(1,3)+MOoccnow*( wfnderv(1,i)*wfnderv(3,i)+wfnhess(1,3,i)*wfnval(i) ) 1493 end if 1494 if (MOtype(i)==2.or.MOtype(i)==0) then 1495 hessrhob(1,1)=hessrhob(1,1)+MOoccnow*( wfnderv(1,i)**2 + wfnval(i)*wfnhess(1,1,i) ) 1496 hessrhob(2,2)=hessrhob(2,2)+MOoccnow*( wfnderv(2,i)**2 + wfnval(i)*wfnhess(2,2,i) ) 1497 hessrhob(3,3)=hessrhob(3,3)+MOoccnow*( wfnderv(3,i)**2 + wfnval(i)*wfnhess(3,3,i) ) 1498 hessrhob(1,2)=hessrhob(1,2)+MOoccnow*( wfnderv(1,i)*wfnderv(2,i)+wfnhess(1,2,i)*wfnval(i) ) 1499 hessrhob(2,3)=hessrhob(2,3)+MOoccnow*( wfnderv(2,i)*wfnderv(3,i)+wfnhess(2,3,i)*wfnval(i) ) 1500 hessrhob(1,3)=hessrhob(1,3)+MOoccnow*( wfnderv(1,i)*wfnderv(3,i)+wfnhess(1,3,i)*wfnval(i) ) 1501 end if 1502 end do 1503 hessrhoa=hessrhoa*2 1504 hessrhob=hessrhob*2 1505 hessrhoa(2,1)=hessrhoa(1,2) 1506 hessrhoa(3,1)=hessrhoa(1,3) 1507 hessrhoa(3,2)=hessrhoa(2,3) 1508 hessrhob(2,1)=hessrhob(1,2) 1509 hessrhob(3,1)=hessrhob(1,3) 1510 hessrhob(3,2)=hessrhob(2,3) 1511! !Calculate Weizsacker functional and its derivatives 1512 Weia=sum(gradrhoa(:)**2)/8D0/rhoa 1513 Weib=sum(gradrhob(:)**2)/8D0/rhob 1514 Wei=Weia+Weib 1515 D=Ts-Wei+corrELF 1516 chi=D/Dh 1517 value=1D0/(1D0+chi**2) 1518 gradWeia(1)= 0.25D0/rhoa*( gradrhoa(1)*hessrhoa(1,1)+gradrhoa(2)*hessrhoa(1,2)+gradrhoa(3)*hessrhoa(1,3) ) - weia/rhoa*gradrhoa(1) 1519 gradWeia(2)= 0.25D0/rhoa*( gradrhoa(2)*hessrhoa(2,2)+gradrhoa(1)*hessrhoa(2,1)+gradrhoa(3)*hessrhoa(2,3) ) - weia/rhoa*gradrhoa(2) 1520 gradWeia(3)= 0.25D0/rhoa*( gradrhoa(3)*hessrhoa(3,3)+gradrhoa(2)*hessrhoa(3,2)+gradrhoa(1)*hessrhoa(3,1) ) - weia/rhoa*gradrhoa(3) 1521 gradWeib(1)= 0.25D0/rhob*( gradrhob(1)*hessrhob(1,1)+gradrhob(2)*hessrhob(1,2)+gradrhob(3)*hessrhob(1,3) ) - weib/rhob*gradrhob(1) 1522 gradWeib(2)= 0.25D0/rhob*( gradrhob(2)*hessrhob(2,2)+gradrhob(1)*hessrhob(2,1)+gradrhob(3)*hessrhob(2,3) ) - weib/rhob*gradrhob(2) 1523 gradWeib(3)= 0.25D0/rhob*( gradrhob(3)*hessrhob(3,3)+gradrhob(2)*hessrhob(3,2)+gradrhob(1)*hessrhob(3,1) ) - weib/rhob*gradrhob(3) 1524 gradWei=gradWeia+gradWeib 1525 chidx=(gradTs(1)-gradWei(1))/Dh - gradDh(1)/Dh**2 *D 1526 chidy=(gradTs(2)-gradWei(2))/Dh - gradDh(2)/Dh**2 *D 1527 chidz=(gradTs(3)-gradWei(3))/Dh - gradDh(3)/Dh**2 *D 1528 grad(1)=-2D0*chi/(1+chi**2)**2 * chidx 1529 grad(2)=-2D0*chi/(1+chi**2)**2 * chidy 1530 grad(3)=-2D0*chi/(1+chi**2)**2 * chidz 1531 else if (funsel=="LOL") then 1532 value=1D0/(1D0+Ts/Dh) 1533 tmp=-1D0/Dh/(1D0+Ts/Dh)**2 1534 grad(1)=tmp*(gradTs(1)-Ts/Dh*gradDh(1)) 1535 grad(2)=tmp*(gradTs(2)-Ts/Dh*gradDh(2)) 1536 grad(3)=tmp*(gradTs(3)-Ts/Dh*gradDh(3)) 1537 end if 1538end if 1539 1540! if (itype==1) return 1541! ! Calculate Hessian for LOL, also need Hessian for rho 1542! hessrho(1,1)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)**2 + wfnval(1:nmo)*wfnhess(1,1,1:nmo) ) ) 1543! hessrho(2,2)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)**2 + wfnval(1:nmo)*wfnhess(2,2,1:nmo) ) ) 1544! hessrho(3,3)=2*sum( MOocc(1:nmo)*( wfnderv(3,1:nmo)**2 + wfnval(1:nmo)*wfnhess(3,3,1:nmo) ) ) 1545! hessrho(1,2)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(2,1:nmo)+wfnhess(1,2,1:nmo)*wfnval(1:nmo) ) ) 1546! hessrho(2,3)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)*wfnderv(3,1:nmo)+wfnhess(2,3,1:nmo)*wfnval(1:nmo) ) ) 1547! hessrho(1,3)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(3,1:nmo)+wfnhess(1,3,1:nmo)*wfnval(1:nmo) ) ) 1548! hessrho(2,1)=hessrho(1,2) 1549! hessrho(3,1)=hessrho(1,3) 1550! hessrho(3,2)=hessrho(2,3) 1551! 1552! hessTs=0D0 1553! do i=1,nmo 1554! hessTs(1,1)=hessTs(1,1)+MOocc(i)*( wfnhess(1,1,i)**2 + wfnderv(1,i)*wfntens3(1,1,1,i) + wfnhess(1,2,i)**2 + wfnderv(2,i)*wfntens3(1,1,2,i) + wfnhess(1,3,i)**2 + wfnderv(3,i)*wfntens3(1,1,3,i) ) 1555! hessTs(2,2)=hessTs(2,2)+MOocc(i)*( wfnhess(2,2,i)**2 + wfnderv(2,i)*wfntens3(2,2,2,i) + wfnhess(2,1,i)**2 + wfnderv(1,i)*wfntens3(2,2,1,i) + wfnhess(2,3,i)**2 + wfnderv(3,i)*wfntens3(2,2,3,i) ) 1556! hessTs(3,3)=hessTs(3,3)+MOocc(i)*( wfnhess(3,3,i)**2 + wfnderv(3,i)*wfntens3(3,3,3,i) + wfnhess(3,2,i)**2 + wfnderv(2,i)*wfntens3(3,3,2,i) + wfnhess(3,1,i)**2 + wfnderv(1,i)*wfntens3(3,3,1,i) ) 1557! hessTs(1,2)=hessTs(1,2)+MOocc(i)*( wfnhess(1,2,i)*wfnhess(1,1,i) + wfnderv(1,i)*wfntens3(1,1,2,i) + wfnhess(2,2,i)*wfnhess(1,2,i) + wfnderv(2,i)*wfntens3(1,2,2,i) + wfnhess(2,3,i)*wfnhess(1,3,i) + wfnderv(3,i)*wfntens3(1,2,3,i) ) 1558! hessTs(2,3)=hessTs(2,3)+MOocc(i)*( wfnhess(2,3,i)*wfnhess(2,2,i) + wfnderv(2,i)*wfntens3(2,2,3,i) + wfnhess(3,3,i)*wfnhess(2,3,i) + wfnderv(3,i)*wfntens3(2,3,3,i) + wfnhess(3,1,i)*wfnhess(2,1,i) + wfnderv(1,i)*wfntens3(2,3,1,i) ) 1559! hessTs(1,3)=hessTs(1,3)+MOocc(i)*( wfnhess(1,3,i)*wfnhess(1,1,i) + wfnderv(1,i)*wfntens3(1,1,3,i) + wfnhess(3,3,i)*wfnhess(1,3,i) + wfnderv(3,i)*wfntens3(1,3,3,i) + wfnhess(3,2,i)*wfnhess(1,2,i) + wfnderv(2,i)*wfntens3(1,3,2,i) ) 1560! end do 1561! hessTs(2,1)=hessTs(1,2) 1562! hessTs(3,1)=hessTs(1,3) 1563! hessTs(3,2)=hessTs(2,3) 1564! 1565! tmp1=10D0/9D0*Fc/rho**(1D0/3D0) 1566! tmp2=5D0/3D0*Fc*rho**(2D0/3D0) 1567! hessDh(1,1)=tmp1*gradrho(1)**2 + tmp2*hessrho(1,1) 1568! hessDh(2,2)=tmp1*gradrho(2)**2 + tmp2*hessrho(2,2) 1569! hessDh(3,3)=tmp1*gradrho(3)**2 + tmp2*hessrho(3,3) 1570! hessDh(1,2)=tmp1*gradrho(1)*gradrho(2) + tmp2*hessrho(1,2) 1571! hessDh(2,3)=tmp1*gradrho(2)*gradrho(3) + tmp2*hessrho(2,3) 1572! hessDh(1,3)=tmp1*gradrho(1)*gradrho(3) + tmp2*hessrho(1,3) 1573! hessDh(2,1)=hessDh(1,2) 1574! hessDh(3,1)=hessDh(1,3) 1575! hessDh(3,2)=hessDh(2,3) 1576! 1577! !Diagonal of Hessian of LOL 1578! apre=1/Dh**2/(1+Ts/Dh)**3 1579! bpre=-1/Dh/(1+Ts/Dh)**2 1580! apartxx=(gradDh(1)+2*gradTs(1)-Ts/Dh*gradDh(1)) * (gradTs(1)-Ts/Dh*gradDh(1)) 1581! bpartxx=hessTs(1,1)-( gradDh(1)*gradTs(1)-Ts/Dh*gradDh(1)**2+Ts*hessDh(1,1) )/Dh 1582! hess(1,1)=apre*apartxx+bpre*bpartxx 1583! apartyy=(gradDh(2)+2*gradTs(2)-Ts/Dh*gradDh(2)) * (gradTs(2)-Ts/Dh*gradDh(2)) 1584! bpartyy=hessTs(2,2)-( gradDh(2)*gradTs(2)-Ts/Dh*gradDh(2)**2+Ts*hessDh(2,2) )/Dh 1585! hess(2,2)=apre*apartyy+bpre*bpartyy 1586! apartzz=(gradDh(3)+2*gradTs(3)-Ts/Dh*gradDh(3)) * (gradTs(3)-Ts/Dh*gradDh(3)) 1587! bpartzz=hessTs(3,3)-( gradDh(3)*gradTs(3)-Ts/Dh*gradDh(3)**2+Ts*hessDh(3,3) )/Dh 1588! hess(3,3)=apre*apartzz+bpre*bpartzz 1589! !Non-diagonal of Hessian of LOL 1590! bpartxy=hessTs(1,2)-( gradDh(1)*gradTs(2)-Ts/Dh*gradDh(1)*gradDh(2)+Ts*hessDh(1,2) )/Dh 1591! hess(1,2)=apre*apartxx+bpre*bpartxy 1592! bpartyz=hessTs(2,3)-( gradDh(2)*gradTs(3)-Ts/Dh*gradDh(2)*gradDh(3)+Ts*hessDh(2,3) )/Dh 1593! hess(2,3)=apre*apartyy+bpre*bpartyz 1594! bpartxz=hessTs(1,3)-( gradDh(1)*gradTs(3)-Ts/Dh*gradDh(1)*gradDh(3)+Ts*hessDh(1,3) )/Dh 1595! hess(1,3)=apre*apartxx+bpre*bpartxz 1596! hess(2,1)=hess(1,2) 1597! hess(3,1)=hess(1,3) 1598! hess(3,2)=hess(2,3) 1599end subroutine 1600 1601 1602!-------- Output average local ionization energy at a point 1603real*8 function avglocion(x,y,z) 1604real*8 x,y,z,wfnval(nmo) 1605call orbderv(1,1,nmo,x,y,z,wfnval) 1606avglocion=0D0 1607rho=0D0 1608do i=1,nmo 1609 avglocion=avglocion+abs(MOene(i))*MOocc(i)*wfnval(i)**2 1610 rho=rho+MOocc(i)*wfnval(i)**2 !Calculate rho 1611end do 1612if (rho==0D0) then 1613 avglocion=0D0 !Avoid at distant region rho becomes zero when exponent cutoff is used 1614else 1615 avglocion=avglocion/rho 1616end if 1617end function 1618 1619!-------- Calculate average local ionization energy at a point and meantime decompose it to occupied orbitals contribution 1620subroutine avglociondecomp(ifileid,x,y,z) 1621real*8 x,y,z,wfnval(nmo) 1622integer ifileid 1623character orbtype*2 1624call orbderv(1,1,nmo,x,y,z,wfnval) 1625totALIE=0D0 1626rho=0D0 1627do i=1,nmo 1628 totALIE=totALIE+abs(MOene(i))*MOocc(i)*wfnval(i)**2 1629 rho=rho+MOocc(i)*wfnval(i)**2 !Calculate rho 1630end do 1631if (rho==0D0) then 1632 totALIE=0D0 !Avoid at distant region rho becomes zero when exponent cutoff is used 1633else 1634 totALIE=totALIE/rho 1635end if 1636write(ifileid,"(' Average local ionization energy:',f16.10,' a.u.')") totALIE 1637write(ifileid,*) "Contribution of each orbital to average local ionization energy (a.u.):" 1638do i=1,nmo 1639 if (MOtype(i)==0) orbtype="AB" 1640 if (MOtype(i)==1) orbtype="A " 1641 if (MOtype(i)==2) orbtype="B " 1642 write(ifileid,"(' Orbital',i6,' Ene:',f12.6,' Occ:',f5.2,' Type:',a,' Contribution:',f12.6)") i,MOene(i),MOocc(i),orbtype,abs(MOene(i))*MOocc(i)*wfnval(i)**2/rho 1643end do 1644end subroutine 1645 1646 1647!-------- Output local electron affinity at a point 1648!Since virtual orbitals are involved, such as .fch/.molden/.gms must be used 1649real*8 function loceleaff(x,y,z) 1650real*8 x,y,z,wfnval(nmo) 1651call orbderv(1,1,nmo,x,y,z,wfnval) 1652loceleaff=0D0 1653rho=0D0 1654do i=1,nmo 1655 if (MOocc(i)==0) then !Only cycles unoccupied orbitals 1656 loceleaff=loceleaff-MOene(i)*wfnval(i)**2 !Don't need to multiply "occupation number", because ROHF is not allowed, so all orbitals have the same type 1657 rho=rho+wfnval(i)**2 !Calculate rho 1658 end if 1659end do 1660if (rho==0D0) then 1661 loceleaff=0D0 !Avoid at distant region rho become zero when exponent cutoff is used 1662else 1663 loceleaff=loceleaff/rho 1664end if 1665end function 1666 1667!!!------ Approximate form of DFT linear response kernel for close-shell, X(r1,r2), see Eq.3 of PCCP,14,3960. Only applied to the case wfntype=0 1668! r1 is taken as reference point and determined by refx,refy,refz. x,y,z in the argument is the coordinate of r2 1669real*8 function linrespkernel(x,y,z) 1670real*8 x,y,z,orbvalr1(nmo),orbvalr2(nmo) 1671call orbderv(1,1,nmo,refx,refy,refz,orbvalr1) 1672call orbderv(1,1,nmo,x,y,z,orbvalr2) 1673linrespkernel=0D0 1674do imo=1,nmo !Cycle occupied MOs 1675 if (nint(MOocc(imo))==2D0) then 1676 do jmo=idxHOMO+1,nmo !Cycle unoccupied MOs 1677 if (nint(MOocc(jmo))==0D0) linrespkernel=linrespkernel+orbvalr1(imo)*orbvalr1(jmo)*orbvalr2(jmo)*orbvalr2(imo)/(MOene(imo)-MOene(jmo)) 1678 end do 1679 end if 1680end do 1681linrespkernel=linrespkernel*4 1682end function 1683 1684!!!------------- Output Exchange-correlation density, correlation hole and correlation factor 1685!rfx,rfy,rfz is reference point (commonly use refx,refy,refz in defvar module), namely r1 1686!x,y,z in the argument is the coordinate of r2 1687!Calculate which function is controlled by "pairfunctype" in settings.ini, correlation type is determined by "paircorrtype" in settings.ini 1688real*8 function pairfunc(rfx,rfy,rfz,x,y,z) 1689real*8 rfx,rfy,rfz,x,y,z,orbvalr1(nmo),orbvalr2(nmo) 1690call orbderv(1,1,nmo,rfx,rfy,rfz,orbvalr1) 1691call orbderv(1,1,nmo,x,y,z,orbvalr2) 1692!Calculate alpha and beta density at r1 and r2 1693adensr1=0.0D0 1694bdensr1=0.0D0 1695adensr2=0.0D0 1696bdensr2=0.0D0 1697do i=1,nmo 1698 if (MOtype(i)==0) then 1699 adensr1=adensr1+MOocc(i)/2*orbvalr1(i)**2 1700 adensr2=adensr2+MOocc(i)/2*orbvalr2(i)**2 1701 bdensr1=bdensr1+MOocc(i)/2*orbvalr1(i)**2 1702 bdensr2=bdensr2+MOocc(i)/2*orbvalr2(i)**2 1703 else if (MOtype(i)==1) then 1704 adensr1=adensr1+MOocc(i)*orbvalr1(i)**2 1705 adensr2=adensr2+MOocc(i)*orbvalr2(i)**2 1706 else if (MOtype(i)==2) then 1707 bdensr1=bdensr1+MOocc(i)*orbvalr1(i)**2 1708 bdensr2=bdensr2+MOocc(i)*orbvalr2(i)**2 1709 end if 1710end do 1711totdensr1=adensr1+bdensr1 1712totdensr2=adensr2+bdensr2 1713 1714ntime=1 1715if (pairfunctype==12) ntime=2 !Will need both alpha and beta information (aXCdens and bXCdens), so process twice 1716do itime=1,ntime 1717 !Calculate exchange-correlation density first, and then calculate correlation hole and correlation factor 1718 !For RHF/ROHF, we calculate them as if present system is open-shell 1719 if (pairfunctype==1.or.pairfunctype==4.or.pairfunctype==7.or.pairfunctype==10.or.(pairfunctype==12.and.itime==1)) then !Set start and end index of alpha orbitals 1720 !Cycle alpha orbitals first to obtain aXCdens 1721 istart=1 1722 if (wfntype==0.or.wfntype==2.or.wfntype==3) then !RHF,ROHF,R-post-HF 1723 iend=nmo 1724 else if (wfntype==1.or.wfntype==4) then !UHF, U-post-HF 1725 do iend=nmo,1,-1 1726 if (MOtype(iend)==1) exit 1727 end do 1728 end if 1729 else if (pairfunctype==2.or.pairfunctype==5.or.pairfunctype==8.or.pairfunctype==11.or.(pairfunctype==12.and.itime==2)) then !Set start and end index of beta orbitals 1730 if (wfntype==0.or.wfntype==3) then !RHF,R-post-HF 1731 istart=1 1732 iend=nmo 1733 else if (wfntype==2) then !ROHF 1734 istart=1 1735 do iend=1,nmo 1736 if (MOtype(iend)==1) exit 1737 end do 1738 iend=iend-1 1739 else if (wfntype==1.or.wfntype==4) then !UHF, U-post-HF 1740 do istart=1,nmo 1741 if (MOtype(istart)==2) exit 1742 end do 1743 iend=nmo 1744 if (nint(nbelec)==0) iend=0 !less than istart, so below cycle will be skipped 1745 end if 1746 end if 1747 1748 XCtmp=0D0 !Really X+C density 1749 Xtmp=0D0 !Only X density 1750 Ctmp=0D0 !Only C density 1751 do i=istart,iend 1752 occi=MOocc(i) 1753 if (MOtype(i)==0) occi=occi/2D0 !Split close-shell orbital to spin orbital 1754 do j=istart,iend 1755 occj=MOocc(j) 1756 if (MOtype(j)==0) occj=occj/2D0 1757 tmpmul=orbvalr1(i)*orbvalr2(j)*orbvalr1(j)*orbvalr2(i) 1758 XCtmp=XCtmp-dsqrt(occi*occj)*tmpmul 1759 Xtmp=Xtmp-occi*occj*tmpmul 1760 Ctmp=Ctmp+(occi*occj-dsqrt(occi*occj))*tmpmul 1761 end do 1762 end do 1763 1764 if (pairfunctype==1.or.pairfunctype==4.or.pairfunctype==7.or.pairfunctype==10.or.(pairfunctype==12.and.itime==1)) then 1765 if (paircorrtype==1) aXCdens=Xtmp 1766 if (paircorrtype==2) aXCdens=Ctmp 1767 if (paircorrtype==3) aXCdens=XCtmp 1768 acorrhole=aXCdens/adensr1 1769 acorrfac=acorrhole/adensr2 1770 if (pairfunctype==1) pairfunc=acorrhole 1771 if (pairfunctype==4) pairfunc=acorrfac 1772 if (pairfunctype==7) pairfunc=aXCdens 1773 if (pairfunctype==10) pairfunc=adensr1*adensr2+aXCdens 1774 else if (pairfunctype==2.or.pairfunctype==5.or.pairfunctype==8.or.pairfunctype==11.or.(pairfunctype==12.and.itime==2)) then 1775 if (paircorrtype==1) bXCdens=Xtmp 1776 if (paircorrtype==2) bXCdens=Ctmp 1777 if (paircorrtype==3) bXCdens=XCtmp 1778 bcorrhole=bXCdens/bdensr1 1779 bcorrfac=bcorrhole/bdensr2 1780 if (pairfunctype==2) pairfunc=bcorrhole 1781 if (pairfunctype==5) pairfunc=bcorrfac 1782 if (pairfunctype==8) pairfunc=bXCdens 1783 if (pairfunctype==11) pairfunc=bdensr1*bdensr2+bXCdens 1784 end if 1785end do 1786if (pairfunctype==12) pairfunc=adensr1*(adensr2+bdensr2)+aXCdens +bdensr1*(adensr2+bdensr2)+bXCdens 1787end function 1788 1789 1790 1791!!!------------- Output source function 1792real*8 function srcfunc(x,y,z,imode) !Default imode=1 1793real*8 x,y,z,denomin 1794integer imode 1795denomin=4*pi*dsqrt((x-refx)**2+(y-refy)**2+(z-refz)**2) 1796if (denomin==0D0) denomin=0.001D0 1797if (imode==1) srcfunc=-flapl(x,y,z,'t')/denomin !Used to study effect of laplacian everywhere on specific point 1798if (imode==2) srcfunc=-flapl(refx,refy,refz,'t')/denomin !Used to study effect of laplacian at specific point on everywhere 1799end function 1800 1801 1802!!!------------------------- Output RDG with promolecular approximation 1803real*8 function RDGprodens(x,y,z) 1804real*8 x,y,z,elerho,elegrad(3) 1805call calchessmat_prodens(x,y,z,elerho,elegrad) 1806elegradnorm=dsqrt(sum(elegrad**2)) 1807if ((RDGprodens_maxrho/=0.0D0.and.elerho>=RDGprodens_maxrho)) then 1808 RDGprodens=100D0 1809else if (elegradnorm==0D0.or.elerho==0D0) then 1810 RDGprodens=999D0 1811else 1812 RDGprodens=0.161620459673995D0*elegradnorm/elerho**(4D0/3D0) 1813end if 1814end function 1815!!!----- Calculate Sign(lambda2(r))*rho(r) with promolecular approximation 1816!!! this is a shell used to convert subroutine to function form 1817real*8 function signlambda2rho_prodens(x,y,z) 1818real*8 x,y,z,sl2r,RDG 1819call signlambda2rho_RDG_prodens(x,y,z,sl2r,RDG) 1820signlambda2rho_prodens=sl2r 1821end function 1822!!!------ Calculate Sign(lambda2(r))*rho(r) and RDG at the same time with promolecular approximation 1823subroutine signlambda2rho_RDG_prodens(x,y,z,sl2r,RDG) 1824use util 1825real*8 x,y,z,elerho,RDG,sl2r,sumgrad2 1826real*8 eigvecmat(3,3),eigval(3),elehess(3,3),elegrad(3) 1827call calchessmat_prodens(x,y,z,elerho,elegrad,elehess) 1828call diagmat(elehess,eigvecmat,eigval,100,1D-6) 1829call sort(eigval) 1830if (eigval(2)/=0.0D0) then 1831 sl2r=elerho*eigval(2)/abs(eigval(2)) !At nuclei of single atom system, hessian returned may be zero matrix 1832else 1833 sl2r=-elerho !Around nuclei, eigval(2)/abs(eigval(2)) always be negative 1834end if 1835! sl2r=elerho !Only obtain promolecular density 1836sumgrad2=sum(elegrad(:)**2) 1837if ((RDGprodens_maxrho/=0.0D0.and.elerho>=RDGprodens_maxrho).or.elerho==0.0D0) then 1838 RDG=100D0 1839else if (sumgrad2==0D0.or.elerho==0D0) then 1840 RDG=999D0 1841else 1842 RDG=0.161620459673995D0*dsqrt(sumgrad2)/elerho**(4D0/3D0) 1843end if 1844end subroutine 1845 1846 1847!!!----- Calculate electron density, its gradient and Hessian matrix at x,y,z with promolecular approximation 1848!Notice that fragment must be properly defined!!! 1849subroutine calchessmat_prodens(xin,yin,zin,elerho,elegrad,elehess) 1850use util 1851real*8 elerho,xin,yin,zin 1852real*8,optional :: elegrad(3),elehess(3,3) 1853real*8 posarr(200),rhoarr(200) 1854elerho=0D0 1855derx=0D0 1856dery=0D0 1857derz=0D0 1858dxx=0D0 1859dyy=0D0 1860dzz=0D0 1861dxy=0D0 1862dyz=0D0 1863dxz=0D0 1864idohess=0 1865if (present(elehess)) idohess=1 1866do i=1,nfragatmnum 1867 iatm=fragatm(i) 1868 ind=a(iatm)%index 1869 rx=a(iatm)%x-xin !Relative x 1870 ry=a(iatm)%y-yin 1871 rz=a(iatm)%z-zin 1872 rx2=rx*rx 1873 ry2=ry*ry 1874 rz2=rz*rz 1875 r2=rx2+ry2+rz2 1876 if (atomdenscut==1) then !Tight cutoff, for CHNO corresponding to cutoff at rho=0.00001 1877 if (ind==1.and.r2>25D0) then !H, 6.63^2=43.9569. But this seems to be unnecessarily large, so I use 5^2=25 1878 cycle 1879 else if (ind==6.and.r2>58.6756D0) then !C, 7.66^2=58.6756 1880 cycle 1881 else if (ind==7.and.r2>43.917129D0) then !N, 6.627^2=43.917129 1882 cycle 1883 else if (ind==8.and.r2>34.9281D0) then !O, 5.91^2=34.9281 1884 cycle 1885 else if (r2>(2.5D0*vdwr(ind))**2) then !Other cases, larger than double of its vdw radius will be skipped 1886 cycle 1887 end if 1888 else if (atomdenscut==2) then !Medium cutoff, the result may be not as accurate as atomdenscut==1, but much more cheaper 1889 if (r2>(2.2D0*vdwr(ind))**2) cycle 1890 else if (atomdenscut==3) then !Loose cutoff, the most inaccurate 1891 if (r2>(1.8D0*vdwr(ind))**2) cycle 1892 else if (atomdenscut==4) then !Foolish cutoff, you need to know what you are doing 1893 if (r2>(1.5D0*vdwr(ind))**2) cycle 1894 end if 1895 r=dsqrt(r2) 1896 if (ind<=18) then !H~Ar 1897 r2_1d5=r2**1.5D0 1898 do j=1,3 1899 if (YWTatomcoeff(ind,j)==0D0) cycle 1900 expterm=YWTatomexp(ind,j) 1901 term=YWTatomcoeff(ind,j)*dexp(-r/expterm) 1902 elerho=elerho+term 1903 if (r==0D0) cycle !Derivative of STO at nuclei is pointless 1904 tmp=term/expterm/r 1905 derx=derx-tmp*rx !Calculating gradient doesn't cost detectable time, so always calculate it 1906 dery=dery-tmp*ry 1907 derz=derz-tmp*rz 1908 if (idohess==1) then 1909 tmp1=1/r2_1d5/expterm 1910 tmp2=1/r2/(expterm*expterm) 1911 dxx=dxx+term*(tmp1*rx2-1/r/expterm+tmp2*rx2) 1912 dyy=dyy+term*(tmp1*ry2-1/r/expterm+tmp2*ry2) 1913 dzz=dzz+term*(tmp1*rz2-1/r/expterm+tmp2*rz2) 1914 tmp=term*(tmp1+tmp2) 1915 dxy=dxy+rx*ry*tmp 1916 dyz=dyz+ry*rz*tmp 1917 dxz=dxz+rx*rz*tmp 1918 end if 1919 end do 1920 else !Heavier than Ar 1921! if (atomdenscut>=1.and.r>3*vdwr(ind)) cycle !Be careful, so use hentai criterion 1922 call genatmraddens(ind,posarr,rhoarr,npt) !Extract spherically averaged radial density of corresponding element 1923 call lagintpol(posarr(1:npt),rhoarr(1:npt),npt,r,term,der1r,der2r,3) 1924 elerho=elerho+term 1925 der1rdr=der1r/r 1926 derx=derx+der1rdr*rx 1927 dery=dery+der1rdr*ry 1928 derz=derz+der1rdr*rz 1929 if (idohess==1) then !I don't know how below code works, but it really works. See promolecular_grid routine in props.f90 of NCIPlot 1930 tmpval=(der2r-der1rdr)/r2 1931 dxx=dxx+der1rdr+tmpval*rx2 1932 dyy=dyy+der1rdr+tmpval*ry2 1933 dzz=dzz+der1rdr+tmpval*rz2 1934 dxy=dxy+tmpval*rx*ry 1935 dyz=dyz+tmpval*ry*rz 1936 dxz=dxz+tmpval*rx*rz 1937 end if 1938 end if 1939end do 1940if (present(elegrad)) then 1941 elegrad(1)=derx 1942 elegrad(2)=dery 1943 elegrad(3)=derz 1944end if 1945if (idohess==1) then 1946 elehess(1,1)=dxx 1947 elehess(2,2)=dyy 1948 elehess(3,3)=dzz 1949 elehess(1,2)=dxy 1950 elehess(2,3)=dyz 1951 elehess(1,3)=dxz 1952 elehess(2,1)=dxy 1953 elehess(3,2)=dyz 1954 elehess(3,1)=dxz 1955end if 1956end subroutine 1957 1958 1959!!---- Calculate atomic density based on STO fitted or radial density 1960!if indSTO==0, then all atom densities will be evaluated based on interpolation. if indSTO=18, then use STO fitted atomic density for element <18 1961real*8 function calcatmdens(iatm,x,y,z,indSTO) 1962use util 1963real*8 rho,x,y,z,posarr(200),rhoarr(200) 1964integer iatm,indSTO 1965calcatmdens=0 1966r=dsqrt( (a(iatm)%x-x)**2 + (a(iatm)%y-y)**2 + (a(iatm)%z-z)**2 ) 1967ind=a(iatm)%index 1968! if (r>6*vdwr(ind)) return !Doesn't improve speed evidently but deteriorate result in rare cases 1969if (ind<=indSTO) then !H~Ar, use STO fitted density. This is faster than using Lagrange interpolation technique, but not normalized to expected electron number 1970 do j=1,3 1971 if (YWTatomcoeff(ind,j)==0D0) cycle 1972 calcatmdens=calcatmdens+YWTatomcoeff(ind,j)*exp(-r/YWTatomexp(ind,j)) 1973 end do 1974else 1975 call genatmraddens(ind,posarr,rhoarr,npt) !Extract spherically averaged radial density of corresponding element 1976 call lagintpol(posarr(1:npt),rhoarr(1:npt),npt,r,calcatmdens,der1r,der2r,1) 1977end if 1978end function 1979!!---- Calculate promolecular density purely based on interpolation of radial density, the promolecular density obtained in this manner is quite accurate 1980real*8 function calcprodens(x,y,z,indSTO) 1981real*8 x,y,z 1982integer indSTO 1983calcprodens=0 1984do iatm=1,ncenter 1985 calcprodens=calcprodens+calcatmdens(iatm,x,y,z,indSTO) 1986end do 1987end function 1988 1989 1990 1991!!!--------------- Output Shannon information entropy function at a point 1992!itype=1 rho/N*ln(rho/N), this is normal definition 1993!itype=2 rho*ln(rho), this is Shannon information density, see J. Chem. Phys., 126, 191107 1994real*8 function infoentro(itype,x,y,z) 1995real*8 x,y,z,rho 1996integer itype 1997if (nelec==0D0) then 1998 infoentro=0D0 1999else 2000 rho=fdens(x,y,z) 2001 if (itype==1) rho=fdens(x,y,z)/nelec 2002 if (rho<=1D-100) then 2003 infoentro=0.0D0 2004 else 2005 infoentro=-rho*log(rho) 2006 end if 2007end if 2008end function 2009 2010 2011!!!------------------------- Output total ESP at a point 2012real*8 function totesp(x,y,z) 2013real*8 x,y,z 2014totesp=eleesp(x,y,z)+nucesp(x,y,z) 2015end function 2016 2017 2018!!!--------------- Output total ESP at a point, but skip the nucleus marked by variable "iskipnuc" 2019real*8 function totespskip(x,y,z,iskip) 2020real*8 x,y,z 2021integer iskip 2022totespskip=0 2023do i=1,ncenter 2024 if (i==iskip) cycle 2025 totespskip=totespskip+a(i)%charge/dsqrt((x-a(i)%x)**2+(y-a(i)%y)**2+(z-a(i)%z)**2) 2026end do 2027totespskip=totespskip+eleesp(x,y,z) 2028end function 2029 2030!!!------------------------- Output ESP from nuclear or atomic charges at a point 2031!At nuclear positions, this function returns 1000 instead of infinity to avoid numerical problems 2032real*8 function nucesp(x,y,z) 2033nucesp=0D0 2034do i=1,nfragatmnum 2035 dist2mpx=(x-a(fragatm(i))%x)**2 2036 dist2mpy=(y-a(fragatm(i))%y)**2 2037 dist2mpz=(z-a(fragatm(i))%z)**2 2038 dist2=dist2mpx+dist2mpy+dist2mpz 2039 if (dist2==0D0) then 2040 nucesp=1D3 2041 return 2042 end if 2043 nucesp=nucesp+a(fragatm(i))%charge/dsqrt(dist2) 2044end do 2045end function 2046 2047 2048!------------------------- Calculate ESP from electrons at a point 2049!Note that the negative sign of electron is already taken into account 2050real*8 function eleesp(Cx,Cy,Cz) 2051implicit none 2052integer,parameter :: narrmax=396 !Enough for h-type GTF, 5+5=10. Alri(0:10,0:5,0:5)-->11*6*6=396 2053real*8 Cx,Cy,Cz,term,ep,Ax,Ay,Az,Bx,By,Bz,Aexp,Bexp,Px,Py,Pz,prefac,tmpval 2054real*8 sqPC,sqAB,expngaPC,PAx,PAy,PAz,PBx,PBy,PBz,PCx,PCy,PCz,fjtmp,addesp 2055real*8 Alri(narrmax),Amsj(narrmax),Antk(narrmax),Fn(0:10) !Enough for h-type GTF, 5+5=10 2056real*8 twoepsqPC,tl,tm,tn,espprivate,espexpcut 2057integer nu,imo,iprim,jprim,maxFn,maplri(narrmax),mapmsj(narrmax),mapntk(narrmax),tmpnuml,tmpnumm,tmpnumn 2058integer Aix,Aiy,Aiz,Bix,Biy,Biz,l,r,i,m,s,j,n,t,k,icen,jcen,sumAi,sumBi 2059eleesp=0.0D0 2060espexpcut=log10(espprecutoff)*3 2061nthreads=getNThreads() 2062!$OMP parallel do private(Aix,Aiy,Aiz,Bix,Biy,Biz,l,r,i,m,s,j,n,t,k,icen,jcen,sumAi,sumBi, & 2063!$OMP nu,imo,iprim,jprim,maxFn,maplri,mapmsj,mapntk,tmpnuml,tmpnumm,tmpnumn,& 2064!$OMP twoepsqPC,tl,tm,tn,Alri,Amsj,Antk,Fn,& 2065!$OMP sqPC,sqAB,expngaPC,PAx,PAy,PAz,PBx,PBy,PBz,PCx,PCy,PCz,fjtmp,addesp,& 2066!$OMP term,ep,Ax,Ay,Az,Bx,By,Bz,Aexp,Bexp,Px,Py,Pz,prefac,tmpval,espprivate) shared(eleesp) schedule(dynamic) NUM_THREADS(nthreads) 2067do iprim=1,nprims 2068 espprivate=0D0 2069 icen=b(iprim)%center 2070 Aexp=b(iprim)%exp 2071 Ax=a(icen)%x 2072 Ay=a(icen)%y 2073 Az=a(icen)%z 2074 Aix=type2ix(b(iprim)%functype) 2075 Aiy=type2iy(b(iprim)%functype) 2076 Aiz=type2iz(b(iprim)%functype) 2077 sumAi=Aix+Aiy+Aiz 2078 do jprim=iprim,nprims 2079 jcen=b(jprim)%center 2080 Bexp=b(jprim)%exp 2081 Bix=type2ix(b(jprim)%functype) 2082 Biy=type2iy(b(jprim)%functype) 2083 Biz=type2iz(b(jprim)%functype) 2084 Bx=a(jcen)%x 2085 By=a(jcen)%y 2086 Bz=a(jcen)%z 2087 sumBi=Bix+Biy+Biz 2088 ep=Aexp+Bexp 2089 Px=(Ax*Aexp+Bx*Bexp)/ep 2090 Py=(Ay*Aexp+By*Bexp)/ep 2091 Pz=(Az*Aexp+Bz*Bexp)/ep 2092 PAx=Px-Ax 2093 PAy=Py-Ay 2094 PAz=Pz-Az 2095 PBx=Px-Bx 2096 PBy=Py-By 2097 PBz=Pz-Bz 2098 sqAB=(Ax-Bx)**2+(Ay-By)**2+(Az-Bz)**2 2099 PCx=Px-Cx 2100 PCy=Py-Cy 2101 PCz=Pz-Cz 2102 sqPC=PCx*PCx+PCy*PCy+PCz*PCz 2103 2104 tmpval=-Aexp*Bexp*sqAB/ep 2105 prefac=2*pi/ep*dexp(tmpval) 2106 2107 if (-ep*sqPC>espexpcut) then 2108 expngaPC=dexp(-ep*sqPC) 2109 else 2110 expngaPC=0 2111 end if 2112 maxFn=sumAi+sumBi 2113 Fn(maxFn)=Fmch(maxFn,ep*sqPC,expngaPC) 2114 nu=maxFn 2115 twoepsqPC=2*ep*sqPC 2116 do while (nu>0) 2117 Fn(nu-1)=(expngaPC+twoepsqPC*Fn(nu))/(2*nu-1) !cook book p280 2118 nu=nu-1 2119 end do 2120 2121 tmpnuml=0 2122 do l=0,Aix+Bix 2123 tl=1.0D0 2124 if (mod(l,2)==1) tl=-1.0D0 2125 fjtmp=fj(l,Aix,Bix,PAx,PBx)*tl*fact(l) 2126 do r=0,l/2.0D0 2127 do i=0,(l-2*r)/2.0D0 2128 tmpnuml=tmpnuml+1 2129 Alri(tmpnuml)=Afac(l,r,i,PCx,ep,fjtmp) 2130 maplri(tmpnuml)=l-2*r-i 2131 end do 2132 end do 2133 end do 2134 2135 tmpnumm=0 2136 do m=0,Aiy+Biy 2137 tm=1.0D0 2138 if (mod(m,2)==1) tm=-1.0D0 2139 fjtmp=fj(m,Aiy,Biy,PAy,PBy)*tm*fact(m) 2140 do s=0,m/2.0D0 2141 do j=0,(m-2*s)/2.0D0 2142 tmpnumm=tmpnumm+1 2143 Amsj(tmpnumm)=Afac(m,s,j,PCy,ep,fjtmp) 2144 mapmsj(tmpnumm)=m-2*s-j 2145 end do 2146 end do 2147 end do 2148 2149 tmpnumn=0 2150 do n=0,Aiz+Biz 2151 tn=1.0D0 2152 if (mod(n,2)==1) tn=-1.0D0 2153 fjtmp=fj(n,Aiz,Biz,PAz,PBz)*tn*fact(n) 2154 do t=0,n/2.0D0 2155 do k=0,(n-2*t)/2.0D0 2156 tmpnumn=tmpnumn+1 2157 Antk(tmpnumn)=Afac(n,t,k,PCz,ep,fjtmp) 2158 mapntk(tmpnumn)=n-2*t-k 2159 end do 2160 end do 2161 end do 2162 2163 term=0.0D0 2164 !Now calc "term"=<psi(iprim)|1/r_Z|psi(jprim)> 2165 do l=1,tmpnuml 2166 do m=1,tmpnumm 2167 do n=1,tmpnumn 2168 term=term+Alri(l)*Amsj(m)*Antk(n)*Fn(maplri(l)+mapmsj(m)+mapntk(n)) 2169 end do 2170 end do 2171 end do 2172 2173 if (iprim/=jprim) term=2.0*term 2174 term=term*prefac 2175 addesp=0.0D0 2176 do imo=1,nmo 2177 addesp=addesp+MOocc(imo)*CO(imo,iprim)*CO(imo,jprim) 2178 end do 2179 espprivate=espprivate+addesp*term 2180 end do !end j primitive 2181!$OMP critical 2182 eleesp=eleesp+espprivate 2183!$OMP end critical 2184end do !end i primitive 2185!$OMP end parallel do 2186eleesp=-eleesp 2187end function 2188 2189 2190!!!------------------------- Calculate ESP in a plane 2191!maxnumgrid is the maximum value of ngridnum1 and ngridnum2, this value determine the size of Alrivec,Amsjvec,Antkvec 2192!We don't allocate Alrivec,Amsjvec,Antkvec dynamically since if we do such thing, this routine will crash in win7-64bit system. 2193!The reason may be that dynamical arrays are not fully compatiable with private property in OpenMP 2194subroutine planeesp(maxnumgrid) 2195implicit none 2196integer maxnumgrid 2197integer,parameter :: narrmax=396 !Enough for h-type GTF, 5+5=10. Alri(0:10,0:5,0:5)-->11*6*6=396 2198real*8 Cx,Cy,Cz,Cxold,Cyold,Czold,term,ep,Ax,Ay,Az,Bx,By,Bz,Aexp,Bexp,Px,Py,Pz,prefac,tmpval 2199real*8 sqPC,sqAB,expngaPC,PAx,PAy,PAz,PBx,PBy,PBz,PCx,PCy,PCz,fjtmp,addesp,espexpcut 2200real*8 Alri(narrmax),Amsj(narrmax),Antk(narrmax),Fnmat(0:ngridnum1-1,0:ngridnum2-1),Fnvec(0:10) !Enough for h-type GTF, 5+5=10 2201real*8:: Alrivec(narrmax,maxnumgrid),Amsjvec(narrmax,maxnumgrid),Antkvec(narrmax,maxnumgrid) 2202real*8 twoepsqPC,tl,tm,tn,pleprivate(ngridnum1,ngridnum2) !Store plane contribution of GTFs in each thread, then sum up 2203integer nu,imo,iprim,jprim,maxFn,maplri(narrmax),mapmsj(narrmax),mapntk(narrmax),tmpnuml,tmpnumm,tmpnumn 2204integer Aix,Aiy,Aiz,Bix,Biy,Biz,l,r,i,m,s,j,n,t,k,icen,jcen,sumAi,sumBi,ii,jj,planetype,numx,numy,numz,ifinish 2205ifinish=0 2206espexpcut=log10(espprecutoff)*3 2207planemat=0D0 2208!Calc ESP of nuclear contribution 2209do ii=0,ngridnum1-1 2210 do jj=0,ngridnum2-1 2211 Cx=orgx2D+ii*v1x+jj*v2x 2212 Cy=orgy2D+ii*v1y+jj*v2y 2213 Cz=orgz2D+ii*v1z+jj*v2z 2214 planemat(ii+1,jj+1)=nucesp(Cx,Cy,Cz) 2215 end do 2216end do 2217 2218nthreads=getNThreads() 2219!$OMP PARALLEL DO private(Aix,Aiy,Aiz,Bix,Biy,Biz,l,r,i,m,s,j,n,t,k,icen,jcen,sumAi,sumBi,ii,jj,planetype,numx,numy,numz,& 2220!$OMP nu,imo,iprim,jprim,maxFn,maplri,mapmsj,mapntk,tmpnuml,tmpnumm,tmpnumn,& 2221!$OMP twoepsqPC,tl,tm,tn,Alrivec,Amsjvec,Antkvec,Alri,Amsj,Antk,Fnmat,Fnvec,& 2222!$OMP sqPC,sqAB,expngaPC,PAx,PAy,PAz,PBx,PBy,PBz,PCx,PCy,PCz,fjtmp,addesp,& 2223!$OMP Cx,Cy,Cz,Cxold,Cyold,Czold,term,ep,Ax,Ay,Az,Bx,By,Bz,Aexp,Bexp,Px,Py,Pz,prefac,tmpval,pleprivate) & 2224!$OMP shared(planemat,ifinish) schedule(dynamic) NUM_THREADS(nthreads) 2225do iprim=1,nprims 2226 pleprivate=0D0 2227 icen=b(iprim)%center 2228 Aexp=b(iprim)%exp 2229 Ax=a(icen)%x 2230 Ay=a(icen)%y 2231 Az=a(icen)%z 2232 Aix=type2ix(b(iprim)%functype) 2233 Aiy=type2iy(b(iprim)%functype) 2234 Aiz=type2iz(b(iprim)%functype) 2235 sumAi=Aix+Aiy+Aiz 2236 do jprim=iprim,nprims 2237 jcen=b(jprim)%center 2238 Bexp=b(jprim)%exp 2239 Bix=type2ix(b(jprim)%functype) 2240 Biy=type2iy(b(jprim)%functype) 2241 Biz=type2iz(b(jprim)%functype) 2242 Bx=a(jcen)%x 2243 By=a(jcen)%y 2244 Bz=a(jcen)%z 2245 ep=Aexp+Bexp 2246 Px=(Ax*Aexp+Bx*Bexp)/ep 2247 Py=(Ay*Aexp+By*Bexp)/ep 2248 Pz=(Az*Aexp+Bz*Bexp)/ep 2249 PAx=Px-Ax 2250 PAy=Py-Ay 2251 PAz=Pz-Az 2252 PBx=Px-Bx 2253 PBy=Py-By 2254 PBz=Pz-Bz 2255 sqAB=(Ax-Bx)**2+(Ay-By)**2+(Az-Bz)**2 2256 tmpval=-Aexp*Bexp*sqAB/ep 2257 prefac=2*pi/ep*dexp(tmpval) 2258 if (prefac<espprecutoff) cycle 2259 2260 Cxold=999.99912D0 !An arbitrary number 2261 Cyold=999.99912D0 2262 Czold=999.99912D0 2263 sumBi=Bix+Biy+Biz 2264 maxFn=sumAi+sumBi 2265 2266 !! Start cycle grid point 2267 do ii=0,ngridnum1-1 2268 do jj=0,ngridnum2-1 2269 Cx=orgx2D+ii*v1x+jj*v2x 2270 Cy=orgy2D+ii*v1y+jj*v2y 2271 Cz=orgz2D+ii*v1z+jj*v2z 2272 PCx=Px-Cx 2273 PCy=Py-Cy 2274 PCz=Pz-Cz 2275 sqPC=PCx*PCx+PCy*PCy+PCz*PCz 2276 twoepsqPC=2*ep*sqPC 2277 term=0.0D0 2278 if (-ep*sqPC>espexpcut) then 2279 expngaPC=dexp(-ep*sqPC) 2280 else 2281 expngaPC=0D0 2282 end if 2283 Fnmat(ii,jj)=Fmch(maxFn,ep*sqPC,expngaPC) 2284 Fnvec(maxFn)=Fnmat(ii,jj) 2285 do nu=maxFn,1,-1 2286 Fnvec(nu-1)=(expngaPC+twoepsqPC*Fnvec(nu))/(2*nu-1) !cook book p280 2287 end do 2288 2289 if (Cx/=Cxold) then 2290 tmpnuml=0 2291 do l=0,Aix+Bix 2292 tl=1.0D0 2293 if (mod(l,2)==1) tl=-1.0D0 2294 fjtmp=fj(l,Aix,Bix,PAx,PBx)*tl*fact(l) 2295 do r=0,l/2.0D0 2296 do i=0,(l-2*r)/2.0D0 2297 tmpnuml=tmpnuml+1 2298 Alri(tmpnuml)=Afac(l,r,i,PCx,ep,fjtmp) 2299 maplri(tmpnuml)=l-2*r-i 2300 end do 2301 end do 2302 end do 2303 Cxold=Cx 2304 end if 2305 if (Cy/=Cyold) then 2306 tmpnumm=0 2307 do m=0,Aiy+Biy 2308 tm=1.0D0 2309 if (mod(m,2)==1) tm=-1.0D0 2310 fjtmp=fj(m,Aiy,Biy,PAy,PBy)*tm*fact(m) 2311 do s=0,m/2.0D0 2312 do j=0,(m-2*s)/2.0D0 2313 tmpnumm=tmpnumm+1 2314 Amsj(tmpnumm)=Afac(m,s,j,PCy,ep,fjtmp) 2315 mapmsj(tmpnumm)=m-2*s-j 2316 end do 2317 end do 2318 end do 2319 Cyold=Cy 2320 end if 2321 if (Cz/=Czold) then 2322 tmpnumn=0 2323 do n=0,Aiz+Biz 2324 tn=1.0D0 2325 if (mod(n,2)==1) tn=-1.0D0 2326 fjtmp=fj(n,Aiz,Biz,PAz,PBz)*tn*fact(n) 2327 do t=0,n/2.0D0 2328 do k=0,(n-2*t)/2.0D0 2329 tmpnumn=tmpnumn+1 2330 Antk(tmpnumn)=Afac(n,t,k,PCz,ep,fjtmp) 2331 mapntk(tmpnumn)=n-2*t-k 2332 end do 2333 end do 2334 end do 2335 Czold=Cz 2336 end if 2337 2338 !Now calc "term"=<psi(iprim)|1/r_Z|psi(jprim)> 2339 do l=1,tmpnuml 2340 do m=1,tmpnumm 2341 do n=1,tmpnumn 2342 term=term+Alri(l)*Amsj(m)*Antk(n)*Fnvec(maplri(l)+mapmsj(m)+mapntk(n)) 2343 end do 2344 end do 2345 end do 2346 2347 if (iprim/=jprim) term=2.0*term 2348 term=term*prefac 2349 addesp=0.0D0 2350 do imo=1,nmo 2351 addesp=addesp+MOocc(imo)*CO(imo,iprim)*CO(imo,jprim) 2352 end do 2353 pleprivate(ii+1,jj+1)=pleprivate(ii+1,jj+1)-addesp*term 2354 end do !end jj cycle 2355 end do !end ii cycle 2356 2357 end do !end j primitive 2358 ifinish=ifinish+1 2359 write(*,"(' Finished: ',i6,' /',i6)") ifinish,nprims 2360!$OMP CRITICAL 2361 planemat=planemat+pleprivate 2362!$OMP END CRITICAL 2363end do !end i primitive 2364!$OMP END PARALLEL DO 2365end subroutine 2366 2367 2368!!!---------------- Calculate grid data of ESP from electrons and accumulate to cubmat (which may already records nuclear ESP) 2369subroutine espcub 2370implicit none 2371integer,parameter :: narrmax=396 !Enough for h-type GTF, 5+5=10. Alri(0:10,0:5,0:5)-->11*6*6=396 2372real*8 Cx,Cy,Cz,term,ep,Ax,Ay,Az,Bx,By,Bz,Aexp,Bexp,Px,Py,Pz,prefac,tmpval 2373real*8 sqPC,sqAB,expngaPC,PAx,PAy,PAz,PBx,PBy,PBz,PCx,PCy,PCz,fjtmp,addesp 2374real*8 Alrivec(narrmax,nx),Amsjvec(narrmax,ny),Antkvec(narrmax,nz),Fnmat(0:nx-1,0:ny-1,0:nz-1),Fnvec(0:10) !Enough for h-type GTF, 5+5=10 2375real*8 twoepsqPC,tl,tm,tn,time_begin,time_end,espexpcut,cubprivate(nx,ny,nz) 2376integer nu,imo,iprim,jprim,maxFn,maplri(narrmax),mapmsj(narrmax),mapntk(narrmax),tmpnuml,tmpnumm,tmpnumn 2377integer Aix,Aiy,Aiz,Bix,Biy,Biz,l,r,i,m,s,j,n,t,k,icen,jcen,sumAi,sumBi,ii,jj,kk,ifinish 2378CALL CPU_TIME ( time_begin ) 2379ifinish=0 2380espexpcut=log10(espprecutoff)*3 2381nthreads=getNThreads() 2382!$OMP PARALLEL DO private(Aix,Aiy,Aiz,Bix,Biy,Biz,l,r,i,m,s,j,n,t,k,icen,jcen,sumAi,sumBi,ii,jj,kk, & 2383!$OMP nu,imo,iprim,jprim,maxFn,maplri,mapmsj,mapntk,tmpnuml,tmpnumm,tmpnumn, & 2384!$OMP twoepsqPC,tl,tm,tn,Alrivec,Amsjvec,Antkvec,Fnmat,Fnvec, & 2385!$OMP sqPC,sqAB,expngaPC,PAx,PAy,PAz,PBx,PBy,PBz,PCx,PCy,PCz,fjtmp,addesp,cubprivate, & 2386!$OMP Cx,Cy,Cz,term,ep,Ax,Ay,Az,Bx,By,Bz,Aexp,Bexp,Px,Py,Pz,tmpval,prefac) & 2387!$OMP shared(cubmat,ifinish) schedule(dynamic) NUM_THREADS(nthreads) 2388do iprim=1,nprims 2389 cubprivate=0D0 2390 icen=b(iprim)%center 2391 Aexp=b(iprim)%exp 2392 Ax=a(icen)%x 2393 Ay=a(icen)%y 2394 Az=a(icen)%z 2395 Aix=type2ix(b(iprim)%functype) 2396 Aiy=type2iy(b(iprim)%functype) 2397 Aiz=type2iz(b(iprim)%functype) 2398 sumAi=Aix+Aiy+Aiz 2399 do jprim=iprim,nprims 2400 jcen=b(jprim)%center 2401 Bexp=b(jprim)%exp 2402 Bix=type2ix(b(jprim)%functype) 2403 Biy=type2iy(b(jprim)%functype) 2404 Biz=type2iz(b(jprim)%functype) 2405 Bx=a(jcen)%x 2406 By=a(jcen)%y 2407 Bz=a(jcen)%z 2408 ep=Aexp+Bexp 2409 Px=(Ax*Aexp+Bx*Bexp)/ep 2410 Py=(Ay*Aexp+By*Bexp)/ep 2411 Pz=(Az*Aexp+Bz*Bexp)/ep 2412 PAx=Px-Ax 2413 PAy=Py-Ay 2414 PAz=Pz-Az 2415 PBx=Px-Bx 2416 PBy=Py-By 2417 PBz=Pz-Bz 2418 sqAB=(Ax-Bx)**2+(Ay-By)**2+(Az-Bz)**2 2419 tmpval=-Aexp*Bexp*sqAB/ep 2420 2421 prefac=2*pi/ep*dexp(tmpval) 2422 if (prefac<espprecutoff) cycle 2423 sumBi=Bix+Biy+Biz 2424 maxFn=sumAi+sumBi 2425 2426 do ii=1,nx 2427 Cx=orgx+(ii-1)*dx 2428 PCx=Px-Cx 2429 tmpnuml=0 2430 do l=0,Aix+Bix 2431 tl=1.0D0 2432 if (mod(l,2)==1) tl=-1.0D0 2433 fjtmp=fj(l,Aix,Bix,PAx,PBx)*tl*fact(l) 2434 do r=0,l/2.0D0 2435 do i=0,(l-2*r)/2.0D0 2436 tmpnuml=tmpnuml+1 2437 Alrivec(tmpnuml,ii)=Afac(l,r,i,PCx,ep,fjtmp) 2438 maplri(tmpnuml)=l-2*r-i 2439 end do 2440 end do 2441 end do 2442 end do 2443 do ii=1,ny 2444 Cy=orgy+(ii-1)*dy 2445 PCy=Py-Cy 2446 tmpnumm=0 2447 do m=0,Aiy+Biy 2448 tm=1.0D0 2449 if (mod(m,2)==1) tm=-1.0D0 2450 fjtmp=fj(m,Aiy,Biy,PAy,PBy)*tm*fact(m) 2451 do s=0,m/2.0D0 2452 do j=0,(m-2*s)/2.0D0 2453 tmpnumm=tmpnumm+1 2454 Amsjvec(tmpnumm,ii)=Afac(m,s,j,PCy,ep,fjtmp) 2455 mapmsj(tmpnumm)=m-2*s-j 2456 end do 2457 end do 2458 end do 2459 end do 2460 do ii=1,nz 2461 Cz=orgz+(ii-1)*dz 2462 PCz=Pz-Cz 2463 tmpnumn=0 2464 do n=0,Aiz+Biz 2465 tn=1.0D0 2466 if (mod(n,2)==1) tn=-1.0D0 2467 fjtmp=fj(n,Aiz,Biz,PAz,PBz)*tn*fact(n) 2468 do t=0,n/2.0D0 2469 do k=0,(n-2*t)/2.0D0 2470 tmpnumn=tmpnumn+1 2471 Antkvec(tmpnumn,ii)=Afac(n,t,k,PCz,ep,fjtmp) 2472 mapntk(tmpnumn)=n-2*t-k 2473 end do 2474 end do 2475 end do 2476 end do 2477 2478 !! Start cycle grid point 2479 do kk=0,nz-1 2480 do jj=0,ny-1 2481 do ii=0,nx-1 2482 Cx=orgx+ii*dx 2483 Cy=orgy+jj*dy 2484 Cz=orgz+kk*dz 2485 PCx=Px-Cx 2486 PCy=Py-Cy 2487 PCz=Pz-Cz 2488 sqPC=PCx*PCx+PCy*PCy+PCz*PCz 2489 twoepsqPC=2*ep*sqPC 2490 term=0.0D0 2491 if (-ep*sqPC>espexpcut) then 2492 expngaPC=dexp(-ep*sqPC) 2493 else 2494 expngaPC=0D0 2495 end if 2496 Fnmat(ii,jj,kk)=Fmch(maxFn,ep*sqPC,expngaPC) 2497 Fnvec(maxFn)=Fnmat(ii,jj,kk) 2498 do nu=maxFn,1,-1 2499 Fnvec(nu-1)=(expngaPC+twoepsqPC*Fnvec(nu))/(2*nu-1) !cook book p280 2500 end do 2501 do l=1,tmpnuml 2502 do m=1,tmpnumm 2503 do n=1,tmpnumn 2504 term=term+Alrivec(l,ii+1)*Amsjvec(m,jj+1)*Antkvec(n,kk+1)*Fnvec(maplri(l)+mapmsj(m)+mapntk(n)) 2505 end do 2506 end do 2507 end do 2508 2509 if (iprim/=jprim) term=2D0*term 2510 term=term*prefac 2511 addesp=0.0D0 2512 do imo=1,nmo 2513 addesp=addesp+MOocc(imo)*CO(imo,iprim)*CO(imo,jprim) 2514 end do 2515 cubprivate(ii+1,jj+1,kk+1)=cubprivate(ii+1,jj+1,kk+1)-addesp*term 2516 end do !end ii cycle 2517 end do !end jj cycle 2518 end do !end kk cycle 2519 2520 end do !end j primitive 2521 ifinish=ifinish+1 2522 write(*,"(' Finished: ',i6,'/',i6)") ifinish,nprims 2523!$OMP CRITICAL 2524 cubmat=cubmat+cubprivate 2525!$OMP END CRITICAL 2526end do !end i primitive 2527!$OMP END PARALLEL DO 2528CALL CPU_TIME ( time_end ) 2529if (isys==1) write(*,"(' ESP calculation took up CPU time',f12.2,' seconds')") time_end-time_begin 2530end subroutine 2531 2532 2533 2534 2535 2536 2537!======================================================================== 2538!============== Utilities routine for function modeule ================== 2539!======================================================================== 2540 2541!!!---------- Generate nuclear attraction potential integral matrix between all GTFs at a given point 2542subroutine genGTFattmat(x,y,z,GTFattmat) 2543use defvar 2544integer,parameter :: narrmax=396 2545real*8 x,y,z,GTFattmat(nprims,nprims),Alri(narrmax),Amsj(narrmax),Antk(narrmax),Fn(0:10) !Enough for h-type GTF, 5+5=10. Alri(0:10,0:5,0:5)-->11*6*6=396 2546integer maxFn,maplri(narrmax),mapmsj(narrmax),mapntk(narrmax),tmpnuml,tmpnumm,tmpnumn,Aix,Aiy,Aiz,Bix,Biy,Biz,r,s,t,sumAi,sumBi 2547espexpcut=log10(espprecutoff)*3 2548nthreads=getNThreads() 2549!$OMP parallel do private(Aix,Aiy,Aiz,Bix,Biy,Biz,l,r,i,m,s,j,n,t,k,icen,jcen,sumAi,sumBi, & 2550!$OMP nu,iprim,jprim,maxFn,maplri,mapmsj,mapntk,tmpnuml,tmpnumm,tmpnumn,& 2551!$OMP twoepsqPC,tl,tm,tn,Alri,Amsj,Antk,Fn,sqPC,sqAB,expngaPC,PAx,PAy,PAz,PBx,PBy,PBz,PCx,PCy,PCz,fjtmp,& 2552!$OMP term,ep,Ax,Ay,Az,Bx,By,Bz,Aexp,Bexp,Px,Py,Pz,prefac,tmpval) shared(GTFattmat) schedule(dynamic) NUM_THREADS(nthreads) 2553do iprim=1,nprims 2554 icen=b(iprim)%center 2555 Aexp=b(iprim)%exp 2556 Ax=a(icen)%x 2557 Ay=a(icen)%y 2558 Az=a(icen)%z 2559 Aix=type2ix(b(iprim)%functype) 2560 Aiy=type2iy(b(iprim)%functype) 2561 Aiz=type2iz(b(iprim)%functype) 2562 sumAi=Aix+Aiy+Aiz 2563 do jprim=iprim,nprims 2564 jcen=b(jprim)%center 2565 Bexp=b(jprim)%exp 2566 Bix=type2ix(b(jprim)%functype) 2567 Biy=type2iy(b(jprim)%functype) 2568 Biz=type2iz(b(jprim)%functype) 2569 Bx=a(jcen)%x 2570 By=a(jcen)%y 2571 Bz=a(jcen)%z 2572 sumBi=Bix+Biy+Biz 2573 ep=Aexp+Bexp 2574 Px=(Ax*Aexp+Bx*Bexp)/ep 2575 Py=(Ay*Aexp+By*Bexp)/ep 2576 Pz=(Az*Aexp+Bz*Bexp)/ep 2577 PAx=Px-Ax 2578 PAy=Py-Ay 2579 PAz=Pz-Az 2580 PBx=Px-Bx 2581 PBy=Py-By 2582 PBz=Pz-Bz 2583 sqAB=(Ax-Bx)**2+(Ay-By)**2+(Az-Bz)**2 2584 PCx=Px-x 2585 PCy=Py-y 2586 PCz=Pz-z 2587 sqPC=PCx*PCx+PCy*PCy+PCz*PCz 2588 tmpval=-Aexp*Bexp*sqAB/ep 2589 prefac=2*pi/ep*dexp(tmpval) 2590 if (-ep*sqPC>espexpcut) then 2591 expngaPC=dexp(-ep*sqPC) 2592 else 2593 expngaPC=0 2594 end if 2595 maxFn=sumAi+sumBi 2596 Fn(maxFn)=Fmch(maxFn,ep*sqPC,expngaPC) 2597 nu=maxFn 2598 twoepsqPC=2*ep*sqPC 2599 do while (nu>0) 2600 Fn(nu-1)=(expngaPC+twoepsqPC*Fn(nu))/(2*nu-1) !cook book p280 2601 nu=nu-1 2602 end do 2603 2604 tmpnuml=0 2605 do l=0,Aix+Bix 2606 tl=1.0D0 2607 if (mod(l,2)==1) tl=-1.0D0 2608 fjtmp=fj(l,Aix,Bix,PAx,PBx)*tl*fact(l) 2609 do r=0,l/2.0D0 2610 do i=0,(l-2*r)/2.0D0 2611 tmpnuml=tmpnuml+1 2612 Alri(tmpnuml)=Afac(l,r,i,PCx,ep,fjtmp) 2613 maplri(tmpnuml)=l-2*r-i 2614 end do 2615 end do 2616 end do 2617 tmpnumm=0 2618 do m=0,Aiy+Biy 2619 tm=1.0D0 2620 if (mod(m,2)==1) tm=-1.0D0 2621 fjtmp=fj(m,Aiy,Biy,PAy,PBy)*tm*fact(m) 2622 do s=0,m/2.0D0 2623 do j=0,(m-2*s)/2.0D0 2624 tmpnumm=tmpnumm+1 2625 Amsj(tmpnumm)=Afac(m,s,j,PCy,ep,fjtmp) 2626 mapmsj(tmpnumm)=m-2*s-j 2627 end do 2628 end do 2629 end do 2630 tmpnumn=0 2631 do n=0,Aiz+Biz 2632 tn=1.0D0 2633 if (mod(n,2)==1) tn=-1.0D0 2634 fjtmp=fj(n,Aiz,Biz,PAz,PBz)*tn*fact(n) 2635 do t=0,n/2.0D0 2636 do k=0,(n-2*t)/2.0D0 2637 tmpnumn=tmpnumn+1 2638 Antk(tmpnumn)=Afac(n,t,k,PCz,ep,fjtmp) 2639 mapntk(tmpnumn)=n-2*t-k 2640 end do 2641 end do 2642 end do 2643 2644 term=0D0 2645 !Now calc "term"=<psi(iprim)|1/r12|psi(jprim)>, r1 is inputted x,y,z, r2 is the integration variable 2646 do l=1,tmpnuml 2647 do m=1,tmpnumm 2648 do n=1,tmpnumn 2649 term=term+Alri(l)*Amsj(m)*Antk(n)*Fn(maplri(l)+mapmsj(m)+mapntk(n)) 2650 end do 2651 end do 2652 end do 2653 term=term*prefac 2654 GTFattmat(iprim,jprim)=term 2655 GTFattmat(jprim,iprim)=term 2656 end do !end j primitive 2657end do !end i primitive 2658!$OMP end parallel do 2659end subroutine 2660 2661!!!------------------------- Calculate A-factor at Cook book p245 2662real*8 function Afac(l,r,i,PC,gamma,fjtmp) 2663integer l,r,i,ti 2664real*8 gamma,PC,comp,PCterm,fjtmp 2665ti=1.0D0 2666if (mod(i,2)==1) ti=-1.0D0 !faster than ti=(-1)**i 2667PCterm=1.0D0 2668if (l-2*r-2*i/=0) PCterm=PC**(l-2*r-2*i) 2669comp=ti*PCterm*(0.25D0/gamma)**(r+i) / ( fact(r)*fact(i)*fact(l-2*r-2*i) ) 2670Afac=fjtmp*comp 2671! comp=(-1)**i*fact(l)*PC**(l-2*r-2*i)*(1/(4*gamma))**(r+i) / ( fact(r)*fact(i)*fact(l-2*r-2*i) ) 2672! Afac=(-1)**l * fj(l,l1,l2,PA,PB)*comp 2673end function 2674 2675!!!------------------------- Calculate fj at Cook book p237 2676real*8 function fj(j,l,m,aa,bb) 2677real*8 aa,bb,pre,at,bt 2678integer j,l,m,k,imin,imax 2679imax=min(j,l) 2680imin=max(0,j-m) 2681fj=0.0D0 2682do k=imin,imax 2683 pre=fact(l)/fact(l-k)/fact(k) * fact(m)/fact(m-j+k)/fact(j-k) 2684 at=1.0D0 2685 bt=1.0D0 2686 if (l-k/=0) at=aa**(l-k) !This determine helps to improve efficient 2687 if (m+k-j/=0) bt=bb**(m+k-j) 2688 fj=fj+pre*at*bt 2689end do 2690end function 2691 2692!!!---- Calculate int('t^(2*m)*exp(-x*t^2)','t',0,1) see Cook book p281 for detail 2693real*8 function Fmch(m,x,expnx) 2694! expnx is input parameter, value should be exp(-x), because calculate the value is time-consuming and in 2695! other place this value also need be calculate, so not recalculate in this subroutine 2696IMPLICIT none 2697integer m,i 2698real*8 x,expnx,a,b,term,partsum,APPROX,xd,FIMULT,NOTRMS,eps,fiprop 2699eps=1.0D-8 !convergence precision 2700Fmch=0D0 2701if (x<=10) then 2702 if (expnx==0D0) RETURN 2703 A=m+0.5D0 2704 term=1.0D0/A 2705 partsum=term 2706 DO I=2,50 2707 A=A+1.0D0 2708 term=term*X/A 2709 partsum=partsum+term 2710 if ( term/partsum < eps) THEN 2711 Fmch = 0.5D0*partsum*expnx 2712 RETURN 2713 END IF 2714 END DO 2715 write(*,*) "Error: Fmch didn't converge" 2716else !x is big, use suitable method for solve this situation 2717 A=M 2718 B=A+0.5D0 2719 A=A-0.5D0 2720 XD=1.0D0/X 2721 APPROX=0.88622692D0*(dsqrt(XD)*XD**m) 2722 DO I=1,m 2723 B=B-1.0D0 2724 APPROX=APPROX*B 2725 END DO 2726 FIMULT=0.5D0*expnx*XD 2727 partsum=0.D0 2728 IF (FIMULT==0.0D0) THEN 2729 Fmch=APPROX-FIMULT*partsum 2730 return 2731 ELSE 2732 FIPROP=FIMULT/APPROX 2733 term=1.0d0 2734 partsum=term 2735 NOTRMS=X 2736 NOTRMS=NOTRMS+M 2737 DO I=2,NOTRMS 2738 term=term*A*XD 2739 partsum=partsum+term 2740 IF (dabs(term*FIPROP/partsum)<eps) THEN 2741 Fmch=APPROX-FIMULT*partsum 2742 RETURN 2743 END IF 2744 A=A-1.0D0 2745 END DO 2746 write(*,*) "Didn't converge" 2747 END IF 2748end if 2749end function 2750 2751 2752 2753 2754!!!------------- Generate Becke weight function, only used by Sobereva 2755!If inp2=0, then return atomic space weight of atom inp1, else return overlap weight of atom inp1 and inp2 2756real*8 function beckewei(x,y,z,inp1,inp2) 2757use defvar 2758real*8 x,y,z 2759integer inp1,inp2 2760! integer :: fraglist(13)=(/ 24,12,23,4,11,3,22,1,10,2,18,20,21 /) 2761real*8 smat(ncenter,ncenter),Pvec(ncenter) 2762smat=1.0D0 2763do ii=1,ncenter 2764 ri=dsqrt( (x-a(ii)%x)**2+(y-a(ii)%y)**2+(z-a(ii)%z)**2 ) 2765 do jj=1,ncenter 2766 if (ii==jj) cycle 2767 rj=dsqrt( (x-a(jj)%x)**2+(y-a(jj)%y)**2+(z-a(jj)%z)**2 ) 2768 rmiu=(ri-rj)/distmat(ii,jj) 2769 2770 chi=covr_tianlu(a(ii)%index)/covr_tianlu(a(jj)%index) !Adjust for heteronuclear 2771 uij=(chi-1)/(chi+1) 2772 aij=uij/(uij**2-1) 2773 if (aij>0.5D0) aij=0.5D0 2774 if (aij<-0.5D0) aij=-0.5D0 2775 rmiu=rmiu+aij*(1-rmiu**2) 2776 tmps=rmiu 2777 do iter=1,3 2778 tmps=1.5D0*(tmps)-0.5D0*(tmps)**3 2779 end do 2780 smat(ii,jj)=0.5D0*(1-tmps) 2781 end do 2782end do 2783 2784Pvec=1.0D0 2785do ii=1,ncenter 2786 Pvec=Pvec*smat(:,ii) 2787end do 2788Pvec=Pvec/sum(Pvec) 2789! beckewei=Pvec(2)*Pvec(3)*Pvec(4) 2790if (inp2==0) then 2791 beckewei=Pvec(inp1) 2792else 2793 beckewei=Pvec(inp1)*Pvec(inp2) 2794end if 2795! beckewei=sum(Pvec(fraglist)) !Get fragmental Becke weight 2796! beckewei=sum(Pvec(1:5)) 2797end function 2798 2799 2800!!!-------- Ellipticity of electron density (itype=1) or eta (itype=2) 2801real*8 function densellip(x,y,z,itype) 2802use util 2803integer itype 2804real*8 x,y,z,dens,grad(3),hess(3,3),eigval(3),eigvecmat(3,3) 2805call calchessmat_dens(2,x,y,z,dens,grad,hess) 2806call diagmat(hess,eigvecmat,eigval,300,1D-12) 2807call sortr8(eigval) 2808eigmax=eigval(3) 2809eigmed=eigval(2) 2810eigmin=eigval(1) 2811if (itype==1) then 2812 densellip=eigmin/eigmed-1 2813else 2814 densellip=abs(eigmin)/eigmax 2815end if 2816end function 2817 2818 2819!!!-------------- Single exponential decay detector (SEDD) 2820!Originally proposed in ChemPhysChem, 13, 3462 (2012), current implementation is based on the new definition in DORI paper (10.1021/ct500490b) 2821real*8 function SEDD(x,y,z) 2822real*8 x,y,z,rho,grad(3),hess(3,3) 2823call calchessmat_dens(2,x,y,z,rho,grad,hess) 2824dersqr=sum(grad**2) 2825tmp1_1=rho*(grad(1)*hess(1,1)+grad(2)*hess(1,2)+grad(3)*hess(1,3)) 2826tmp1_2=grad(1)*dersqr 2827tmp2_1=rho*(grad(1)*hess(1,2)+grad(2)*hess(2,2)+grad(3)*hess(2,3)) 2828tmp2_2=grad(2)*dersqr 2829tmp3_1=rho*(grad(1)*hess(1,3)+grad(2)*hess(2,3)+grad(3)*hess(3,3)) 2830tmp3_2=grad(3)*dersqr 2831eps=4/rho**8*( (tmp1_1-tmp1_2)**2 + (tmp2_1-tmp2_2)**2 + (tmp3_1-tmp3_2)**2 ) 2832SEDD=dlog(1+eps) 2833end function 2834 2835 2836!!!-------------- Density Overlap Regions Indicator (DORI) DOI:10.1021/ct500490b 2837real*8 function DORI(x,y,z) 2838real*8 x,y,z,rho,grad(3),hess(3,3) 2839call calchessmat_dens(2,x,y,z,rho,grad,hess) 2840dersqr=sum(grad**2) 2841tmp1_1=rho*(grad(1)*hess(1,1)+grad(2)*hess(1,2)+grad(3)*hess(1,3)) 2842tmp1_2=grad(1)*dersqr 2843tmp2_1=rho*(grad(1)*hess(1,2)+grad(2)*hess(2,2)+grad(3)*hess(2,3)) 2844tmp2_2=grad(2)*dersqr 2845tmp3_1=rho*(grad(1)*hess(1,3)+grad(2)*hess(2,3)+grad(3)*hess(3,3)) 2846tmp3_2=grad(3)*dersqr 2847theta=4/dersqr**3*( (tmp1_1-tmp1_2)**2 + (tmp2_1-tmp2_2)**2 + (tmp3_1-tmp3_2)**2 ) 2848DORI=theta/(1+theta) 2849end function 2850 2851 2852!!!----- Integrand of LSDA exchange functional 2853real*8 function xLSDA(x,y,z) 2854real*8 x,y,z 2855call gendensgradab(x,y,z,adens,bdens,tdens,agrad,bgrad,tgrad,abgrad) 2856xLSDA=-3D0/2D0*(3D0/4D0/pi)**(1D0/3D0)*(adens**(4D0/3D0)+bdens**(4D0/3D0)) 2857end function 2858 2859!!!----- Integrand of Becke88 exchange functional 2860real*8 function xBecke88(x,y,z) 2861real*8 x,y,z 2862call gendensgradab(x,y,z,adens,bdens,tdens,agrad,bgrad,tgrad,abgrad) 2863adens4d3=adens**(4D0/3D0) 2864bdens4d3=bdens**(4D0/3D0) 2865slatercoeff=-3D0/2D0*(3D0/4D0/pi)**(1D0/3D0) 2866slaterxa=slatercoeff*adens4d3 2867slaterxb=slatercoeff*bdens4d3 2868slaterx=slaterxa+slaterxb 2869redagrad=agrad/adens4d3 !Reduced density gradient 2870redbgrad=bgrad/bdens4d3 2871arshredagrad=log(redagrad+dsqrt(redagrad**2+1)) 2872Beckexa=adens4d3*redagrad**2/(1+6*0.0042D0*redagrad*arshredagrad) 2873arshredbgrad=log(redbgrad+dsqrt(redbgrad**2+1)) 2874Beckexb=bdens4d3*redbgrad**2/(1+6*0.0042D0*redbgrad*arshredbgrad) 2875Beckex=-0.0042D0*(Beckexa+Beckexb) 2876xBecke88=slaterx+Beckex 2877end function 2878 2879!!!----- Integrand of LYP corelation functional 2880real*8 function cLYP(x,y,z) 2881real*8 x,y,z 2882call gendensgradab(x,y,z,adens,bdens,tdens,agrad,bgrad,tgrad,abgrad) 2883parma=0.04918D0 2884parmb=0.132D0 2885parmc=0.2533D0 2886parmd=0.349D0 2887densn1d3=tdens**(-1D0/3D0) 2888parmw=exp(-parmc*densn1d3) / (1+parmd*densn1d3) * tdens**(-11D0/3D0) 2889parmdel=parmc*densn1d3+parmd*densn1d3/(1+parmd*densn1d3) 2890parmCf=3D0/10D0*(3*pi*pi)**(2D0/3D0) 2891tmp1=-parma*4D0/(1+parmd*densn1d3)*adens*bdens/tdens 2892tmp2a=2**(11D0/3D0)*parmCf*(adens**(8D0/3D0)+bdens**(8D0/3D0)) 2893tmp2b=(47D0/18D0-7D0/18D0*parmdel)*tgrad**2 2894tmp2c=-(2.5D0-parmdel/18D0)*(agrad**2+bgrad**2) 2895tmp2d=-(parmdel-11D0)/9D0*(adens/tdens*agrad**2+bdens/tdens*bgrad**2) 2896tmp2=adens*bdens*(tmp2a+tmp2b+tmp2c+tmp2d) 2897tmp3=-2D0/3D0*tdens**2*tgrad**2+(2D0/3D0*tdens**2-adens**2)*bgrad**2+(2D0/3D0*tdens**2-bdens**2)*agrad**2 2898cLYP=tmp1-parma*parmb*parmw*(tmp2+tmp3) 2899end function 2900 2901 2902!!!--- Integrand of Weizsacker (steric energy) 2903! DO NOT consider EDF, because the result outputted by quantum chemistry programs is always for only valence electrons! 2904real*8 function weizsacker(x,y,z) 2905real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo),gradrho(3) !,EDFgrad(3) 2906call orbderv(2,1,nmo,x,y,z,wfnval,wfnderv) 2907rho=0D0 2908gradrho=0D0 2909do i=1,nmo 2910 rho=rho+MOocc(i)*wfnval(i)**2 2911 gradrho(:)=gradrho(:)+MOocc(i)*wfnval(i)*wfnderv(:,i) 2912end do 2913gradrho=2*gradrho 2914! if (allocated(b_EDF)) then 2915! call EDFrho(2,x,y,z,EDFdens,EDFgrad) 2916! rho=rho+EDFdens 2917! gradrho=gradrho+EDFgrad 2918! end if 2919if (rho<1D-30) then 2920 weizsacker=0 2921else 2922 weizsacker=sum(gradrho(:)**2)/8/rho 2923end if 2924end function 2925!!!--- Weizsacker potential 2926real*8 function weizpot(x,y,z) 2927real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo),gradrho(3),laplx,laply,laplz,lapltot 2928rho=0D0 2929gradrho=0D0 2930laplx=0D0 2931laply=0D0 2932laplz=0D0 2933call orbderv(3,1,nmo,x,y,z,wfnval,wfnderv,wfnhess) 2934do i=1,nmo 2935 rho=rho+MOocc(i)*wfnval(i)**2 2936 gradrho(:)=gradrho(:)+MOocc(i)*wfnval(i)*wfnderv(:,i) 2937 laplx=laplx+MOocc(i)*( wfnderv(1,i)**2 + wfnval(i)*wfnhess(1,1,i) ) 2938 laply=laply+MOocc(i)*( wfnderv(2,i)**2 + wfnval(i)*wfnhess(2,2,i) ) 2939 laplz=laplz+MOocc(i)*( wfnderv(3,i)**2 + wfnval(i)*wfnhess(3,3,i) ) 2940end do 2941gradrho=2*gradrho 2942lapltot=2*(laplx+laply+laplz) 2943weizpot=sum(gradrho(:)**2)/8D0/rho**2-lapltot/4D0/rho 2944end function 2945 2946 2947!!!--- Steric potential (J. Chem. Phys., 126, 244103), which negative value is one-electron potential 2948real*8 function stericpot(x,y,z) 2949real*8 x,y,z,gradrho(3),hessrho(3,3),lapltot 2950call calchessmat_dens(2,x,y,z,rho,gradrho,hessrho) 2951lapltot=hessrho(1,1)+hessrho(2,2)+hessrho(3,3) 2952if (rho<steric_potcutrho) then 2953 stericpot=steric_potcons 2954 return 2955end if 2956stericpot=sum(gradrho**2)/(rho+steric_addminimal)**2/8D0-lapltot/(rho+steric_addminimal)/4D0 2957end function 2958 2959!!!----- Calculate the first-order derivative of steric potential 2960subroutine stericderv(x,y,z,derv) 2961real*8 x,y,z,derv(3) 2962real*8 eleval,elegrad(3),elehess(3,3),laplval,laplgrad(3),rhotens3(3,3,3) 2963real*8 wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo),wfntens3(3,3,3,nmo) 2964real*8 EDFval,EDFgrad(3),EDFhess(3,3),EDFtens3(3,3,3) 2965call orbderv(5,1,nmo,x,y,z,wfnval,wfnderv,wfnhess,wfntens3) 2966eleval=sum( MOocc(1:nmo)*wfnval(1:nmo)*wfnval(1:nmo) ) 2967elegrad(1)=2*sum( MOocc(1:nmo)*wfnval(1:nmo)*wfnderv(1,1:nmo) ) 2968elegrad(2)=2*sum( MOocc(1:nmo)*wfnval(1:nmo)*wfnderv(2,1:nmo) ) 2969elegrad(3)=2*sum( MOocc(1:nmo)*wfnval(1:nmo)*wfnderv(3,1:nmo) ) 2970elehess(1,1)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)**2 + wfnval(1:nmo)*wfnhess(1,1,1:nmo) ) ) 2971elehess(2,2)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)**2 + wfnval(1:nmo)*wfnhess(2,2,1:nmo) ) ) 2972elehess(3,3)=2*sum( MOocc(1:nmo)*( wfnderv(3,1:nmo)**2 + wfnval(1:nmo)*wfnhess(3,3,1:nmo) ) ) 2973elehess(1,2)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(2,1:nmo)+wfnhess(1,2,1:nmo)*wfnval(1:nmo) ) ) 2974elehess(2,3)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)*wfnderv(3,1:nmo)+wfnhess(2,3,1:nmo)*wfnval(1:nmo) ) ) 2975elehess(1,3)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(3,1:nmo)+wfnhess(1,3,1:nmo)*wfnval(1:nmo) ) ) 2976elehess(2,1)=elehess(1,2) 2977elehess(3,2)=elehess(2,3) 2978elehess(3,1)=elehess(1,3) 2979laplval=elehess(1,1)+elehess(2,2)+elehess(3,3) 2980rhotens3=0D0 2981do i=1,nmo 2982 rhotens3(1,1,1)=rhotens3(1,1,1)+MOocc(i)*( 3*wfnderv(1,i)*wfnhess(1,1,i)+wfnval(i)*wfntens3(1,1,1,i) ) 2983 rhotens3(2,2,2)=rhotens3(2,2,2)+MOocc(i)*( 3*wfnderv(2,i)*wfnhess(2,2,i)+wfnval(i)*wfntens3(2,2,2,i) ) 2984 rhotens3(3,3,3)=rhotens3(3,3,3)+MOocc(i)*( 3*wfnderv(3,i)*wfnhess(3,3,i)+wfnval(i)*wfntens3(3,3,3,i) ) 2985 rhotens3(1,1,2)=rhotens3(1,1,2)+MOocc(i)*( 2*wfnderv(1,i)*wfnhess(1,2,i)+wfnderv(2,i)*wfnhess(1,1,i)+wfnval(i)*wfntens3(1,1,2,i) ) 2986 rhotens3(1,1,3)=rhotens3(1,1,3)+MOocc(i)*( 2*wfnderv(1,i)*wfnhess(1,3,i)+wfnderv(3,i)*wfnhess(1,1,i)+wfnval(i)*wfntens3(1,1,3,i) ) 2987 rhotens3(2,2,3)=rhotens3(2,2,3)+MOocc(i)*( 2*wfnderv(2,i)*wfnhess(2,3,i)+wfnderv(3,i)*wfnhess(2,2,i)+wfnval(i)*wfntens3(2,2,3,i) ) 2988 rhotens3(1,2,2)=rhotens3(1,2,2)+MOocc(i)*( 2*wfnderv(2,i)*wfnhess(2,1,i)+wfnderv(1,i)*wfnhess(2,2,i)+wfnval(i)*wfntens3(2,2,1,i) ) 2989 rhotens3(1,3,3)=rhotens3(1,3,3)+MOocc(i)*( 2*wfnderv(3,i)*wfnhess(3,1,i)+wfnderv(1,i)*wfnhess(3,3,i)+wfnval(i)*wfntens3(3,3,1,i) ) 2990 rhotens3(2,3,3)=rhotens3(2,3,3)+MOocc(i)*( 2*wfnderv(3,i)*wfnhess(3,2,i)+wfnderv(2,i)*wfnhess(3,3,i)+wfnval(i)*wfntens3(3,3,2,i) ) 2991end do 2992rhotens3=rhotens3*2D0 2993laplgrad(1)=rhotens3(1,1,1)+rhotens3(1,2,2)+rhotens3(1,3,3) 2994laplgrad(2)=rhotens3(1,1,2)+rhotens3(2,2,2)+rhotens3(2,3,3) 2995laplgrad(3)=rhotens3(1,1,3)+rhotens3(2,2,3)+rhotens3(3,3,3) 2996if (allocated(b_EDF)) then 2997 call EDFrho(5,x,y,z,EDFval,EDFgrad,EDFhess,EDFtens3) 2998 eleval=eleval+EDFval 2999 elegrad=elegrad+EDFgrad 3000 elehess=elehess+EDFhess 3001 laplgrad(1)=laplgrad(1)+EDFtens3(1,1,1)+EDFtens3(1,2,2)+EDFtens3(1,3,3) 3002 laplgrad(2)=laplgrad(2)+EDFtens3(1,1,2)+EDFtens3(2,2,2)+EDFtens3(2,3,3) 3003 laplgrad(3)=laplgrad(3)+EDFtens3(1,1,3)+EDFtens3(2,2,3)+EDFtens3(3,3,3) 3004end if 3005! Above codes can be simplified as below two lines, but will consume additional 1/3 time 3006! call calchessmat_dens(2,x,y,z,eleval,elegrad,elehess) 3007! call calchessmat_lapl(1,x,y,z,laplval,laplgrad,laplhess) 3008 3009eleval=eleval+steric_addminimal 3010elenorm2=sum(elegrad**2) !Square of norm of electron density gradient 3011do i=1,3 !x,y,z 3012 tmp1=(elegrad(1)*elehess(1,i)+elegrad(2)*elehess(2,i)+elegrad(3)*elehess(3,i)) / eleval**2 -elenorm2/eleval**3*elegrad(i) 3013 tmp2=-laplgrad(i)/eleval+laplval/eleval**2*elegrad(i) 3014 derv(i)=(tmp1+tmp2)/4D0 3015end do 3016end subroutine 3017 3018!---- Magnitude of steric force 3019real*8 function stericforce(x,y,z) 3020real*8 x,y,z,derv(3) 3021call stericderv(x,y,z,derv) 3022stericforce=dsqrt(sum(derv**2)) 3023end function 3024 3025 3026!memo: Shubin reported that steric potential/force is quite sensitive to steric_addminimal, & 3027!so 2016-Sep-30 I devised another solution to solve the diverse behavior of steric potential/force via Becke's damping function 3028!!!--- Steric potential with damping function to zero 3029real*8 function stericpot_damp(x,y,z) 3030real*8 x,y,z,gradrho(3),hessrho(3,3),lapltot 3031call calchessmat_dens(2,x,y,z,rho,gradrho,hessrho) 3032lapltot=hessrho(1,1)+hessrho(2,2)+hessrho(3,3) 3033stericpotorg=sum(gradrho**2)/rho**2/8D0-lapltot/rho/4D0 3034weiwidth=2 !e.g. weiwidth=2 and steric_potcutrho=-13 means the Becke damping function of [1,0] corresponds to 1D-11~1D-15 3035tmps=-(dlog(rho)-steric_potcutrho)/weiwidth 3036if (tmps<-1) then 3037 consorg=1 3038else if (tmps>1) then 3039 consorg=0 3040else 3041 do iter=1,2 3042 tmps=1.5D0*(tmps)-0.5D0*(tmps)**3 3043 end do 3044 consorg=0.5D0*(1-tmps) !The weight to switch to constant value steric_potcons 3045end if 3046stericpot_damp=stericpotorg*consorg+steric_potcons*(1-consorg) 3047end function 3048!!!--- Steric force based on damped steric potential 3049real*8 function stericforce_damp(x,y,z) 3050real*8 x,y,z,derv(3) 3051diffstep=2D-5 3052derv(1)=(stericpot_damp(x+diffstep,y,z)-stericpot_damp(x-diffstep,y,z))/(2*diffstep) 3053derv(2)=(stericpot_damp(x,y+diffstep,z)-stericpot_damp(x,y-diffstep,z))/(2*diffstep) 3054derv(3)=(stericpot_damp(x,y,z+diffstep)-stericpot_damp(x,y,z-diffstep))/(2*diffstep) 3055stericforce_damp=dsqrt(sum(derv**2)) 3056end function 3057!!!--- Steric force directly damped to zero rather than based on damped steric potential 3058real*8 function stericforce_directdamp(x,y,z) 3059real*8 x,y,z 3060weiwidth=2 3061tmps=-(dlog(fdens(x,y,z))-steric_potcutrho)/weiwidth 3062if (tmps<-1) then 3063 consorg=1 3064else if (tmps>1) then 3065 consorg=0 3066else 3067 do iter=1,2 3068 tmps=1.5D0*(tmps)-0.5D0*(tmps)**3 3069 end do 3070 consorg=0.5D0*(1-tmps) 3071end if 3072steric_addminimalold=steric_addminimal 3073steric_addminimal=0 3074stericforce_directdamp=stericforce(x,y,z)*consorg 3075steric_addminimal=steric_addminimalold 3076end function 3077 3078 3079 3080!!!----- Steric charge, =lapl(steric potential)/(-4*pi) 3081! Based on analytical first-order derivative, using finite difference to obtain d2v/dx2, d2v/dy2 and d2v/dz2 3082real*8 function stericcharge(x,y,z) 3083real*8 x,y,z,derv1add(3),derv1min(3) 3084if (fdens(x,y,z)<steric_potcutrho) then 3085 stericcharge=0D0 3086 return 3087end if 3088diffstep=2D-5 3089call stericderv(x+diffstep,y,z,derv1add) 3090call stericderv(x-diffstep,y,z,derv1min) 3091derv2x=(derv1add(1)-derv1min(1))/(2*diffstep) !d2v/dx2 3092call stericderv(x,y+diffstep,z,derv1add) 3093call stericderv(x,y-diffstep,z,derv1min) 3094derv2y=(derv1add(2)-derv1min(2))/(2*diffstep) !d2v/dy2 3095call stericderv(x,y,z+diffstep,derv1add) 3096call stericderv(x,y,z-diffstep,derv1min) 3097derv2z=(derv1add(3)-derv1min(3))/(2*diffstep) !d2v/dz2 3098stericcharge=-(derv2x+derv2y+derv2z)/4D0/pi 3099end function 3100 3101 3102!!!----- Calculate Fisher information density, see JCP,126,191107 for example 3103!itype=1 Normal definition 3104!itype=2 Second Fisher information density, Eq.5 of JCP,126,191107 3105real*8 function Fisherinfo(itype,x,y,z) 3106real*8 x,y,z,eleval,elegrad(3),elehess(3,3) 3107integer itype 3108Fisherinfo=0 3109eleval=fdens(x,y,z) 3110if (eleval<=1D-30) return 3111if (itype==1) Fisherinfo=fgrad(x,y,z,'t')**2/eleval 3112if (itype==2) Fisherinfo=-flapl(x,y,z,'t')*log(eleval) 3113end function 3114 3115 3116!!!----- Ghosh entropy density, PNAS, 81, 8028 3117!If itype==1, G(r) will be used as kinetic energy density 3118!If itype==2, G(r)-der2rho/8 will be used instead, which is the kinetic energy density exactly corresponding to Eq. 22 of PNAS, 81, 8028. 3119real*8 function Ghoshentro(x,y,z,itype) 3120integer itype 3121real*8 kintot,x,y,z,wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo) 3122if (itype==1) call orbderv(2,1,nmo,x,y,z,wfnval,wfnderv) 3123if (itype==2) call orbderv(3,1,nmo,x,y,z,wfnval,wfnderv,wfnhess) !If K(r) is used, use this!!! 3124rho=0D0 3125do i=1,nmo 3126 rho=rho+MOocc(i)*wfnval(i)**2 3127end do 3128ck=2.871234D0 3129TFkin=ck*rho**(5D0/3D0) 3130kintot=0D0 3131! If we use Hamiltonian kinetic density 3132! hamx=sum( MOocc(1:nmo)*wfnhess(1,1,1:nmo)*wfnval(1:nmo) ) 3133! hamy=sum( MOocc(1:nmo)*wfnhess(2,2,1:nmo)*wfnval(1:nmo) ) 3134! hamz=sum( MOocc(1:nmo)*wfnhess(3,3,1:nmo)*wfnval(1:nmo) ) 3135! kintot=-(hamx+hamy+hamz)/2 3136! If we use Lagrangian kinetic density G(r) 3137do i=1,nmo 3138 kintot=kintot+MOocc(i)*sum(wfnderv(:,i)**2) 3139end do 3140kintot=kintot/2D0 3141if (itype==2) then 3142 xlapl=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)**2 + wfnval(1:nmo)*wfnhess(1,1,1:nmo) ) ) 3143 ylapl=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)**2 + wfnval(1:nmo)*wfnhess(2,2,1:nmo) ) ) 3144 zlapl=2*sum( MOocc(1:nmo)*( wfnderv(3,1:nmo)**2 + wfnval(1:nmo)*wfnhess(3,3,1:nmo) ) ) 3145 kintot=kintot-(xlapl+ylapl+zlapl)/8 3146end if 3147! if (kintot<0) then 3148! write(*,"(5f16.10)") kintot,TFkin,x,y,z 3149! read(*,*) 3150! end if 3151rlambda=5D0/3D0+log(4D0*pi*ck/3D0) 3152if (kintot<0) then 3153 rlogterm=0 3154else 3155 rlogterm=log(kintot/TFkin) 3156end if 3157Ghoshentro=1.5D0*rho*(rlambda+rlogterm) 3158end function 3159 3160 3161 3162!!!----- Pauli potential, see Shubin's paper: Comp. Theor. Chem., 1006, 92-99 3163!Only suitable for close-shell cases 3164real*8 function paulipot(x,y,z) 3165real*8 x,y,z 3166paulipot=totesp(x,y,z)-DFTxcpot(x,y,z)-weizpot(x,y,z) !Note that the sign of ESP in shubin's CTC paper is inversed 3167end function 3168 3169!!!----- The magnitude of Pauli force, namely the gradient norm of negative Pauli potential 3170real*8 function pauliforce(x,y,z) 3171real*8 x,y,z 3172diff=2D-5 3173forcex=-(paulipot(x+diff,y,z)-paulipot(x-diff,y,z))/(2*diff) 3174forcey=-(paulipot(x,y+diff,z)-paulipot(x,y-diff,z))/(2*diff) 3175forcez=-(paulipot(x,y,z+diff)-paulipot(x,y,z-diff))/(2*diff) 3176pauliforce=dsqrt(forcex**2+forcey**2+forcez**2) 3177end function 3178 3179!!!----- Pauli charge, =lapl(Pauli potential)/(-4*pi) 3180real*8 function paulicharge(x,y,z) 3181real*8 x,y,z 3182diff=2D-4 !Should not be smaller, otherwise some dirty points will be presented 3183value=paulipot(x,y,z) 3184valuexaddadd=paulipot(x+2*diff,y,z) 3185valuexminmin=paulipot(x-2*diff,y,z) 3186valueyaddadd=paulipot(x,y+2*diff,z) 3187valueyminmin=paulipot(x,y-2*diff,z) 3188valuezaddadd=paulipot(x,y,z+2*diff) 3189valuezminmin=paulipot(x,y,z-2*diff) 3190xcomp=(valuexaddadd-2*value+valuexminmin)/(2*diff)**2 3191ycomp=(valueyaddadd-2*value+valueyminmin)/(2*diff)**2 3192zcomp=(valuezaddadd-2*value+valuezminmin)/(2*diff)**2 3193paulicharge=(xcomp+ycomp+zcomp)/(-4*pi) 3194end function 3195 3196!!!----- Quantum potential 3197real*8 function quantumpot(x,y,z) 3198real*8 x,y,z 3199quantumpot=totesp(x,y,z)-weizpot(x,y,z) 3200end function 3201 3202!!!----- The magnitude of quantum force, namely the gradient norm of quantum potential 3203real*8 function quantumforce(x,y,z) 3204real*8 x,y,z 3205diff=2D-5 3206forcex=-(quantumpot(x+diff,y,z)-quantumpot(x-diff,y,z))/(2*diff) 3207forcey=-(quantumpot(x,y+diff,z)-quantumpot(x,y-diff,z))/(2*diff) 3208forcez=-(quantumpot(x,y,z+diff)-quantumpot(x,y,z-diff))/(2*diff) 3209quantumforce=dsqrt(forcex**2+forcey**2+forcez**2) 3210end function 3211 3212!!!----- Quantum charge 3213real*8 function quantumcharge(x,y,z) 3214real*8 x,y,z 3215diff=2D-4 !Should not be smaller, otherwise some dirty points will be presented 3216value=quantumpot(x,y,z) 3217valuexaddadd=quantumpot(x+2*diff,y,z) 3218valuexminmin=quantumpot(x-2*diff,y,z) 3219valueyaddadd=quantumpot(x,y+2*diff,z) 3220valueyminmin=quantumpot(x,y-2*diff,z) 3221valuezaddadd=quantumpot(x,y,z+2*diff) 3222valuezminmin=quantumpot(x,y,z-2*diff) 3223xcomp=(valuexaddadd-2*value+valuexminmin)/(2*diff)**2 3224ycomp=(valueyaddadd-2*value+valueyminmin)/(2*diff)**2 3225zcomp=(valuezaddadd-2*value+valuezminmin)/(2*diff)**2 3226quantumcharge=(xcomp+ycomp+zcomp)/(-4*pi) 3227end function 3228 3229!!!----- The magnitude of electrostatic force, namely the gradient norm of electrostatic potential 3230real*8 function elestatforce(x,y,z) 3231real*8 x,y,z 3232diff=2D-5 3233forcex=-(totesp(x+diff,y,z)-totesp(x-diff,y,z))/(2*diff) 3234forcey=-(totesp(x,y+diff,z)-totesp(x,y-diff,z))/(2*diff) 3235forcez=-(totesp(x,y,z+diff)-totesp(x,y,z-diff))/(2*diff) 3236elestatforce=dsqrt(forcex**2+forcey**2+forcez**2) 3237end function 3238 3239!!!----- Electrostatic charge 3240real*8 function elestatcharge(x,y,z) 3241real*8 x,y,z 3242diff=2D-4 !Should not be smaller, otherwise some dirty points will be presented 3243value=totesp(x,y,z) 3244valuexaddadd=totesp(x+2*diff,y,z) 3245valuexminmin=totesp(x-2*diff,y,z) 3246valueyaddadd=totesp(x,y+2*diff,z) 3247valueyminmin=totesp(x,y-2*diff,z) 3248valuezaddadd=totesp(x,y,z+2*diff) 3249valuezminmin=totesp(x,y,z-2*diff) 3250xcomp=(valuexaddadd-2*value+valuexminmin)/(2*diff)**2 3251ycomp=(valueyaddadd-2*value+valueyminmin)/(2*diff)**2 3252zcomp=(valuezaddadd-2*value+valuezminmin)/(2*diff)**2 3253elestatcharge=(xcomp+ycomp+zcomp)/(-4*pi) 3254end function 3255 3256 3257 3258!!!---- Use trilinear interpolation to obtain value at a given point by using cubmat 3259!itype==1: interpolate from cubmat, =2: from cubmattmp 3260real*8 function linintp3d(x,y,z,itype) 3261real*8 x,y,z 3262integer itype 3263character*80 c80tmp 3264do ix=1,nx 3265 x1=orgx+(ix-1)*dx 3266 x2=orgx+ix*dx 3267 if (x>=x1.and.x<x2) exit !1D-10 is used to avoid numerical uncertainty 3268end do 3269do iy=1,ny 3270 y1=orgy+(iy-1)*dy 3271 y2=orgy+iy*dy 3272 if (y>=y1.and.y<y2) exit 3273end do 3274do iz=1,nz 3275 z1=orgz+(iz-1)*dz 3276 z2=orgz+iz*dz 3277 if (z>=z1.and.z<z2) exit 3278end do 3279if (ix>=nx.or.iy>=ny.or.iz>=nz) then !Out of grid data range 3280 linintp3d=0D0 3281else 3282 if (itype==1) then 3283 valxy1=( cubmat(ix,iy,iz )*(x2-x)*(y2-y) + cubmat(ix+1,iy,iz )*(x-x1)*(y2-y) + & 3284 cubmat(ix,iy+1,iz )*(x2-x)*(y-y1) + cubmat(ix+1,iy+1,iz )*(x-x1)*(y-y1) ) /dx/dy 3285 valxy2=( cubmat(ix,iy,iz+1)*(x2-x)*(y2-y) + cubmat(ix+1,iy,iz+1)*(x-x1)*(y2-y) + & 3286 cubmat(ix,iy+1,iz+1)*(x2-x)*(y-y1) + cubmat(ix+1,iy+1,iz+1)*(x-x1)*(y-y1) ) /dx/dy 3287 else 3288 valxy1=( cubmattmp(ix,iy,iz )*(x2-x)*(y2-y) + cubmattmp(ix+1,iy,iz )*(x-x1)*(y2-y) + & 3289 cubmattmp(ix,iy+1,iz )*(x2-x)*(y-y1) + cubmattmp(ix+1,iy+1,iz )*(x-x1)*(y-y1) ) /dx/dy 3290 valxy2=( cubmattmp(ix,iy,iz+1)*(x2-x)*(y2-y) + cubmattmp(ix+1,iy,iz+1)*(x-x1)*(y2-y) + & 3291 cubmattmp(ix,iy+1,iz+1)*(x2-x)*(y-y1) + cubmattmp(ix+1,iy+1,iz+1)*(x-x1)*(y-y1) ) /dx/dy 3292 end if 3293 linintp3d=valxy1+(z-z1)*(valxy2-valxy1)/dz 3294end if 3295end function 3296 3297 3298!!!---- Trilinear interpolation of 3D-vector field by using cubmatvec 3299subroutine linintp3dvec(x,y,z,vecintp) 3300real*8 x,y,z,vecintp(3),valxy1(3),valxy2(3) 3301do ix=1,nx 3302 x1=orgx+(ix-1)*dx 3303 x2=orgx+ix*dx 3304 if (x>=x1.and.x<x2) exit 3305end do 3306do iy=1,ny 3307 y1=orgy+(iy-1)*dy 3308 y2=orgy+iy*dy 3309 if (y>=y1.and.y<y2) exit 3310end do 3311do iz=1,nz 3312 z1=orgz+(iz-1)*dz 3313 z1=orgz+iz*dz 3314 if (z>=z1.and.z<z2) exit 3315end do 3316if (ix>nx.or.iy>ny.or.iz>nz) then !Out of grid data range 3317 vecintp=0D0 3318else 3319 valxy1(:)=( cubmatvec(:,ix,iy,iz )*(x2-x)*(y2-y) + cubmatvec(:,ix+1,iy,iz )*(x-x1)*(y2-y) + & 3320 cubmatvec(:,ix,iy+1,iz )*(x2-x)*(y-y1) + cubmatvec(:,ix+1,iy+1,iz )*(x-x1)*(y-y1) ) /dx/dy 3321 valxy2(:)=( cubmatvec(:,ix,iy,iz+1)*(x2-x)*(y2-y) + cubmatvec(:,ix+1,iy,iz+1)*(x-x1)*(y2-y) + & 3322 cubmatvec(:,ix,iy+1,iz+1)*(x2-x)*(y-y1) + cubmatvec(:,ix+1,iy+1,iz+1)*(x-x1)*(y-y1) ) /dx/dy 3323 vecintp=valxy1+(z-z1)*(valxy2-valxy1)/dz 3324end if 3325end subroutine 3326 3327 3328 3329 3330!!!---- Calculate various kinds of integrand of DFT exchange-correlation functionals 3331!The routines are provided by DFT repository (ftp://ftp.dl.ac.uk/qcg/dft_library/index.html) 3332!The global variable "iDFTxcsel" is used to select the XC functional, see manual 3333!Note that the inner core density represented by EDF field is not taken into account 3334real*8 function DFTxcfunc(x,y,z) 3335real*8 x,y,z 3336if (wfntype==0.or.wfntype==3) then !Close-shell 3337 DFTxcfunc=DFTxcfunc_close(x,y,z) 3338else !Open-shell 3339 DFTxcfunc=DFTxcfunc_open(x,y,z) 3340end if 3341end function 3342!---- Calculate various kinds of DFT exchange-correlation potentials, see the comment of DFTxcfunc 3343real*8 function DFTxcpot(x,y,z) 3344real*8 x,y,z 3345if (wfntype==0.or.wfntype==3) then !Close-shell 3346 DFTxcpot=DFTxcpot_close(x,y,z) 3347else !Open-shell 3348 DFTxcpot=DFTxcpot_open(x,y,z) 3349 write(*,*) "XC potential for open-shell has not been supported yet!" 3350 read(*,*) 3351end if 3352end function 3353 3354 3355 3356!!---Close-shell form of DFTxcfunc routine 3357real*8 function DFTxcfunc_close(x,y,z) 3358real*8 x,y,z 3359call gendensgradab(x,y,z,adens,bdens,tdens,agrad,bgrad,tgrad,abgrad) 3360call getXCdata_close(0,tdens,tgrad**2,DFTxcfunc_close,rnouse,rnouse,rnouse,rnouse,rnouse) 3361end function 3362 3363!!---Close-shell form of DFTxcpot routine 3364real*8 function DFTxcpot_close(x,y,z) 3365real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo),wfnhess(3,3,nmo),gradrho(3),hessrho(3,3),tmparr(3,1),tmpval(1,1),lapltot 3366rho=0D0 3367gradrho=0D0 3368call orbderv(4,1,nmo,x,y,z,wfnval,wfnderv,wfnhess) 3369do i=1,nmo 3370 rho=rho+MOocc(i)*wfnval(i)**2 3371 gradrho(:)=gradrho(:)+MOocc(i)*wfnval(i)*wfnderv(:,i) 3372end do 3373gradrho=2*gradrho 3374hessrho(1,1)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)**2 + wfnval(1:nmo)*wfnhess(1,1,1:nmo) ) ) 3375hessrho(2,2)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)**2 + wfnval(1:nmo)*wfnhess(2,2,1:nmo) ) ) 3376hessrho(3,3)=2*sum( MOocc(1:nmo)*( wfnderv(3,1:nmo)**2 + wfnval(1:nmo)*wfnhess(3,3,1:nmo) ) ) 3377hessrho(1,2)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(2,1:nmo)+wfnhess(1,2,1:nmo)*wfnval(1:nmo) ) ) 3378hessrho(2,3)=2*sum( MOocc(1:nmo)*( wfnderv(2,1:nmo)*wfnderv(3,1:nmo)+wfnhess(2,3,1:nmo)*wfnval(1:nmo) ) ) 3379hessrho(1,3)=2*sum( MOocc(1:nmo)*( wfnderv(1,1:nmo)*wfnderv(3,1:nmo)+wfnhess(1,3,1:nmo)*wfnval(1:nmo) ) ) 3380hessrho(2,1)=hessrho(1,2) 3381hessrho(3,2)=hessrho(2,3) 3382hessrho(3,1)=hessrho(1,3) 3383lapltot=hessrho(1,1)+hessrho(2,2)+hessrho(3,3) 3384sigma=sum(gradrho(:)**2) 3385tmparr(:,1)=gradrho 3386tmpval=matmul(2*matmul(transpose(tmparr),hessrho),tmparr) !dot product between grad(sigma) and grad(rho) 3387call getXCdata_close(2,rho,sigma,value,d1rho,d1sig,d2rho,d2rhosig,d2sig) 3388DFTxcpot_close=d1rho-2*(d2rhosig*sigma+d2sig*tmpval(1,1)+d1sig*lapltot) 3389end function 3390 3391 3392 3393!!!For close-shell cases. Input rho and gradrho^2, return the value and its derivative of selected XC (or X/C only) functional 3394!The global variable "iDFTxcsel" is used to select the XC functional, see manual 3395!ixcderv=0: only get value, =1: also get d1rho and d1sig, =2: also get d2rho, d2rhosig and d2sig 3396!rho, sigma: inputted rho and gradrho^2 3397!value: outputted integrand of functional 3398!d1rho, d1sig: 1st derivative of functional w.r.t. rho and sigma, respectively 3399!d2rho, d2sig: 2nd derivative of functional w.r.t. rho and sigma, respectively 3400!d2rhosig: 1st derv w.r.t. rho and 1st derv w.r.t. sigma 3401subroutine getXCdata_close(ixcderv,rho,sigma,value,d1rho,d1sig,d2rho,d2rhosig,d2sig) 3402integer ixcderv 3403real*8 rho,sigma,value,d1rho,d1sig,d2rho,d2rhosig,d2sig,rhoa1(1),sigmaaa1(1) 3404real*8 XCzk(1),Xzk(1),Czk(1),XCvrhoa(1),Xvrhoa(1),Cvrhoa(1),XCvsigmaaa(1),Xvsigmaaa(1),Cvsigmaaa(1) 3405real*8 XCv2rhoa2(1),Xv2rhoa2(1),Cv2rhoa2(1),XCv2rhoasigmaaa(1),Xv2rhoasigmaaa(1),Cv2rhoasigmaaa(1) 3406real*8 XCv2sigmaaa2(1),Xv2sigmaaa2(1),Cv2sigmaaa2(1) 3407rhoa1(1)=rho 3408sigmaaa1(1)=sigma 3409!X part 3410if (iDFTxcsel==0.or.iDFTxcsel==80) then 3411 call rks_x_lda(ixcderv,1,rhoa1,sigmaaa1,Xzk,Xvrhoa,Xvsigmaaa,Xv2rhoa2,Xv2rhoasigmaaa,Xv2sigmaaa2) 3412else if (iDFTxcsel==1.or.iDFTxcsel==81.or.iDFTxcsel==82.or.iDFTxcsel==83) then 3413 call rks_x_b88(ixcderv,1,rhoa1,sigmaaa1,Xzk,Xvrhoa,Xvsigmaaa,Xv2rhoa2,Xv2rhoasigmaaa,Xv2sigmaaa2) 3414else if (iDFTxcsel==2.or.iDFTxcsel==84) then 3415 call rks_x_pbe(ixcderv,1,rhoa1,sigmaaa1,Xzk,Xvrhoa,Xvsigmaaa,Xv2rhoa2,Xv2rhoasigmaaa,Xv2sigmaaa2) 3416else if (iDFTxcsel==3.or.iDFTxcsel==85) then 3417 call rks_x_pw91(ixcderv,1,rhoa1,sigmaaa1,Xzk,Xvrhoa,Xvsigmaaa,Xv2rhoa2,Xv2rhoasigmaaa,Xv2sigmaaa2) 3418end if 3419!C part 3420if (iDFTxcsel==30.or.iDFTxcsel==80) then 3421 call rks_c_vwn5(ixcderv,1,rhoa1,sigmaaa1,Czk,Cvrhoa,Cvsigmaaa,Cv2rhoa2,Cv2rhoasigmaaa,Cv2sigmaaa2) 3422else if (iDFTxcsel==31.or.iDFTxcsel==81) then 3423 call rks_c_p86(ixcderv,1,rhoa1,sigmaaa1,Czk,Cvrhoa,Cvsigmaaa,Cv2rhoa2,Cv2rhoasigmaaa,Cv2sigmaaa2) 3424else if (iDFTxcsel==32.or.iDFTxcsel==82) then 3425 call rks_c_lyp(ixcderv,1,rhoa1,sigmaaa1,Czk,Cvrhoa,Cvsigmaaa,Cv2rhoa2,Cv2rhoasigmaaa,Cv2sigmaaa2) 3426else if (iDFTxcsel==33.or.iDFTxcsel==83.or.iDFTxcsel==85) then 3427 call rks_c_pw91(ixcderv,1,rhoa1,sigmaaa1,Czk,Cvrhoa,Cvsigmaaa,Cv2rhoa2,Cv2rhoasigmaaa,Cv2sigmaaa2) 3428else if (iDFTxcsel==34.or.iDFTxcsel==84) then 3429 call rks_c_pbe(ixcderv,1,rhoa1,sigmaaa1,Czk,Cvrhoa,Cvsigmaaa,Cv2rhoa2,Cv2rhoasigmaaa,Cv2sigmaaa2) 3430end if 3431!Whole XC 3432if (iDFTxcsel==70) then 3433 call rks_xc_b97(ixcderv,1,rhoa1,sigmaaa1,XCzk,XCvrhoa,XCvsigmaaa,XCv2rhoa2,XCv2rhoasigmaaa,XCv2sigmaaa2) 3434else if (iDFTxcsel==71) then 3435 call rks_xc_hcth407(ixcderv,1,rhoa1,sigmaaa1,XCzk,XCvrhoa,XCvsigmaaa,XCv2rhoa2,XCv2rhoasigmaaa,XCv2sigmaaa2) 3436else if (iDFTxcsel>=80.and.iDFTxcsel<99) then 3437 XCzk=Xzk+Czk 3438 XCvrhoa=Xvrhoa+Cvrhoa 3439 XCvsigmaaa=Xvsigmaaa+Cvsigmaaa 3440 XCv2rhoa2=Xv2rhoa2+Cv2rhoa2 3441 XCv2rhoasigmaaa=Xv2rhoasigmaaa+Cv2rhoasigmaaa 3442 XCv2sigmaaa2=Xv2sigmaaa2+Cv2sigmaaa2 3443end if 3444!Note that the problem of derivative of the DFT repository is revised by dividing a factor, similarly hereinafter 3445if (iDFTxcsel<30) then 3446 value=Xzk(1) 3447 d1rho=Xvrhoa(1) 3448 d1sig=Xvsigmaaa(1)/4D0 3449 d2rho=Xv2rhoa2(1)/2D0 3450 d2rhosig=Xv2rhoasigmaaa(1)/4D0 3451 d2sig=Xv2sigmaaa2(1)/16D0 3452else if (iDFTxcsel<70) then 3453 value=Czk(1) 3454 d1rho=Cvrhoa(1) 3455 d1sig=Cvsigmaaa(1)/4D0 3456 d2rho=Cv2rhoa2(1)/2D0 3457 d2rhosig=Cv2rhoasigmaaa(1)/4D0 3458 d2sig=Cv2sigmaaa2(1)/16D0 3459else if (iDFTxcsel<100) then 3460 value=XCzk(1) 3461 d1rho=XCvrhoa(1) 3462 d1sig=XCvsigmaaa(1)/4D0 3463 d2rho=XCv2rhoa2(1)/2D0 3464 d2rhosig=XCv2rhoasigmaaa(1)/4D0 3465 d2sig=XCv2sigmaaa2(1)/16D0 3466end if 3467end subroutine 3468 3469 3470 3471 3472 3473 3474!!---Open-shell form of DFTxcfunc routine 3475real*8 function DFTxcfunc_open(x,y,z) 3476real*8 x,y,z 3477call gendensgradab(x,y,z,adens,bdens,tdens,agrad,bgrad,tgrad,abgrad) 3478call getXCdata_open(0,adens,bdens,agrad**2,bgrad**2,abgrad,DFTxcfunc_open,& 3479sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb,sb) 3480end function 3481 3482!---Open-shell form of DFTxcpot routine 3483real*8 function DFTxcpot_open(x,y,z) 3484real*8 x,y,z 3485DFTxcpot_open=0D0 3486end function 3487 3488!!!For open-shell cases. Input rho and gradrho^2, return the value and its derivative of selected XC (or X/C only) functional 3489!The global variable "iDFTxcsel" is used to select the XC functional, see manual 3490!ixcderv=0: only get value, =1: also get 1st derv., =2: also get 2nd derv. 3491!Input quantities: 3492!rhoa=rho_alpha rhob=rho_beta 3493!sigaa=|gradrho_alpha|^2 sigbb=|gradrho_beta|^2 3494!sigab=Dot product between gradrho_alpha and gradrho_beta 3495subroutine getXCdata_open(ixcderv,rhoa,rhob,sigaa,sigbb,sigab,value,d1rhoa,d1rhob,d1sigaa,d1sigbb,d1sigab,d2rhoaa,d2rhobb,d2rhoab,& 3496d2rhoasigaa,d2rhoasigab,d2rhoasigbb,d2rhobsigbb,d2rhobsigab,d2rhobsigaa,d2sigaaaa,d2sigaaab,d2sigaabb,d2sigabab,d2sigabbb,d2sigbbbb) 3497integer ixcderv 3498!Input arguments 3499real*8 rhoa,rhob,sigaa,sigbb,sigab,value,d1rhoa,d1rhob,d1sigaa,d1sigbb,d1sigab,d2rhoaa,d2rhobb,d2rhoab,& 3500d2rhoasigaa,d2rhoasigab,d2rhoasigbb,d2rhobsigbb,d2rhobsigab,d2rhobsigaa,d2sigaaaa,d2sigaaab,d2sigaabb,d2sigabab,d2sigabbb,d2sigbbbb 3501!Inputted information 3502real*8 rhoa1(1),rhob1(1),sigmaaa1(1),sigmabb1(1),sigmaab1(1) 3503!Returned information 3504real*8 Xzk(1),Xvrhoa(1),Xvrhob(1),Xvsigmaaa(1),Xvsigmabb(1),Xvsigmaab(1),Xv2rhoa2(1),Xv2rhob2(1),Xv2rhoab(1)& 3505,Xv2rhoasigmaaa(1),Xv2rhoasigmaab(1),Xv2rhoasigmabb(1),Xv2rhobsigmabb(1),Xv2rhobsigmaab(1),Xv2rhobsigmaaa(1)& 3506,Xv2sigmaaa2(1),Xv2sigmaaaab(1),Xv2sigmaaabb(1),Xv2sigmaab2(1),Xv2sigmaabbb(1),Xv2sigmabb2(1) 3507real*8 Czk(1),Cvrhoa(1),Cvrhob(1),Cvsigmaaa(1),Cvsigmabb(1),Cvsigmaab(1),Cv2rhoa2(1),Cv2rhob2(1),Cv2rhoab(1)& 3508,Cv2rhoasigmaaa(1),Cv2rhoasigmaab(1),Cv2rhoasigmabb(1),Cv2rhobsigmabb(1),Cv2rhobsigmaab(1),Cv2rhobsigmaaa(1)& 3509,Cv2sigmaaa2(1),Cv2sigmaaaab(1),Cv2sigmaaabb(1),Cv2sigmaab2(1),Cv2sigmaabbb(1),Cv2sigmabb2(1) 3510real*8 XCzk(1),XCvrhoa(1),XCvrhob(1),XCvsigmaaa(1),XCvsigmabb(1),XCvsigmaab(1),XCv2rhoa2(1),XCv2rhob2(1),XCv2rhoab(1)& 3511,XCv2rhoasigmaaa(1),XCv2rhoasigmaab(1),XCv2rhoasigmabb(1),XCv2rhobsigmabb(1),XCv2rhobsigmaab(1),XCv2rhobsigmaaa(1)& 3512,XCv2sigmaaa2(1),XCv2sigmaaaab(1),XCv2sigmaaabb(1),XCv2sigmaab2(1),XCv2sigmaabbb(1),XCv2sigmabb2(1) 3513 3514rhoa1(1)=rhoa 3515rhob1(1)=rhob 3516sigmaaa1(1)=sigaa 3517sigmabb1(1)=sigbb 3518sigmaab1(1)=sigab 3519!X part 3520if (iDFTxcsel==0.or.iDFTxcsel==80) then 3521 call uks_x_lda(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,Xzk,Xvrhoa,Xvrhob,Xvsigmaaa,Xvsigmabb,Xvsigmaab,Xv2rhoa2,Xv2rhob2,Xv2rhoab,& 3522 Xv2rhoasigmaaa,Xv2rhoasigmaab,Xv2rhoasigmabb,Xv2rhobsigmabb,Xv2rhobsigmaab,Xv2rhobsigmaaa,& 3523 Xv2sigmaaa2,Xv2sigmaaaab,Xv2sigmaaabb,Xv2sigmaab2,Xv2sigmaabbb,Xv2sigmabb2) 3524else if (iDFTxcsel==1.or.iDFTxcsel==81.or.iDFTxcsel==82.or.iDFTxcsel==83) then 3525 call uks_x_b88(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,Xzk,Xvrhoa,Xvrhob,Xvsigmaaa,Xvsigmabb,Xvsigmaab,Xv2rhoa2,Xv2rhob2,Xv2rhoab,& 3526 Xv2rhoasigmaaa,Xv2rhoasigmaab,Xv2rhoasigmabb,Xv2rhobsigmabb,Xv2rhobsigmaab,Xv2rhobsigmaaa,& 3527 Xv2sigmaaa2,Xv2sigmaaaab,Xv2sigmaaabb,Xv2sigmaab2,Xv2sigmaabbb,Xv2sigmabb2) 3528else if (iDFTxcsel==2.or.iDFTxcsel==84) then 3529 call uks_x_pbe(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,Xzk,Xvrhoa,Xvrhob,Xvsigmaaa,Xvsigmabb,Xvsigmaab,Xv2rhoa2,Xv2rhob2,Xv2rhoab,& 3530 Xv2rhoasigmaaa,Xv2rhoasigmaab,Xv2rhoasigmabb,Xv2rhobsigmabb,Xv2rhobsigmaab,Xv2rhobsigmaaa,& 3531 Xv2sigmaaa2,Xv2sigmaaaab,Xv2sigmaaabb,Xv2sigmaab2,Xv2sigmaabbb,Xv2sigmabb2) 3532else if (iDFTxcsel==3.or.iDFTxcsel==85) then 3533 call uks_x_pw91(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,Xzk,Xvrhoa,Xvrhob,Xvsigmaaa,Xvsigmabb,Xvsigmaab,Xv2rhoa2,Xv2rhob2,Xv2rhoab,& 3534 Xv2rhoasigmaaa,Xv2rhoasigmaab,Xv2rhoasigmabb,Xv2rhobsigmabb,Xv2rhobsigmaab,Xv2rhobsigmaaa,& 3535 Xv2sigmaaa2,Xv2sigmaaaab,Xv2sigmaaabb,Xv2sigmaab2,Xv2sigmaabbb,Xv2sigmabb2) 3536end if 3537!C part 3538if (iDFTxcsel==30.or.iDFTxcsel==80) then 3539 call uks_c_vwn5(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,Czk,Cvrhoa,Cvrhob,Cvsigmaaa,Cvsigmabb,Cvsigmaab,Cv2rhoa2,Cv2rhob2,Cv2rhoab,& 3540 Cv2rhoasigmaaa,Cv2rhoasigmaab,Cv2rhoasigmabb,Cv2rhobsigmabb,Cv2rhobsigmaab,Cv2rhobsigmaaa,& 3541 Cv2sigmaaa2,Cv2sigmaaaab,Cv2sigmaaabb,Cv2sigmaab2,Cv2sigmaabbb,Cv2sigmabb2) 3542else if (iDFTxcsel==31.or.iDFTxcsel==81) then 3543 call uks_c_p86(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,Czk,Cvrhoa,Cvrhob,Cvsigmaaa,Cvsigmabb,Cvsigmaab,Cv2rhoa2,Cv2rhob2,Cv2rhoab,& 3544 Cv2rhoasigmaaa,Cv2rhoasigmaab,Cv2rhoasigmabb,Cv2rhobsigmabb,Cv2rhobsigmaab,Cv2rhobsigmaaa,& 3545 Cv2sigmaaa2,Cv2sigmaaaab,Cv2sigmaaabb,Cv2sigmaab2,Cv2sigmaabbb,Cv2sigmabb2) 3546else if (iDFTxcsel==32.or.iDFTxcsel==82) then 3547 call uks_c_lyp(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,Czk,Cvrhoa,Cvrhob,Cvsigmaaa,Cvsigmabb,Cvsigmaab,Cv2rhoa2,Cv2rhob2,Cv2rhoab,& 3548 Cv2rhoasigmaaa,Cv2rhoasigmaab,Cv2rhoasigmabb,Cv2rhobsigmabb,Cv2rhobsigmaab,Cv2rhobsigmaaa,& 3549 Cv2sigmaaa2,Cv2sigmaaaab,Cv2sigmaaabb,Cv2sigmaab2,Cv2sigmaabbb,Cv2sigmabb2) 3550else if (iDFTxcsel==33.or.iDFTxcsel==83.or.iDFTxcsel==85) then 3551 call uks_c_pw91(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,Czk,Cvrhoa,Cvrhob,Cvsigmaaa,Cvsigmabb,Cvsigmaab,Cv2rhoa2,Cv2rhob2,Cv2rhoab,& 3552 Cv2rhoasigmaaa,Cv2rhoasigmaab,Cv2rhoasigmabb,Cv2rhobsigmabb,Cv2rhobsigmaab,Cv2rhobsigmaaa,& 3553 Cv2sigmaaa2,Cv2sigmaaaab,Cv2sigmaaabb,Cv2sigmaab2,Cv2sigmaabbb,Cv2sigmabb2) 3554else if (iDFTxcsel==34.or.iDFTxcsel==84) then 3555 call uks_c_pbe(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,Czk,Cvrhoa,Cvrhob,Cvsigmaaa,Cvsigmabb,Cvsigmaab,Cv2rhoa2,Cv2rhob2,Cv2rhoab,& 3556 Cv2rhoasigmaaa,Cv2rhoasigmaab,Cv2rhoasigmabb,Cv2rhobsigmabb,Cv2rhobsigmaab,Cv2rhobsigmaaa,& 3557 Cv2sigmaaa2,Cv2sigmaaaab,Cv2sigmaaabb,Cv2sigmaab2,Cv2sigmaabbb,Cv2sigmabb2) 3558end if 3559!Whole XC 3560if (iDFTxcsel==70) then 3561 call uks_xc_b97(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,XCzk,XCvrhoa,XCvrhob,XCvsigmaaa,XCvsigmabb,XCvsigmaab,XCv2rhoa2,XCv2rhob2,XCv2rhoab,& 3562 XCv2rhoasigmaaa,XCv2rhoasigmaab,XCv2rhoasigmabb,XCv2rhobsigmabb,XCv2rhobsigmaab,XCv2rhobsigmaaa,& 3563 XCv2sigmaaa2,XCv2sigmaaaab,XCv2sigmaaabb,XCv2sigmaab2,XCv2sigmaabbb,XCv2sigmabb2) 3564else if (iDFTxcsel==71) then 3565 call uks_xc_hcth407(ixcderv,1,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,XCzk,XCvrhoa,XCvrhob,XCvsigmaaa,XCvsigmabb,XCvsigmaab,XCv2rhoa2,XCv2rhob2,XCv2rhoab,& 3566 XCv2rhoasigmaaa,XCv2rhoasigmaab,XCv2rhoasigmabb,XCv2rhobsigmabb,XCv2rhobsigmaab,XCv2rhobsigmaaa,& 3567 XCv2sigmaaa2,XCv2sigmaaaab,XCv2sigmaaabb,XCv2sigmaab2,XCv2sigmaabbb,XCv2sigmabb2) 3568else if (iDFTxcsel>=80.and.iDFTxcsel<99) then 3569 XCzk=Xzk+Czk 3570 XCvrhoa=Xvrhoa+Cvrhoa 3571 XCvrhob=Xvrhob+Cvrhob 3572 XCvsigmaaa=Xvsigmaaa+Cvsigmaaa 3573 XCvsigmabb=Xvsigmabb+Cvsigmabb 3574 XCvsigmaab=Xvsigmaab+Cvsigmaab 3575 XCv2rhoa2=Xv2rhoa2+Cv2rhoa2 3576 XCv2rhob2=Xv2rhob2+Cv2rhob2 3577 XCv2rhoab=Xv2rhoab+Cv2rhoab 3578 XCv2rhoasigmaaa=Xv2rhoasigmaaa+Cv2rhoasigmaaa 3579 XCv2rhoasigmaab=Xv2rhoasigmaab+Cv2rhoasigmaab 3580 XCv2rhoasigmabb=Xv2rhoasigmabb+Cv2rhoasigmabb 3581 XCv2rhobsigmabb=Xv2rhobsigmabb+Cv2rhobsigmabb 3582 XCv2rhobsigmaab=Xv2rhobsigmaab+Cv2rhobsigmaab 3583 XCv2rhobsigmaaa=Xv2rhobsigmaaa+Cv2rhobsigmaaa 3584 XCv2sigmaaa2= Xv2sigmaaa2+Cv2sigmaaa2 3585 XCv2sigmaaaab=Xv2sigmaaaab+Cv2sigmaaaab 3586 XCv2sigmaaabb=Xv2sigmaaabb+Cv2sigmaaabb 3587 XCv2sigmaab2= Xv2sigmaab2+Cv2sigmaab2 3588 XCv2sigmaabbb=Xv2sigmaabbb+Cv2sigmaabbb 3589 XCv2sigmabb2= Xv2sigmabb2+Cv2sigmabb2 3590end if 3591 3592if (iDFTxcsel<30) then 3593 value=Xzk(1) 3594 d1rhoa=Xvrhoa(1) 3595 d1rhob=Xvrhob(1) 3596 d1sigaa=Xvsigmaaa(1) 3597 d1sigbb=Xvsigmabb(1) 3598 d1sigab=Xvsigmaab(1) 3599 d2rhoaa=Xv2rhoa2(1) 3600 d2rhobb=Xv2rhob2(1) 3601 d2rhoab=Xv2rhoab(1) 3602 d2rhoasigaa=Xv2rhoasigmaaa(1) 3603 d2rhoasigab=Xv2rhoasigmaab(1) 3604 d2rhoasigbb=Xv2rhoasigmabb(1) 3605 d2rhobsigbb=Xv2rhobsigmabb(1) 3606 d2rhobsigab=Xv2rhobsigmaab(1) 3607 d2rhobsigaa=Xv2rhobsigmaaa(1) 3608 d2sigaaaa=Xv2sigmaaa2(1) 3609 d2sigaaab=Xv2sigmaaaab(1) 3610 d2sigaabb=Xv2sigmaaabb(1) 3611 d2sigabab=Xv2sigmaab2(1) 3612 d2sigabbb=Xv2sigmaabbb(1) 3613 d2sigbbbb=Xv2sigmabb2(1) 3614else if (iDFTxcsel<70) then 3615 value=Czk(1) 3616 d1rhoa=Cvrhoa(1) 3617 d1rhob=Cvrhob(1) 3618 d1sigaa=Cvsigmaaa(1) 3619 d1sigbb=Cvsigmabb(1) 3620 d1sigab=Cvsigmaab(1) 3621 d2rhoaa=Cv2rhoa2(1) 3622 d2rhobb=Cv2rhob2(1) 3623 d2rhoab=Cv2rhoab(1) 3624 d2rhoasigaa=Cv2rhoasigmaaa(1) 3625 d2rhoasigab=Cv2rhoasigmaab(1) 3626 d2rhoasigbb=Cv2rhoasigmabb(1) 3627 d2rhobsigbb=Cv2rhobsigmabb(1) 3628 d2rhobsigab=Cv2rhobsigmaab(1) 3629 d2rhobsigaa=Cv2rhobsigmaaa(1) 3630 d2sigaaaa=Cv2sigmaaa2(1) 3631 d2sigaaab=Cv2sigmaaaab(1) 3632 d2sigaabb=Cv2sigmaaabb(1) 3633 d2sigabab=Cv2sigmaab2(1) 3634 d2sigabbb=Cv2sigmaabbb(1) 3635 d2sigbbbb=Cv2sigmabb2(1) 3636else if (iDFTxcsel<100) then 3637 value=XCzk(1) 3638 d1rhoa=XCvrhoa(1) 3639 d1rhob=XCvrhob(1) 3640 d1sigaa=XCvsigmaaa(1) 3641 d1sigbb=XCvsigmabb(1) 3642 d1sigab=XCvsigmaab(1) 3643 d2rhoaa=XCv2rhoa2(1) 3644 d2rhobb=XCv2rhob2(1) 3645 d2rhoab=XCv2rhoab(1) 3646 d2rhoasigaa=XCv2rhoasigmaaa(1) 3647 d2rhoasigab=XCv2rhoasigmaab(1) 3648 d2rhoasigbb=XCv2rhoasigmabb(1) 3649 d2rhobsigbb=XCv2rhobsigmabb(1) 3650 d2rhobsigab=XCv2rhobsigmaab(1) 3651 d2rhobsigaa=XCv2rhobsigmaaa(1) 3652 d2sigaaaa=XCv2sigmaaa2(1) 3653 d2sigaaab=XCv2sigmaaaab(1) 3654 d2sigaabb=XCv2sigmaaabb(1) 3655 d2sigabab=XCv2sigmaab2(1) 3656 d2sigabbb=XCv2sigmaabbb(1) 3657 d2sigbbbb=XCv2sigmabb2(1) 3658end if 3659end subroutine 3660 3661 3662 3663!!---- The distance from a point (x,y,z) to the nearest atom in the array 3664real*8 function surfana_di(x,y,z,nlen,atmlist) 3665real*8 x,y,z 3666integer nlen,atmlist(nlen) 3667dist2min=1D100 3668do iatm=1,ncenter 3669 if (any(atmlist==iatm)) then !The atom is in the list 3670 dist2=(a(iatm)%x-x)**2+(a(iatm)%y-y)**2+(a(iatm)%z-z)**2 3671 if (dist2<dist2min) dist2min=dist2 3672 end if 3673end do 3674surfana_di=dsqrt(dist2min) 3675end function 3676 3677!!---- The distance from a point (x,y,z) to the nearest atom not in the array 3678real*8 function surfana_de(x,y,z,nlen,atmlist) 3679real*8 x,y,z 3680integer nlen,atmlist(nlen) 3681dist2min=1D100 3682do iatm=1,ncenter 3683 if (all(atmlist/=iatm)) then !The atom is not in the list 3684 dist2=(a(iatm)%x-x)**2+(a(iatm)%y-y)**2+(a(iatm)%z-z)**2 3685 if (dist2<dist2min) dist2min=dist2 3686 end if 3687end do 3688surfana_de=dsqrt(dist2min) 3689end function 3690 3691!!---- Normalized contact distance, defined in terms of de, di and the vdW radii of the atoms 3692real*8 function surfana_norm(x,y,z,nlen,atmlist) 3693real*8 x,y,z 3694integer nlen,atmlist(nlen) 3695dist2minin=1D100 !The nearest distance to atoms inside 3696dist2minext=1D100 !The nearest distance to atoms outside 3697do iatm=1,ncenter 3698 dist2=(a(iatm)%x-x)**2+(a(iatm)%y-y)**2+(a(iatm)%z-z)**2 3699 if (any(atmlist==iatm)) then !Atoms inside 3700 if (dist2<dist2minin) then 3701 dist2minin=dist2 3702 iminin=iatm 3703 end if 3704 else !Atoms outside 3705 if (dist2<dist2minext) then 3706 dist2minext=dist2 3707 iminext=iatm 3708 end if 3709 end if 3710end do 3711di=dsqrt(dist2minin) 3712de=dsqrt(dist2minext) 3713rvdwin=vdwr(a(iminin)%index) 3714rvdwext=vdwr(a(iminext)%index) 3715surfana_norm=(di-rvdwin)/rvdwin+(de-rvdwext)/rvdwext 3716end function 3717 3718 3719 3720 3721!!-------- PAEM, potential acting on one electron in a molecule, defined by Zhongzhi Yang in JCC,35,965(2014) 3722!If itype=1, evaluate the XC potential based on pair density; if =2, it will be equivalent to DFT XC potential 3723real*8 function PAEM(x,y,z,itype) 3724integer itype 3725real*8 x,y,z,wfnval(nmo),GTFint(nprims,nprims) 3726!Evaluate electron contribution to ESP 3727call genGTFattmat(x,y,z,GTFint) 3728rhopot=0 3729do imo=1,nmo 3730 do iprim=1,nprims 3731 do jprim=1,nprims 3732 rhopot=rhopot+MOocc(imo)*CO(imo,iprim)*CO(imo,jprim)*GTFint(iprim,jprim) 3733 end do 3734 end do 3735end do 3736!Evaluate XC potential 3737if (itype==1) then !Based on Muller approximation form 3738 call orbderv(1,1,nmo,x,y,z,wfnval) 3739 rho=sum(MOocc(1:nmo)*wfnval(1:nmo)**2) 3740 xcpot=0 3741 if (wfntype==0.or.wfntype==3) then !Close-shell 3742 do imo=1,nmo 3743 if (MOocc(imo)==0D0) cycle 3744 do jmo=1,nmo 3745 if (MOocc(jmo)==0D0) cycle 3746 tmpval=dsqrt(MOocc(imo)*MOocc(jmo))*wfnval(imo)*wfnval(jmo) 3747 do iprim=1,nprims 3748 do jprim=1,nprims 3749 xcpot=xcpot-tmpval*CO(imo,iprim)*CO(jmo,jprim)*GTFint(iprim,jprim) 3750 end do 3751 end do 3752 end do 3753 end do 3754 else if (wfntype==1.or.wfntype==4) then !Unrestricted open-shell 3755 do ialphaend=nmo,1,-1 !Find the ending index of alpha MO 3756 if (MOtype(ialphaend)==1) exit 3757 end do 3758 do imo=1,ialphaend !Alpha part 3759 do jmo=1,ialphaend 3760 tmpval=dsqrt(MOocc(imo)*MOocc(jmo))*wfnval(imo)*wfnval(jmo) 3761 do iprim=1,nprims 3762 do jprim=1,nprims 3763 xcpot=xcpot-tmpval*CO(imo,iprim)*CO(jmo,jprim)*GTFint(iprim,jprim) 3764 end do 3765 end do 3766 end do 3767 end do 3768 do imo=ialphaend+1,nmo !Beta part 3769 do jmo=ialphaend+1,nmo 3770 tmpval=dsqrt(MOocc(imo)*MOocc(jmo))*wfnval(imo)*wfnval(jmo) 3771 do iprim=1,nprims 3772 do jprim=1,nprims 3773 xcpot=xcpot-tmpval*CO(imo,iprim)*CO(jmo,jprim)*GTFint(iprim,jprim) 3774 end do 3775 end do 3776 end do 3777 end do 3778 else if (wfntype==2) then !Restricted open-shell 3779 do imo=1,nmo !Alpha part 3780 if (MOocc(imo)==0) cycle 3781 do jmo=1,nmo 3782 if (MOocc(jmo)==0) cycle 3783 tmpval=wfnval(imo)*wfnval(jmo) !Every occupied ROHF MOs contributes one alpha electron 3784 do iprim=1,nprims 3785 do jprim=1,nprims 3786 xcpot=xcpot-tmpval*CO(imo,iprim)*CO(jmo,jprim)*GTFint(iprim,jprim) 3787 end do 3788 end do 3789 end do 3790 end do 3791 do imo=1,nmo !Beta part 3792 if (MOocc(imo)/=2D0) cycle 3793 do jmo=1,nmo 3794 if (MOocc(jmo)/=2D0) cycle 3795 tmpval=wfnval(imo)*wfnval(jmo) !Every doubly occupied ROHF MOs contributes one beta electron 3796 do iprim=1,nprims 3797 do jprim=1,nprims 3798 xcpot=xcpot-tmpval*CO(imo,iprim)*CO(jmo,jprim)*GTFint(iprim,jprim) 3799 end do 3800 end do 3801 end do 3802 end do 3803 end if 3804 xcpot=xcpot/rho 3805 3806else if (itype==2) then !Directly using DFT XC potential 3807 xcpot=DFTxcpot(x,y,z) 3808end if 3809 3810PAEM=-nucesp(x,y,z)+rhopot+xcpot 3811end function 3812 3813 3814 3815!!----- Angle between the eigenvectors of rho and the plane defined by option 4 of main function 1000 3816! The plane is represented by global variables pleA,pleB,pleC,pleD 3817! ivec=1/2/3 means the eigenvector corresponding to the first/second/third highest eigenvalue is calculated 3818real*8 function Ang_rhoeigvec_ple(x,y,z,ivec) 3819use util 3820integer ivec 3821real*8 x,y,z,eigvecmat(3,3),eigval(3),eigvaltmp(3),elehess(3,3),elegrad(3) 3822call calchessmat_dens(2,x,y,z,elerho,elegrad,elehess) 3823call diagmat(elehess,eigvecmat,eigval,100,1D-10) 3824eigvaltmp=eigval 3825call sort(eigvaltmp) !From small to large 3826call invarr(eigvaltmp,1,3) !1/2/3=large to small 3827do i=1,3 3828 if (eigval(i)==eigvaltmp(ivec)) exit 3829end do 3830Ang_rhoeigvec_ple=vecang(eigvecmat(1,i),eigvecmat(2,i),eigvecmat(3,i),pleA,pleB,pleC) 3831end function 3832 3833 3834!!------ Local electron correlation function (DOI: 10.1021/acs.jctc.7b00293) 3835!itype=1: Local total electron correlation function 3836!itype=2: Local dynamic electron correlation function 3837!itype=3: Local nondynamic electron correlation function 3838real*8 function localcorr(x,y,z,itype) 3839integer itype 3840real*8 x,y,z,wfnval(nmo),occ(nmo) 3841call orbderv(1,1,nmo,x,y,z,wfnval) 3842localcorr=0D0 3843if (wfntype==3) then 3844 occ=MOocc/2 3845 where(occ>1) occ=1 !Remove unphysical larger than unity occupation number 3846 where(occ<0) occ=0 !Remove unphysical negative occupation number 3847 if (itype==1) then 3848 do i=1,nmo 3849 localcorr=localcorr+ dsqrt(occ(i)*(1-occ(i)))*wfnval(i)**2 3850 end do 3851 localcorr=localcorr/4 3852 else if (itype==2) then 3853 do i=1,nmo 3854 localcorr=localcorr+ ( dsqrt(occ(i)*(1-occ(i))) - 2*occ(i)*(1-occ(i)) ) *wfnval(i)**2 3855 end do 3856 localcorr=localcorr/4 3857 else if (itype==3) then 3858 do i=1,nmo 3859 localcorr=localcorr+ occ(i)*(1-occ(i))*wfnval(i)**2 3860 end do 3861 localcorr=localcorr/2 3862 end if 3863 localcorr=localcorr*2 3864else if (wfntype==4) then 3865 occ=MOocc 3866 where(occ>1) occ=1 3867 where(occ<0) occ=0 3868 if (itype==1) then 3869 do i=1,nmo 3870 localcorr=localcorr+ dsqrt(occ(i)*(1-occ(i)))*wfnval(i)**2 3871 end do 3872 localcorr=localcorr/4 3873 else if (itype==2) then 3874 do i=1,nmo 3875 localcorr=localcorr+ ( dsqrt(occ(i)*(1-occ(i))) - 2*occ(i)*(1-occ(i)) ) *wfnval(i)**2 3876 end do 3877 localcorr=localcorr/4 3878 else if (itype==3) then 3879 do i=1,nmo 3880 localcorr=localcorr+ occ(i)*(1-occ(i))*wfnval(i)**2 3881 end do 3882 localcorr=localcorr/2 3883 end if 3884end if 3885end function 3886 3887 3888!For visually examine functions used in DFRT2.0 project 3889real*8 function funcvalLSB(x,y,z,itype) 3890integer itype 3891integer,parameter :: nfunc=7 3892real*8 x,y,z,valarr(nfunc),rho,gradrho(3) 3893iexpcutoffold=expcutoff 3894expcutoff=1 3895call valaryyLSB(x,y,z,valarr,rho,gradrho) 3896funcvalLSB=valarr(itype) 3897expcutoff=iexpcutoffold 3898end function 3899 3900 3901 3902 3903 3904!!=================================================================================!! 3905!! Below codes are contributed by Arshad Mehmood for calculating EDR(r;d) and D(r) !! 3906!!=================================================================================!! 3907 3908!!----- For three-point numerical fit to evaluate D(r) 3909subroutine three_point_interpolation(n,x,y,xmax,ymax) 3910integer, intent(in) :: n 3911real*8, intent(in) :: x(max_edr_exponents),y(max_edr_exponents) 3912real*8, intent(out) :: xmax,ymax 3913integer i , imax 3914real*8 x1,x2,x3, y1,y2,y3, a,b 3915100 format ('XXX ',3F9.5) 3916ymax = -1.0d0 3917imax = -1 3918do i=1,n 3919 if(y(i) .gt. ymax) then 3920 ymax = y(i) 3921 imax = i 3922 endif 3923 end do 3924if(imax<1 .or. imax>n) then 3925 write(*,*) "Error: Bad imax" 3926 call EXIT() 3927end if 3928if(imax .eq. 1 .or. imax.eq.n) then 3929 xmax = x(imax)**(-0.5d0) 3930 return 3931endif 3932x1 = x(imax-1)**(-0.5d0) 3933x2 = x(imax )**(-0.5d0) 3934x3 = x(imax+1)**(-0.5d0) 3935y1 = y(imax-1) 3936y2 = y(imax ) 3937y3 = y(imax+1) 3938a = ( (y3-y2)/(x3-x2) -(y2-y1)/(x2-x1) )/(x3-x1) 3939b = ( (y3-y2)/(x3-x2)*(x2-x1) + (y2-y1)/(x2-x1)*(x3-x2) )& 3940 /(x3-x1) 3941xmax = x2 - b/(2d0*a) 3942ymax = y2 - b**(2d0)/(4d0*a) 3943end subroutine 3944 3945!!----- Function to calculate EDR(r,d) 3946!J. Chem. Phys. 141, 144104(2014), J. Chem. Theory Comput. 12, 79(2016), Angew. Chem. Int. Ed. 56, 6878(2017) 3947real*8 function edr(x,y,z) 3948real*8 :: ed(max_edr_exponents),edrval(max_edr_exponents) 3949nedr=1 3950ed(1)=dedr**(-2.0d0) 3951call EDRcal(1,x,y,z,nedr,ed,edrval,edrdmaxval) 3952edr=edrval(1) 3953end function 3954 3955!!----- Function to calculate D(r) 3956!J. Chem. Theory Comput. 12, 3185(2016), Phys. Chem. Chem. Phys. 17, 18305(2015) 3957real*8 function edrdmax(x,y,z) 3958real*8 :: ed(max_edr_exponents),edrdmaxval,edrval(max_edr_exponents) 3959real*8 :: edrexponent 3960integer iedr 3961edrexponent = edrastart 3962do iedr=1,nedr 3963 ed(iedr)=edrexponent 3964 edrexponent=edrexponent/edrainc 3965end do 3966call EDRcal(2,x,y,z,nedr,ed,edrval,edrdmaxval) 3967edrdmax=edrdmaxval 3968end function 3969 3970!!----- Working routine used to evaluate EDR(r;d) and D(r) 3971subroutine EDRcal(runtype,x,y,z,nedr,ed,edrval,edrdmaxval) 3972real*8, intent(in) :: x,y,z,ed(max_edr_exponents) 3973integer, intent(in) :: nedr 3974real*8, intent(out):: edrdmaxval,edrval(max_edr_exponents) 3975real*8 :: rho,dmaxdummy 3976real*8 :: psi(nmo),AMUVal(max_edr_exponents),Bint(nmo,max_edr_exponents) 3977real*8 :: xamu(3,max_edr_exponents),amu0(max_edr_exponents) 3978integer :: j,ixyz,i,iedr,runtype 3979edrval = 0D0 3980if (runtype==2) then 3981 edrdmaxval = 0D0 3982end if 3983 3984rho=fdens(x,y,z) 3985 3986if(rho.gt.1D-10) then 3987 psi = 0d0 3988 Bint=0d0 3989 do j=1,nprims 3990 ix=type2ix(b(j)%functype) 3991 iy=type2iy(b(j)%functype) 3992 iz=type2iz(b(j)%functype) 3993 ep=b(j)%exp 3994 sftx=x-a(b(j)%center)%x 3995 sfty=y-a(b(j)%center)%y 3996 sftz=z-a(b(j)%center)%z 3997 sftx2=sftx*sftx 3998 sfty2=sfty*sfty 3999 sftz2=sftz*sftz 4000 rr=sftx2+sfty2+sftz2 4001 expterm=0.0 4002 amu0 = 0d0 4003 if (expcutoff>0.or.-ep*rr>expcutoff) then 4004 expterm=exp(-ep*rr) 4005 do iedr=1,nedr 4006 amu0(iedr)=(2d0*ed(iedr)/pi)**(3d0/4d0) *(pi/(ep+ed(iedr)))**(3d0/2d0)& 4007 * exp(-ep*ed(iedr)/(ep+ed(iedr))*rr) 4008 end do 4009 end if 4010 4011 if (expterm==0D0) cycle 4012 4013 GTFval=sftx**ix *sfty**iy *sftz**iz *expterm 4014 do iedr=1,nedr 4015 do ixyz=1,3 4016 ival=ix 4017 sftval=sftx 4018 if(ixyz.eq.2) then 4019 ival=iy 4020 sftval=sfty 4021 else if(ixyz.eq.3) then 4022 ival=iz 4023 sftval=sftz 4024 end if 4025 If(ival.eq.0) then 4026 xamu(ixyz,iedr)=1d0 4027 else if(ival.eq.1) then 4028 xamu(ixyz,iedr)=sftval*ed(iedr)/(ed(iedr)+ep) 4029 else If(ival.eq.2) then 4030 xamu(ixyz,iedr)=(sftval*ed(iedr)/(ed(iedr)+ep))**2d0 + 1d0/(2d0*(ed(iedr)+ep)) 4031 else If(ival.eq.3) then 4032 xamu(ixyz,iedr)=(sftval*ed(iedr)/(ed(iedr)+ep))**3d0 + sftval*3d0*ed(iedr)/(2d0*(ed(iedr)+ep)**2d0) 4033 else If(ival.eq.4) then 4034 xamu(ixyz,iedr)=sftval**4d0* (ed(iedr)/(ed(iedr)+ep))**4d0 + sftval**2d0 *3d0*ed(iedr)**2d0/(ed(iedr)+ep)**3d0 & 4035 + 3d0/(4d0*(ed(iedr)+ep)**2d0) 4036 else If(ival.eq.5) then 4037 xamu(ixyz,iedr)=sftval**5d0* ed(iedr)**5d0/(ed(iedr)+ep)**5d0 + sftval**3d0* 5d0*ed(iedr)**3d0/(ed(iedr)+ep)**4d0 & 4038 + sftval *15d0*ed(iedr)/(4d0*(ed(iedr)+ep)**3d0) 4039 else 4040 write(*,*) "Angular momentum out of range" 4041 Call EXIT() 4042 end if 4043 end do 4044 end do 4045 do iedr=1,nedr 4046 AMUVal(iedr)=amu0(iedr)*xamu(1,iedr)*xamu(2,iedr)*xamu(3,iedr) 4047 end do 4048 4049 do i=1,nmo 4050 if (nint(MOocc(i)).GE.1D0) then 4051 psi(i)=psi(i)+co(i,j)*GTFval 4052 end if 4053 end do 4054 4055 do iedr=1,nedr 4056 do i=1,nmo 4057 if (nint(MOocc(i)).GE.1D0) then 4058 Bint(i,iedr)=Bint(i,iedr)+co(i,j)*AMUVal(iedr) 4059 end if 4060 end do 4061 end do 4062 end do 4063 4064 edrval = 0d0 4065 do i=1,nmo 4066 do iedr=1,nedr 4067 edrval(iedr)=edrval(iedr)+psi(i)*Bint(i,iedr) 4068 end do 4069 end do 4070 if (runtype==2) then 4071 call three_point_interpolation(nedr,ed,edrval,edmax,dmaxdummy) 4072 edrdmaxval=edmax 4073 else if (runtype==1) then 4074 do iedr=1,nedr 4075 edrval(iedr)=edrval(iedr)*rho**(-0.5D0) 4076 end do 4077 else 4078 write(*,*) "EDRcal runtype out of range" 4079 call EXIT() 4080 end if 4081end if 4082end subroutine 4083 4084 4085 4086end module 4087