1subroutine popana 2use defvar 3use util 4implicit real*8 (a-h,o-z) 5character c2000tmp*2000 6 7if (ifragcontri==1) then 8 write(*,*) "Population analysis function couldn't be used combining with self-defined fragment" 9else 10 do while(.true.) 11 imodwfnold=imodwfn 12 write(*,*) " ============== Population analysis ==============" 13 write(*,*) "-2 Calculate interaction energy between fragments based on atomic charges" 14 if (.not.allocated(frag1)) then 15 write(*,*) "-1 Define fragment" 16 else 17 write(*,"(a,i5)") "-1 Redefine fragment, current atoms:",size(frag1) 18 end if 19 write(*,*) "0 Return" 20 write(*,*) "1 Hirshfeld population" 21 write(*,*) "2 Voronoi deformation density (VDD) population" 22 !Not available because integration error of below two methods by means of Becke integration are too large 23 ! write(*,*) "3 Integrate electron density in voronoi cell" 24 ! write(*,*) "4 Adjusted method 3 by Rousseau et al." 25 if (allocated(cobasa)) then 26 write(*,*) "5 Mulliken atom & basis function population analysis" 27 write(*,*) "6 Lowdin population analysis" 28 write(*,*) "7 Modified Mulliken population defined by Ros & Schuit (SCPA)" 29 write(*,*) "8 Modified Mulliken population defined by Stout & Politzer" 30 write(*,*) "9 Modified Mulliken population defined by Bickelhaupt" 31 end if 32 write(*,*) "10 Becke atomic charge with atomic dipole moment correction" 33 write(*,*) "11 Atomic dipole corrected Hirshfeld population (ADCH)" 34 write(*,*) "12 CHELPG ESP fitting charge" 35 write(*,*) "13 Merz-Kollmann (MK) ESP fitting charge" 36 write(*,*) "14 AIM charge" 37 write(*,*) "15 Hirshfeld-I charge" 38 write(*,*) "16 CM5 charge" 39 write(*,*) "17 Electronegativity Equalization Method (EEM) charge" 40 read(*,*) ipopsel 41 42 if (ipopsel==0) then 43 if (allocated(frag1)) deallocate(frag1) 44 return 45 else if (ipopsel==-2) then 46 if (ifiletype/=4) then 47 write(*,*) "Error: You must use .chg file for this function!" 48 write(*,*) "Press ENTER to continue" 49 read(*,*) 50 cycle 51 end if 52 call coulint_atmchg 53 else if (ipopsel==-1) then 54 if (allocated(frag1)) then 55 write(*,*) "Atoms in current fragment:" 56 write(*,"(13i6)") frag1 57 write(*,"(a)") " Input 0 to keep unchanged, or redefine fragment, e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute the fragment" 58 else 59 write(*,"(a)") " Input atomic indices to define the fragment. e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will constitute the fragment" 60 end if 61 read(*,"(a)") c2000tmp 62 if (c2000tmp(1:1)/='0') then 63 if (allocated(frag1)) deallocate(frag1) 64 call str2arr(c2000tmp,nfragatm) 65 allocate(frag1(nfragatm)) 66 call str2arr(c2000tmp,nfragatm,frag1) 67 if (any(frag1>ncenter)) then 68 write(*,*) "Error: Some atomic indices exceeded valid range! Please define again" 69 write(*,*) 70 deallocate(frag1) 71 end if 72 end if 73 else if (ipopsel==1) then 74 write(*,*) "Citation: Theor. Chim. Acta. (Berl), 44, 129-138 (1977)" 75 call spacecharge(1) 76 else if (ipopsel==2) then 77 write(*,*) "Citation: Organomet., 15, 2923-2931" 78 write(*,*) "Citation: J. Comput. Chem., 25, 189-210 (2004)" 79 call spacecharge(2) 80 else if (ipopsel==3) then 81 call spacecharge(3) 82 else if (ipopsel==4) then 83 write(*,*) "Citation: J. Mol. Struct.(Theochem), 538, 235-238 (2001)" 84 call spacecharge(4) 85 else if (ipopsel==5) then 86 write(*,*) "0 Return" 87 write(*,*) "1 Output Mulliken population/charges and decompose them to MO's contribution" 88 write(*,*) "2 Output gross atomic population matrix and decompose it" 89 write(*,*) "3 Output gross basis function population matrix and decompose it" 90 read(*,*) ipopsel2 91 if (ipopsel2==0) cycle 92 call MPA(ipopsel2) 93 else if (ipopsel==6) then 94 call symmortho 95 call MPA(1) 96 call dealloall 97 call readinfile(firstfilename,1) !Current wavefunction has been altered, recover the initial state 98 else if (ipopsel==7) then 99 call MMPA(1) 100 else if (ipopsel==8) then 101 call MMPA(2) 102 else if (ipopsel==9) then 103 call Bickelhaupt 104 else if (ipopsel==10) then 105 call spacecharge(5) 106 else if (ipopsel==11) then 107 write(*,*) "Citation of ADCH: Tian Lu, Feiwu Chen, J. Theor. Comput. Chem., 11, 163 (2011)" 108 call spacecharge(6) 109 else if (ipopsel==12) then 110 call fitESP(1) 111 else if (ipopsel==13) then 112 call fitESP(2) 113 else if (ipopsel==14) then 114 write(*,"(a)") " NOTE: AIM charges cannot be calculated in present module but can be calculated in basin analysis module, & 115 please check the example given in Section 4.17.1 of the manual on how to do this" 116 else if (ipopsel==15) then 117 if (iautointgrid==1) then 118 radpot=30 119 sphpot=170 120 if (any(a%index>18)) radpot=40 121 if (any(a%index>36)) radpot=50 122 if (any(a%index>54)) radpot=60 123 end if 124 call Hirshfeld_I_wrapper(1) 125 else if (ipopsel==16) then 126 call spacecharge(7) 127 else if (ipopsel==17) then 128 call EEM 129 end if 130 if (imodwfnold==1.and.(ipopsel==1.or.ipopsel==2.or.ipopsel==6.or.ipopsel==11)) then !1,2,6,11 are the methods need to reload the initial wavefunction 131 write(*,"(a)") " Note: The wavefunction file has been reloaded, your previous modifications on occupation number will be ignored" 132 end if 133 end do 134end if 135end subroutine 136 137 138!!---------- Calculate Coulomb interaction between two fragment based on atomic charges in .chg file 139subroutine coulint_atmchg 140use defvar 141use util 142character c2000tmp*2000 143do while(.true.) 144 write(*,*) "Input atom list for fragment 1, e.g. 1,4,6-9" 145 write(*,*) "Input 0 can exit" 146 read(*,"(a)") c2000tmp 147 if (c2000tmp(1:1)=='0') exit 148 call str2arr(c2000tmp,nfrag1) 149 if (allocated(frag1)) deallocate(frag1) 150 allocate(frag1(nfrag1)) 151 call str2arr(c2000tmp,nfrag1,frag1) 152 write(*,*) "Input atom list for fragment 2, e.g. 1,4,6-9" 153 read(*,"(a)") c2000tmp 154 call str2arr(c2000tmp,nfrag2) 155 if (allocated(frag2)) deallocate(frag2) 156 allocate(frag2(nfrag2)) 157 call str2arr(c2000tmp,nfrag2,frag2) 158 eleint=0 159 do iatmidx=1,nfrag1 160 iatm=frag1(iatmidx) 161 do jatmidx=1,nfrag2 162 jatm=frag2(jatmidx) 163 eleint=eleint+a(iatm)%charge*a(jatm)%charge/distmat(iatm,jatm) 164 end do 165 end do 166 write(*,"(' Electrostatic interaction energy:',f12.2,' kJ/mol',f12.2,' kcal/mol',/)") eleint*au2kJ,eleint*au2kcal 167end do 168end subroutine 169 170 171 172!!---------- Modified Mulliken population analysis defined by Bickelhaupt 173subroutine Bickelhaupt 174use defvar 175use util 176implicit real*8 (a-h,o-z) 177character selectyn,chgfilename*200 178real*8,target :: atmeletot(ncenter),atmelea(ncenter) 179real*8,pointer :: tmpmat(:,:),tmpele(:) 180do itime=1,2 181 if (itime==1) then 182 tmpele=>atmeletot 183 tmpmat=>Ptot 184 else if (itime==2) then 185 tmpele=>atmelea 186 tmpmat=>Palpha 187 end if 188 tmpele=0D0 189 do i=1,ncenter 190 do j=basstart(i),basend(i) 191 cross=0D0 192 do k=1,nbasis !This method equalvalent to use diagonal element of density matrix to partition nondiagonal element of P*S matrix 193 if (k/=j) then 194 if (tmpmat(j,j)+tmpmat(k,k)/=0D0) then 195 cross=cross+tmpmat(j,j)/(tmpmat(j,j)+tmpmat(k,k))*tmpmat(j,k)*Sbas(j,k) 196 else 197 cross=cross+tmpmat(j,k)*Sbas(j,k)/2D0 !Use equivalent partition, when denominator is zero 198 end if 199 end if 200 end do 201 tmpele(i)=tmpele(i)+tmpmat(j,j)+2*cross !Plus electrons localized in basis and partitioned cross term 202 end do 203 end do 204 if (wfntype==0.or.wfntype==3) exit 205end do 206 207if (wfntype==0.or.wfntype==3) then 208 do iatm=1,ncenter 209 write(*,"(' Atom',i6,'(',a2,')',' Population:',f10.5,' Atomic charge:',f10.5)") iatm,a(iatm)%name,atmeletot(iatm),a(iatm)%charge-atmeletot(iatm) 210 end do 211 write(*,"(' Total net charge:',f10.5)") sum(a(:)%charge)-sum(atmeletot(:)) 212else if (wfntype==1.or.wfntype==2.or.wfntype==4) then 213 write(*,*) " Atom Alpha_pop. Beta_pop. Spin_pop. Atomic charge" 214 totbetapop=0D0 215 do iatm=1,ncenter 216 betapop=atmeletot(iatm)-atmelea(iatm) 217 totbetapop=totbetapop+betapop 218 write(*,"(i6,'(',a2,')',3f13.5,f16.5)") iatm,a(iatm)%name,atmelea(iatm),betapop,atmelea(iatm)-betapop,a(iatm)%charge-atmeletot(iatm) 219 end do 220 write(*,"(' Total net charge:',f10.5,' Total spin electrons:',f10.5)") sum(a(:)%charge)-sum(atmeletot(:)),sum(atmelea(:))-totbetapop 221end if 222 223!Show fragment charge 224if (allocated(frag1)) write(*,"(/,' Fragment charge:',f12.6)") sum(a(frag1)%charge-atmeletot(frag1)) 225 226call path2filename(firstfilename,chgfilename) 227write(*,*) 228write(*,"(a)") " If output atom with charges to "//trim(chgfilename)//".chg in current folder? (y/n)" 229read(*,*) selectyn 230if (selectyn=="y".or.selectyn=="Y") then 231 open(10,file=trim(chgfilename)//".chg",status="replace") 232 do i=1,ncenter 233 write(10,"(a4,4f12.6)") a(i)%name,a(i)%x*b2a,a(i)%y*b2a,a(i)%z*b2a,a(i)%charge-atmeletot(i) 234 end do 235 close(10) 236 write(*,"(a)") " Result have been saved to "//trim(chgfilename)//".chg in current folder" 237 write(*,"(a)") " Columns ranging from 1 to 5 are name,X,Y,Z,charge, respectively. The ith row corresponds to the ith atom. The unit is Angstrom" 238end if 239end subroutine 240 241 242 243!!---------- Modified Mulliken population analysis 244! isel=1 :Defined by Ros & Schuit (SCPA)" 245! isel=2 :Defined by Stout & Politzer 246subroutine MMPA(isel) 247use defvar 248use util 249implicit real*8 (a-h,o-z) 250integer isel 251character selectyn,chgfilename*200 252real*8,target :: atmelea(ncenter),atmeleb(ncenter) 253real*8,pointer :: tmpmat(:,:),tmpele(:) 254atmelea=0D0 255atmeleb=0D0 256do itime=1,2 257 do imo=1,nbasis 258 if (itime==1) then 259 irealmo=imo 260 tmpele=>atmelea 261 tmpmat=>cobasa 262 else if (itime==2) then 263 irealmo=imo+nbasis 264 tmpele=>atmeleb 265 tmpmat=>cobasb 266 end if 267 if (MOocc(irealmo)==0D0) cycle 268 if (isel==1) allbassqr=sum(tmpmat(:,imo)**2) 269 do i=1,ncenter 270 atmbassqr=sum(tmpmat(basstart(i):basend(i),imo)**2) 271 if (isel==1) then 272 tmpele(i)=tmpele(i)+MOocc(irealmo)*atmbassqr/allbassqr 273 else if (isel==2) then 274 cross=0D0 275 do ii=basstart(i),basend(i) 276 do jj=1,nbasis 277 denomin=tmpmat(ii,imo)**2+tmpmat(jj,imo)**2 278 if (jj/=ii.and.denomin>=1D-120) cross=cross+tmpmat(ii,imo)**2/denomin*tmpmat(ii,imo)*tmpmat(jj,imo)*Sbas(ii,jj) 279 end do 280 end do 281 tmpele(i)=tmpele(i)+MOocc(irealmo)*(atmbassqr+2*cross) 282 end if 283 end do 284 end do 285 if (wfntype==0.or.wfntype==3) exit 286end do 287if (wfntype==0.or.wfntype==3) then 288 do iatm=1,ncenter 289 write(*,"(' Atom',i6,'(',a2,')',' Population:',f10.5,' Atomic charge:',f10.5)") iatm,a(iatm)%name,atmelea(iatm),a(iatm)%charge-atmelea(iatm) 290 end do 291 write(*,"(' Total net charge:',f10.5)") sum(a(:)%charge)-sum(atmelea(:)) 292 if (allocated(frag1)) write(*,"(/,' Fragment charge:',f12.6)") sum(a(frag1)%charge-atmelea(frag1)) 293else if (wfntype==1.or.wfntype==2.or.wfntype==4) then 294 write(*,*) " Atom Alpha_pop. Beta_pop. Spin_pop. Atomic charge" 295 do iatm=1,ncenter 296 write(*,"(i6,'(',a2,')',3f13.5,f16.5)") iatm,a(iatm)%name,atmelea(iatm),atmeleb(iatm),atmelea(iatm)-atmeleb(iatm),a(iatm)%charge-atmelea(iatm)-atmeleb(iatm) 297 end do 298 write(*,"(' Total net charge:',f10.5,' Total spin electrons:',f10.5)") sum(a(:)%charge)-sum(atmelea(:))-sum(atmeleb(:)),sum(atmelea(:))-sum(atmeleb(:)) 299 if (allocated(frag1)) write(*,"(/,' Fragment charge:',f12.6)") sum(a(frag1)%charge-atmelea(frag1)-atmeleb(frag1)) 300end if 301 302call path2filename(firstfilename,chgfilename) 303write(*,*) 304write(*,"(a)") " If output atom coordinates with charges to "//trim(chgfilename)//".chg file in current folder? (y/n)" 305read(*,*) selectyn 306if (selectyn=="y".or.selectyn=="Y") then 307 open(10,file=trim(chgfilename)//".chg",status="replace") 308 do i=1,ncenter 309 if (wfntype==0.or.wfntype==3) write(10,"(a4,4f12.6)") a(i)%name,a(i)%x*b2a,a(i)%y*b2a,a(i)%z*b2a,a(i)%charge-atmelea(i) 310 if (wfntype==1.or.wfntype==2.or.wfntype==4) write(10,"(a4,4f12.6)") a(i)%name,a(i)%x*b2a,a(i)%y*b2a,a(i)%z*b2a,a(i)%charge-atmelea(i)--atmeleb(i) 311 end do 312 close(10) 313 write(*,"(a)") "Result have been saved to "//trim(chgfilename)//".chg in current folder" 314 write(*,"(a)") "Columns ranging from 1 to 5 are name,X,Y,Z,charge respectively, unit is Angstrom" 315end if 316end subroutine 317 318 319 320!!--------- Mulliken/Lowdin population analysis & decompose to MO contribution 321! isel=1 Output Mulliken/Lowdin charge and decompose it to MO contribution" 322! isel=2 Output gross atomic population matrix and decompose it" 323! isel=3 Output gross basis function population matrix and decompose it" 324!Note: If doing Lowdin population, density matrix and overlap matrix should be transformed first before invoking this routine 325subroutine MPA(isel) 326use defvar 327use util 328implicit real*8 (a-h,o-z) 329real*8 MOcenmat(nbasis,ncenter),groatmmat(ncenter+1,ncenter),atmele(ncenter),charge(ncenter) 330real*8,pointer :: ptmat(:,:) 331real*8,allocatable :: tmpmat(:,:),basmata(:,:),angorbpop(:,:),angorbpopa(:,:),angorbpopb(:,:) 332character selectyn,corbnum*6,cOcc*12,chgfilename*200 333integer isel 334 335if (isel==1.or.isel==2) then 336 allocate(tmpmat(nbasis,nbasis),basmata(nbasis,nbasis)) !basmata stores basis gross population of alpha part 337 do itime=1,3 !total elec, alpha elec, beta elec 338 if (itime==1) tmpmat=Ptot*Sbas 339 if (itime==2) then 340 tmpmat=Palpha*Sbas 341 basmata=tmpmat !Backup 342 end if 343 if (itime==3) tmpmat=Pbeta*Sbas 344 !Calculate gross atomic population matrix 345 do i=1,ncenter 346 do j=1,ncenter 347 accum=0D0 348 do ii=basstart(i),basend(i) 349 do jj=basstart(j),basend(j) 350 accum=accum+tmpmat(ii,jj) 351 end do 352 end do 353 groatmmat(i,j)=accum 354 end do 355 end do 356 do i=1,ncenter !Stored atom populations to the last row 357 groatmmat(ncenter+1,i)=sum(groatmmat(1:ncenter,i)) 358 end do 359 totelec=0D0 360 361 if (isel==2) then !Output gross atomic population matrix 362 if (itime==1) call showmatgau(groatmmat,"Total gross atomic population matrix",0,"f14.8") 363 if (itime==2) call showmatgau(groatmmat,"Alpha gross atomic population matrix",0,"f14.8") 364 if (itime==3) call showmatgau(groatmmat,"Beta gross atomic population matrix",0,"f14.8") 365 else if (isel==1) then !Contract gross atomic population matrix and output population in each basis function/shell/atom 366 if (wfntype==0.or.wfntype==3) then !Notice that only perform once (itime=1) 367 allocate(angorbpop(ncenter,0:5)) !Record the population number in each angular moment orbitals, up to H 368 angorbpop=0D0 369 write(*,*) "Population of basis functions:" 370 write(*,"(' Basis Type Atom Shell Population')") 371 do ibas=1,nbasis 372 write(*,"(i6,3x,a,i5,a,i5,f13.5)") ibas,GTFtype2name(bastype(ibas)),bascen(ibas),'('//a(bascen(ibas))%name//')',basshell(ibas),sum(tmpmat(ibas,:)) 373 end do 374 write(*,*) 375 write(*,*) "Population of shells of basis functions:" 376 do ish=1,nshell 377 shellpop=0D0 378 do ibas=1,nbasis 379 if (basshell(ibas)==ish) then 380 iatm=bascen(ibas) !Which atom this shell attribute to 381 shellpop=shellpop+sum(tmpmat(ibas,:)) 382 end if 383 end do 384 write(*,"(' Shell',i6,' Type: ',a,' in atom',i5,'(',a,') :',f9.5)") ish,shtype2name(shtype(ish)),iatm,a(iatm)%name,shellpop 385 iangtmp=abs(shtype(ish)) 386 angorbpop(iatm,iangtmp)=angorbpop(iatm,iangtmp)+shellpop 387 end do 388 write(*,*) 389 write(*,*) "Population of each type of angular moment orbitals:" 390 do iatm=1,ncenter 391 write(*,"(' Atom',i6,'(',a2,')',' s:',f7.4,' p:',f7.4,' d:',f7.4,' f:',f7.4,' g:',f7.4,' h:',f7.4)") iatm,a(iatm)%name,angorbpop(iatm,:) 392 end do 393 write(*,"(' Sum s:',f9.4,' p:',f9.4,' d:',f9.4,' f:',f9.4,' g:',f9.4,' h:',f9.4)") & 394 sum(angorbpop(:,0)),sum(angorbpop(:,1)),sum(angorbpop(:,2)),sum(angorbpop(:,3)),sum(angorbpop(:,4)),sum(angorbpop(:,5)) 395 write(*,*) 396 write(*,*) "Population of atoms:" 397 do iatm=1,ncenter 398 charge(iatm)=a(iatm)%charge-groatmmat(ncenter+1,iatm) 399 write(*,"(' Atom',i6,'(',a2,')',' Population:',f11.5,' Net charge:',f11.5)") iatm,a(iatm)%name,groatmmat(ncenter+1,iatm),charge(iatm) 400 end do 401 write(*,"(' Total net charge:',f10.5)") sum(a(:)%charge)-sum(groatmmat(ncenter+1,:)) 402 else if ((wfntype==1.or.wfntype==2.or.wfntype==4).and.itime==3) then !For unrestrict wfn, at last "itime" cycle print result 403 allocate(angorbpopa(ncenter,0:5),angorbpopb(ncenter,0:5)) 404 angorbpopa=0D0 405 angorbpopb=0D0 406 write(*,*) "Population of basis functions:" 407 write(*,"(' Basis Type Atom Shell Alpha_pop. Beta_pop. Total_pop. Spin_pop.')") 408 do ibas=1,nbasis !Note: Currently, tmpmat stores basis gross population of beta part, basmata stores alpha part 409 baspopa=sum(basmata(ibas,:)) 410 baspopb=sum(tmpmat(ibas,:)) 411 write(*,"(i6,3x,a,i5,a,i5,1x,4f12.5)") ibas,GTFtype2name(bastype(ibas)),bascen(ibas),& 412 '('//a(bascen(ibas))%name//')',basshell(ibas),baspopa,baspopb,baspopa+baspopb,baspopa-baspopb 413 end do 414 write(*,*) 415 write(*,*) "Population of shells:" 416 write(*,*) "Shell Type Atom Alpha_pop. Beta_pop. Total_pop. Spin_pop." 417 do ish=1,nshell 418 shellpopa=0D0 419 shellpopb=0D0 420 do ibas=1,nbasis 421 if (basshell(ibas)==ish) then 422 iatm=bascen(ibas) !Which atom this shell attribute to 423 shellpopa=shellpopa+sum(basmata(ibas,:)) 424 shellpopb=shellpopb+sum(tmpmat(ibas,:)) 425 end if 426 end do 427 write(*,"(i5,5x,a,i7,'(',a,')' ,4f12.5)") ish,shtype2name(shtype(ish)),iatm,a(iatm)%name,shellpopa,shellpopb,shellpopa+shellpopb,shellpopa-shellpopb 428 iangtmp=abs(shtype(ish)) 429 angorbpopa(iatm,iangtmp)=angorbpopa(iatm,iangtmp)+shellpopa 430 angorbpopb(iatm,iangtmp)=angorbpopb(iatm,iangtmp)+shellpopb 431 end do 432 write(*,*) 433 write(*,*) "Population of each type of angular moment orbitals:" 434 write(*,*) " Atom Type Alpha_pop. Beta_pop. Total_pop. Spin_pop." 435 do iatm=1,ncenter 436 if (angorbpopa(iatm,0)/=0D0.or.angorbpopb(iatm,0)/=0D0) write(*,"(i6,'(',a2,') s',4f13.5)") & 437 iatm,a(iatm)%name,angorbpopa(iatm,0),angorbpopb(iatm,0),angorbpopa(iatm,0)+angorbpopb(iatm,0),angorbpopa(iatm,0)-angorbpopb(iatm,0) 438 if (angorbpopa(iatm,1)/=0D0.or.angorbpopb(iatm,1)/=0D0) write(*,"(' p',4f13.5)") & 439 angorbpopa(iatm,1),angorbpopb(iatm,1),angorbpopa(iatm,1)+angorbpopb(iatm,1),angorbpopa(iatm,1)-angorbpopb(iatm,1) 440 if (angorbpopa(iatm,2)/=0D0.or.angorbpopb(iatm,2)/=0D0) write(*,"(' d',4f13.5)") & 441 angorbpopa(iatm,2),angorbpopb(iatm,2),angorbpopa(iatm,2)+angorbpopb(iatm,2),angorbpopa(iatm,2)-angorbpopb(iatm,2) 442 if (angorbpopa(iatm,3)/=0D0.or.angorbpopb(iatm,3)/=0D0) write(*,"(' f',4f13.5)") & 443 angorbpopa(iatm,3),angorbpopb(iatm,3),angorbpopa(iatm,3)+angorbpopb(iatm,3),angorbpopa(iatm,3)-angorbpopb(iatm,3) 444 if (angorbpopa(iatm,4)/=0D0.or.angorbpopb(iatm,4)/=0D0) write(*,"(' g',4f13.5)") & 445 angorbpopa(iatm,4),angorbpopb(iatm,4),angorbpopa(iatm,4)+angorbpopb(iatm,4),angorbpopa(iatm,4)-angorbpopb(iatm,4) 446 if (angorbpopa(iatm,5)/=0D0.or.angorbpopb(iatm,5)/=0D0) write(*,"(' h',4f13.5)") & 447 angorbpopa(iatm,5),angorbpopb(iatm,5),angorbpopa(iatm,5)+angorbpopb(iatm,5),angorbpopa(iatm,5)-angorbpopb(iatm,5) 448 end do 449 write(*,*) 450 if (sum(angorbpopa(:,0))/=0D0.or.sum(angorbpopb(:,0))/=0D0) write(*,"(' Total s',4f13.5)") & 451 sum(angorbpopa(:,0)),sum(angorbpopb(:,0)),sum(angorbpopa(:,0))+sum(angorbpopb(:,0)),sum(angorbpopa(:,0))-sum(angorbpopb(:,0)) 452 if (sum(angorbpopa(:,1))/=0D0.or.sum(angorbpopb(:,1))/=0D0) write(*,"(' p',4f13.5)") & 453 sum(angorbpopa(:,1)),sum(angorbpopb(:,1)),sum(angorbpopa(:,1))+sum(angorbpopb(:,1)),sum(angorbpopa(:,1))-sum(angorbpopb(:,1)) 454 if (sum(angorbpopa(:,2))/=0D0.or.sum(angorbpopb(:,2))/=0D0) write(*,"(' d',4f13.5)") & 455 sum(angorbpopa(:,2)),sum(angorbpopb(:,2)),sum(angorbpopa(:,2))+sum(angorbpopb(:,2)),sum(angorbpopa(:,2))-sum(angorbpopb(:,2)) 456 if (sum(angorbpopa(:,3))/=0D0.or.sum(angorbpopb(:,3))/=0D0) write(*,"(' f',4f13.5)") & 457 sum(angorbpopa(:,3)),sum(angorbpopb(:,3)),sum(angorbpopa(:,3))+sum(angorbpopb(:,3)),sum(angorbpopa(:,3))-sum(angorbpopb(:,3)) 458 if (sum(angorbpopa(:,4))/=0D0.or.sum(angorbpopb(:,4))/=0D0) write(*,"(' g',4f13.5)") & 459 sum(angorbpopa(:,4)),sum(angorbpopb(:,4)),sum(angorbpopa(:,4))+sum(angorbpopb(:,4)),sum(angorbpopa(:,4))-sum(angorbpopb(:,4)) 460 if (sum(angorbpopa(:,5))/=0D0.or.sum(angorbpopb(:,5))/=0D0) write(*,"(' h',4f13.5)") & 461 sum(angorbpopa(:,5)),sum(angorbpopb(:,5)),sum(angorbpopa(:,5))+sum(angorbpopb(:,5)),sum(angorbpopa(:,5))-sum(angorbpopb(:,5)) 462 write(*,*) 463 write(*,*) "Population of atoms:" 464 write(*,*) " Atom Alpha_pop. Beta_pop. Spin_pop. Atomic charge" 465 totspinpop=0D0 466 do iatm=1,ncenter 467 alphaele=atmele(iatm) 468 betaele=groatmmat(ncenter+1,iatm) 469 charge(iatm)=a(iatm)%charge-(alphaele+betaele) 470 write(*,"(i6,'(',a2,')',3f13.5,f16.5)") iatm,a(iatm)%name,alphaele,betaele,alphaele-betaele,charge(iatm) 471 totspinpop=totspinpop+alphaele-betaele 472 end do 473 write(*,"(' Total net charge:',f10.5,' Total spin electrons:',f10.5)") sum(charge),totspinpop 474 end if 475 if (itime==2) atmele(:)=groatmmat(ncenter+1,:) !Store alpha occupation of each atom to a temporary array 476 end if 477 if (wfntype==0.or.wfntype==3) exit !RHF or ROHF, don't continue to process alpha & beta respectively 478 end do 479 480 !Show fragment charge 481 if (allocated(frag1)) write(*,"(/,' Fragment charge:',f12.6)") sum(charge(frag1)) 482 483 !Asking user if decompose result 484 write(*,*) 485 if (isel==1) then 486 write(*,*) "Decompose atomic charges to MO contribution? (y/n)" 487 else if (isel==2) then 488 write(*,*) "The last row is the sum of corresponding column elements (Atomic population)" 489 write(*,*) 490 write(*,"(a)") "Decompose to MO contribution and write to groatmdcp.txt in current folder? (y/n)" 491 end if 492 read(*,*) selectyn 493 if (selectyn=='y'.or.selectyn=='Y') then 494 if (isel==1) then 495 write(*,*) "The (i,j) elements means the contribution to the jth atoms from the ith MO" 496 else if (isel==2) then 497 open(10,file="groatmdcp.txt",status="replace") 498 write(10,*) "The last row is the sum of corresponding column elements" 499 write(10,*) 500 end if 501 do itime=1,2 502 MOcenmat=0D0 503 if (itime==1) ptmat=>cobasa 504 if (itime==2) ptmat=>cobasb 505 do imo=1,nbasis 506 if (itime==1) irealmo=imo 507 if (itime==2) irealmo=imo+nbasis 508 write(corbnum,"(i6)") imo 509 write(cOcc,"(f12.8)") MOocc(irealmo) 510 if (MOocc(irealmo)==0D0) cycle 511 512 do i=1,ncenter 513 do j=1,ncenter 514 accum=0D0 515 do ii=basstart(i),basend(i) 516 do jj=basstart(j),basend(j) 517 accum=accum+MOocc(irealmo)*ptmat(ii,imo)*ptmat(jj,imo)*Sbas(ii,jj) 518 end do 519 end do 520 groatmmat(i,j)=accum 521 end do 522 end do 523 do i=1,ncenter 524 groatmmat(ncenter+1,i)=sum(groatmmat(1:ncenter,i)) 525 end do 526 if (isel==1) then 527 MOcenmat(imo,:)=groatmmat(ncenter+1,:) 528 else if (isel==2) then 529 if (wfntype==0.or.wfntype==2.or.wfntype==3) then 530 call showmatgau(groatmmat,"Orbital"//corbnum//" Occ:"//cOcc,0,"f14.8",10) 531 else if (itime==1.and.(wfntype==1.or.wfntype==4)) then 532 call showmatgau(groatmmat,"Alpha Orbital"//corbnum//" Occ:"//cOcc,0,"f14.8",10) 533 else if (itime==2.and.(wfntype==1.or.wfntype==4)) then 534 call showmatgau(groatmmat,"Beta Orbital"//corbnum//" Occ:"//cOcc,0,"f14.8",10) 535 end if 536 write(10,*) 537 end if 538 end do 539 if (isel==1) then 540 if (wfntype==0.or.wfntype==2) call showmatgau(MOcenmat,"",0,"f14.8",6,nint(naelec)) 541 if (itime==1.and.wfntype==1) call showmatgau(MOcenmat,"Alpha part",0,"f14.8",6,nint(naelec)) 542 if (itime==2.and.wfntype==1) call showmatgau(MOcenmat,"Beta part",0,"f14.8",6,nint(nbelec)) 543 if (wfntype==3) call showmatgau(MOcenmat,"",0,"f14.8") 544 if (itime==1.and.wfntype==4) call showmatgau(MOcenmat,"Alpha part",0,"f14.8") 545 if (itime==2.and.wfntype==4) call showmatgau(MOcenmat,"Beta part",0,"f14.8") 546 end if 547 if (wfntype==0.or.wfntype==2.or.wfntype==3) exit !ROHF needn't to separate to alpha and beta 548 end do 549 if (isel==2) close(10) 550 if (isel==2) write(*,*) "Done!" 551 end if 552 553 !If calculating atomic charges, asking user if output it 554 if (isel==1) then 555 call path2filename(firstfilename,chgfilename) 556 write(*,*) 557 write(*,"(a)") " If output atom with charges to "//trim(chgfilename)//".chg in current folder? (y/n)" 558 read(*,*) selectyn 559 if (selectyn=="y".or.selectyn=="Y") then 560 open(10,file=trim(chgfilename)//".chg",status="replace") 561 do i=1,ncenter 562 if (wfntype==0.or.wfntype==3) write(10,"(a4,4f12.6)") a(i)%name,a(i)%x*b2a,a(i)%y*b2a,a(i)%z*b2a,charge(i) 563 if (wfntype==1.or.wfntype==2.or.wfntype==4) write(10,"(a4,4f12.6)") a(i)%name,a(i)%x*b2a,a(i)%y*b2a,a(i)%z*b2a,charge(i) 564 end do 565 close(10) 566 write(*,"(a)") " Result have been saved to "//trim(chgfilename)//".chg in current folder" 567 write(*,"(a)") " Columns ranging from 1 to 5 are name,X,Y,Z,charge respectively, unit is Angstrom" 568 end if 569 end if 570 571else if (isel==3) then 572 call showmatgau(Ptot*Sbas,"Total gross basis function population matrix",0,"f14.8") 573 if (wfntype==1.or.wfntype==2.or.wfntype==4) then 574 call showmatgau(Palpha*Sbas,"Alpha gross basis function population matrix",0,"f14.8") 575 call showmatgau(Pbeta*Sbas,"Beta gross basis function population matrix",0,"f14.8") 576 end if 577 write(*,*) "The (i,j) element means ��[imo] Occ(imo)*C(i,imo)*C(j,imo)*S(i,j)" 578 write(*,*) 579 write(*,"(a)") "Decompose to MO contribution and write to grobasdcp.txt in current folder? (y/n)" 580 read(*,*) selectyn 581 if (selectyn=='y'.or.selectyn=='Y') then 582 open(10,file="grobasdcp.txt",status="replace") 583 write(10,*) "Notes:" 584 write(10,*) "The (i,j) element means C(i,imo)*C(j,imo)*S(i,j)" 585 write(10,*) "The last row is the sum of corresponding column elements" 586 write(10,*) 587 allocate(tmpmat(nbasis+1,nbasis)) 588 do itime=1,2 589 if (itime==1) ptmat=>cobasa 590 if (itime==2) ptmat=>cobasb 591 do imo=1,nbasis 592 if (itime==1) irealmo=imo 593 if (itime==2) irealmo=imo+nbasis 594 write(corbnum,"(i6)") imo 595 write(cOcc,"(f12.8)") MOocc(irealmo) 596 if (MOocc(irealmo)==0D0) cycle 597 598 do ii=1,nbasis 599 do jj=1,nbasis 600 tmpmat(ii,jj)=MOocc(irealmo)*ptmat(ii,imo)*ptmat(jj,imo)*Sbas(ii,jj) 601 end do 602 end do 603 do j=1,nbasis 604 tmpmat(nbasis+1,j)=sum(tmpmat(1:nbasis,j)) 605 end do 606 if (wfntype==0.or.wfntype==2.or.wfntype==3) then 607 call showmatgau(tmpmat,"Orbital"//corbnum//" Occ:"//cOcc,0,"f14.8",10) 608 else if (itime==1.and.(wfntype==1.or.wfntype==4)) then 609 call showmatgau(tmpmat,"Alpha Orbital"//corbnum//" Occ:"//cOcc,0,"f14.8",10) 610 else if (itime==2.and.(wfntype==1.or.wfntype==4)) then 611 call showmatgau(tmpmat,"Beta Orbital"//corbnum//" Occ:"//cOcc,0,"f14.8",10) 612 end if 613 write(10,*) 614 end do 615 if (wfntype==0.or.wfntype==2.or.wfntype==3) exit 616 end do 617 close(10) 618 write(*,*) "Done!" 619 end if 620end if 621end subroutine 622 623 624 625 626!!-------------- Calculate charge based on space partition method 627!1=Hirshfeld, 2=VDD, 3=Integrate electron density in voronoi cell 628!4=Adjusted method 3 by Rousseau et al., 5= Becke with/without ADC, 6= ADCH, 7= CM5 629subroutine spacecharge(chgtype) 630use defvar 631use function 632use util 633implicit real*8(a-h,o-z) 634integer chgtype 635real*8 molrho(radpot*sphpot),promol(radpot*sphpot),tmpdens(radpot*sphpot),beckeweigrid(radpot*sphpot),selfdens(radpot*sphpot) 636type(content) gridatm(radpot*sphpot),gridatmorg(radpot*sphpot) 637real*8 atmdipx(ncenter),atmdipy(ncenter),atmdipz(ncenter),charge(ncenter) 638real*8 :: covr_becke(0:nelesupp) !covalent radii used for Becke population 639character selectyn,chgfilename*200 640character radfilename*800 641integer :: nbeckeiter=3 642 643if (chgtype==5) then !Select atomic radii for Becke population 644 covr_becke=covr_TianLu 645 iraddefine=2 646 do while(.true.) 647 write(*,*) "-1 Return" 648 write(*,*) "0 Start calculation of Becke charge!" 649 if (iraddefine==0) write(*,*) "1 Select the definition of atomic radii, current: Custom" 650 if (iraddefine==1) write(*,*) "1 Select the definition of atomic radii, current: CSD" 651 if (iraddefine==2) write(*,*) "1 Select the definition of atomic radii, current: Modified CSD" 652 if (iraddefine==3) write(*,*) "1 Select the definition of atomic radii, current: Pyykko" 653 if (iraddefine==4) write(*,*) "1 Select the definition of atomic radii, current: Suresh" 654 if (iraddefine==5) write(*,*) "1 Select the definition of atomic radii, current: Hugo" 655 write(*,"(a,i2)") " 2 Set the number of iterations for defining Becke atomic space, current:",nbeckeiter 656 write(*,*) "10 Read radii from external file" 657 write(*,*) "11 Modify current radii by manual input" 658 write(*,*) "12 Print current radii list" 659 read(*,*) isel 660 if (isel==-1) then 661 return 662 else if (isel==0) then 663 exit 664 else if (isel==1) then 665 write(*,*) "1 Use CSD radii (Dalton Trans., 2008, 2832-2838)" 666 write(*,*) "2 Use the modified version of CSD radii defined by Tian Lu (Recommended)" 667 write(*,*) "3 Use Pyykko radii (Chem. Eur.-J., 15, 186-197)" 668 write(*,*) "4 Use Suresh radii (J. Phys. Chem. A, 105, 5940-5944)" 669 write(*,*) "5 Use Hugo radii (Chem. Phys. Lett., 480, 127-131)" 670 read(*,*) iselrad 671 if (iselrad==1) then 672 covr_becke=covr 673 iraddefine=1 674 else if (iselrad==2) then 675 covr_becke=covr_TianLu 676 iraddefine=2 677 else if (iselrad==3) then 678 covr_becke=covr_pyy 679 iraddefine=3 680 else if (iselrad==4) then 681 covr_becke=covr_Suresh 682 iraddefine=4 683 else if (iselrad==5) then 684 covr_becke=radii_hugo 685 iraddefine=5 686 end if 687 else if (isel==2) then 688 write(*,*) "Input a number, e.g. 3" 689 read(*,*) nbeckeiter 690 else if (isel==10) then 691 iraddefine=0 692 write(*,"(a)") " About the file format:" 693 write(*,"(a)") " The first line should be the number of elements you want to modify, followed by element indices and radii (in Angstrom), for example:" 694 write(*,"(a)") " 4" 695 write(*,"(a)") " 1 0.35" 696 write(*,"(a)") " 4 1.2" 697 write(*,"(a)") " 5 1.12" 698 write(*,"(a)") " 14 1.63" 699 write(*,*) 700 write(*,*) "Input filename" 701 read(*,"(a)") radfilename 702 inquire(file=radfilename,exist=alive) 703 if (alive.eqv..true.) then 704 open(10,file=radfilename,status="old") 705 read(10,*) nmodrad 706 do irad=1,nmodrad 707 read(10,*) indtmp,radtmp 708 covr_becke(indtmp)=radtmp/b2a 709 end do 710 close(10) 711 write(*,*) "Done!" 712 else 713 write(*,*) "Error: File cannot be found" 714 end if 715 else if (isel==11) then 716 iraddefine=0 717 write(*,*) "Input element index and radius (in Angstrom), e.g. 5,0.84" 718 read(*,*) indtmp,radtmp 719 covr_becke(indtmp)=radtmp/b2a 720 write(*,*) "Done!" 721 else if (isel==12) then 722 do irad=0,nelesupp 723 write(*,"(' Element:',i5,'(',a,') Radius:',f8.3,' Angstrom')") irad,ind2name(irad),covr_becke(irad)*b2a 724 end do 725 write(*,*) 726 end if 727 end do 728end if 729 730!Generate quadrature point and weighs by combination of Gauss-Chebyshev and Lebedev grids 731call gen1cintgrid(gridatmorg,iradcut) 732 733!***** 1=Hirshfeld, 2=VDD, 6=ADCH, 7=CM5 734if (chgtype==1.or.chgtype==2.or.chgtype==6.or.chgtype==7) then 735 write(*,*) "This task requests atomic densities, please select how to obtain them" 736 write(*,*) "1 Use build-in sphericalized atomic densities in free-states (more convenient)" 737 write(*,"(a)") " 2 Provide wavefunction file of involved elements by yourself or invoke Gaussian to automatically calculate them" 738 read(*,*) iatmdensmode 739 if (iatmdensmode==2) call setpromol !In this routine reload first molecule at the end 740 write(*,"(' Radial grids:',i5,' Angular grids:',i5,' Total:',i10)") radpot,sphpot,radpot*sphpot 741 write(*,*) "Calculating, please wait..." 742 write(*,*) 743 call walltime(nwalltime1) 744 do iatm=1,ncenter !Cycle each atom to calculate their charges and dipole 745 call delvirorb(0) !For faster calculation, remove virtual MOs in whole system, will not affect result 746 atmx=a(iatm)%x 747 atmy=a(iatm)%y 748 atmz=a(iatm)%z 749 gridatm%value=gridatmorg%value !Weight in this grid point 750 gridatm%x=gridatmorg%x+atmx !Move quadrature point to actual position in molecule 751 gridatm%y=gridatmorg%y+atmy 752 gridatm%z=gridatmorg%z+atmz 753 !Calculate molecular density first 754nthreads=getNThreads() 755!$OMP parallel do shared(molrho) private(i) num_threads(nthreads) 756 do i=1,radpot*sphpot 757 molrho(i)=fdens(gridatm(i)%x,gridatm(i)%y,gridatm(i)%z) 758 end do 759!$OMP end parallel do 760 !Calc free atomic density to obtain promolecule density 761 promol=0D0 762 if (iatmdensmode==1) then 763 do jatm=1,ncenter 764nthreads=getNThreads() 765!$OMP parallel do shared(tmpdens) private(ipt) num_threads(nthreads) 766 do ipt=1,radpot*sphpot 767 tmpdens(ipt)=calcatmdens(jatm,gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z,0) 768 end do 769!$OMP end parallel do 770 promol=promol+tmpdens 771 if (jatm==iatm) selfdens=tmpdens 772 end do 773 else if (iatmdensmode==2) then 774 do jatm=1,ncenter 775 call dealloall 776 call readwfn(custommapname(jatm),1) 777nthreads=getNThreads() 778!$OMP parallel do shared(tmpdens) private(ipt) num_threads(nthreads) 779 do ipt=1,radpot*sphpot 780 tmpdens(ipt)=fdens(gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z) 781 end do 782!$OMP end parallel do 783 promol=promol+tmpdens 784 if (jatm==iatm) selfdens=tmpdens 785 end do 786 call dealloall 787 call readinfile(firstfilename,1) !Retrieve to first loaded file(whole molecule) to calc real rho again 788 end if 789 !Now we have needed data in hand, calculate atomic charges and atomic dipole moments 790 tmpcharge=0D0 791 dipx=0D0 792 dipy=0D0 793 dipz=0D0 794 if (chgtype==1.or.chgtype==6.or.chgtype==7) then !Hirshfeld, ADCH charge, CM5 charge 795 do i=1,radpot*sphpot 796 if (promol(i)/=0D0) then 797 tmpv=selfdens(i)/promol(i)*molrho(i)*gridatm(i)%value 798 tmpcharge=tmpcharge-tmpv 799 dipx=dipx-(gridatm(i)%x-atmx)*tmpv 800 dipy=dipy-(gridatm(i)%y-atmy)*tmpv 801 dipz=dipz-(gridatm(i)%z-atmz)*tmpv 802 end if 803 end do 804 if (nEDFelec==0) charge(iatm)=a(iatm)%charge+tmpcharge 805 if (nEDFelec>0) charge(iatm)=a(iatm)%index+tmpcharge !EDF is provided 806 else if (chgtype==2) then !VDD charge 807 do i=1,radpot*sphpot !Cycle each grid point of iatm, if the distance between the grid point and other atom is shorter than iatm, weight=0 808 vddwei=1D0 809 discen2=(gridatm(i)%x-atmx)**2+(gridatm(i)%y-atmy)**2+(gridatm(i)%z-atmz)**2 !Distance between this grid and current center atom 810 do jatm=1,ncenter_org !Note: Current wfn is atomic wfn, so use _org suffix 811 if (jatm==iatm) cycle 812 disother2=(gridatm(i)%x-a_org(jatm)%x)**2+(gridatm(i)%y-a_org(jatm)%y)**2+(gridatm(i)%z-a_org(jatm)%z)**2 813 if (disother2<discen2) then 814 vddwei=0D0 !Using this weight is equivalent to using Voronoi cell 815 exit 816 end if 817 end do 818 tmpv=vddwei*(molrho(i)-promol(i))*gridatm(i)%value 819 tmpcharge=tmpcharge-tmpv 820 dipx=dipx-(gridatm(i)%x-atmx)*tmpv 821 dipy=dipy-(gridatm(i)%y-atmy)*tmpv 822 dipz=dipz-(gridatm(i)%z-atmz)*tmpv 823 end do 824 charge(iatm)=tmpcharge 825 end if 826 atmdipx(iatm)=dipx 827 atmdipy(iatm)=dipy 828 atmdipz(iatm)=dipz 829 if (chgtype==1.or.chgtype==6.or.chgtype==7) write(*,"(' Hirshfeld charge of atom ',i5,'(',a2,')',' is',f12.6)") iatm,a_org(iatm)%name,charge(iatm) 830 if (chgtype==2) write(*,"(' VDD charge of atom ',i5,'(',a2,')',' is',f12.6)") iatm,a_org(iatm)%name,charge(iatm) 831 end do 832 833!***** 3=Integrate electron density in Voronoi cell, 4=Adjusted method 3 by Rousseau et al 834else if (chgtype==3.or.chgtype==4) then 835 write(*,"(' Radial grids:',i5,' Angular grids:',i5,' Total:',i10)") radpot,sphpot,radpot*sphpot 836 write(*,*) "Calculating, please wait..." 837 write(*,*) 838 call walltime(nwalltime1) 839 if (chgtype==4) then !vdW radius From J.Mol.Stru.(Theo.) 538,235-238 is not identical to original definition 840 vdwr(1)=0.68D0/b2a 841 !B,C,N,O,F 842 vdwr(5)=1.46D0/b2a 843 vdwr(6)=1.46D0/b2a 844 vdwr(7)=1.39D0/b2a 845 vdwr(8)=1.35D0/b2a 846 vdwr(9)=1.29D0/b2a 847 !P S Cl 848 vdwr(15)=1.78D0/b2a 849 vdwr(16)=1.74D0/b2a 850 vdwr(17)=1.69D0/b2a 851 end if 852 do iatm=1,ncenter 853 tmpcharge=0D0 854 dipx=0D0 855 dipy=0D0 856 dipz=0D0 857 atmx=a(iatm)%x 858 atmy=a(iatm)%y 859 atmz=a(iatm)%z 860 gridatm%value=gridatmorg%value 861 gridatm%x=gridatmorg%x+atmx !Move quadrature point with center of current atom 862 gridatm%y=gridatmorg%y+atmy 863 gridatm%z=gridatmorg%z+atmz 864 do i=1,radpot*sphpot 865 vorwei=1.0D0 866 discen2=(gridatm(i)%x-atmx)**2+(gridatm(i)%y-atmy)**2+(gridatm(i)%z-atmz)**2 !Distance between this grid and current center atom 867 do jatm=1,ncenter !Determine the boundary of cell 868 if (jatm==iatm) cycle 869 disother2=(gridatm(i)%x-a(jatm)%x)**2+(gridatm(i)%y-a(jatm)%y)**2+(gridatm(i)%z-a(jatm)%z)**2 870 if (chgtype==3) then 871 if (disother2<discen2) then 872 vorwei=0.0D0 !Use this weights equivalent to use voronoi cell 873 exit 874 end if 875 else if (chgtype==4) then !Adjusted voronoi 876 vdwra=vdwr(a(iatm)%index) 877 vdwrb=vdwr(a(jatm)%index) 878 RAB=distmat(iatm,jatm) 879 rhoval=(RAB**2+discen2-disother2)/2.0D0/RAB 880 rhoa=vdwra/(vdwra+vdwrb)*RAB 881 if (rhoval>rhoa) then 882 vorwei=0.0D0 883 exit 884 end if 885 end if 886 end do 887 if (vorwei/=0.0D0) then 888 tmpv=vorwei*fdens(gridatm(i)%x,gridatm(i)%y,gridatm(i)%z)*gridatm(i)%value 889 tmpcharge=tmpcharge-tmpv 890 dipx=dipx-(gridatm(i)%x-atmx)*tmpv 891 dipy=dipy-(gridatm(i)%y-atmy)*tmpv 892 dipz=dipz-(gridatm(i)%z-atmz)*tmpv 893 end if 894 end do 895 charge(iatm)=tmpcharge+a(iatm)%charge 896 atmdipx(iatm)=dipx 897 atmdipy(iatm)=dipy 898 atmdipz(iatm)=dipz 899 write(*,"(' The charge of atom ',i5,'(',a2,')',' is',f12.6)") iatm,a(iatm)%name,charge(iatm) 900 end do 901 902!***** Becke population 903else if (chgtype==5) then 904 write(*,"(' Radial grids:',i5,' Angular grids:',i5,' Total:',i10)") radpot,sphpot,radpot*sphpot 905 write(*,*) "Calculating, please wait..." 906 write(*,*) 907 call walltime(nwalltime1) 908 do iatm=1,ncenter !Cycle each atom to calculate their charges and dipole 909 gridatm%value=gridatmorg%value !Weight in this grid point 910 gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule 911 gridatm%y=gridatmorg%y+a(iatm)%y 912 gridatm%z=gridatmorg%z+a(iatm)%z 913nthreads=getNThreads() 914!$OMP parallel do shared(tmpdens) private(i) num_threads(nthreads) 915 do i=1,radpot*sphpot !Calc molecular density first 916 tmpdens(i)=fdens(gridatm(i)%x,gridatm(i)%y,gridatm(i)%z) 917 end do 918!$OMP end parallel do 919 call gen1cbeckewei(iatm,iradcut,gridatm,beckeweigrid) 920 tmpcharge=0D0 921 dipx=0D0 922 dipy=0D0 923 dipz=0D0 924 do i=1+iradcut*sphpot,radpot*sphpot 925 tmpv=tmpdens(i)*beckeweigrid(i)*gridatm(i)%value 926 tmpcharge=tmpcharge-tmpv 927 dipx=dipx-(gridatm(i)%x-a(iatm)%x)*tmpv 928 dipy=dipy-(gridatm(i)%y-a(iatm)%y)*tmpv 929 dipz=dipz-(gridatm(i)%z-a(iatm)%z)*tmpv 930 end do 931 if (nEDFelec==0) charge(iatm)=a(iatm)%charge+tmpcharge 932 if (nEDFelec>0) charge(iatm)=a(iatm)%index+tmpcharge !EDF is provided 933 atmdipx(iatm)=dipx 934 atmdipy(iatm)=dipy 935 atmdipz(iatm)=dipz 936 write(*,"(' Becke charge of atom ',i5,'(',a2,')',' is',f12.6)") iatm,a(iatm)%name,charge(iatm) 937 end do 938end if 939 940write(*,"(' Summing up all charges:',f15.8)") sum(charge) 941write(*,*) 942write(*,*) "Atomic dipole moments (a.u.):" 943do iatm=1,ncenter 944 write(*,"(' Atom ',i5,'(',a2,')',' in X/Y/Z:',3f11.6,' Norm:',f11.6)") iatm,a(iatm)%name,atmdipx(iatm),atmdipy(iatm),atmdipz(iatm),dsqrt(atmdipx(iatm)**2+atmdipy(iatm)**2+atmdipz(iatm)**2) 945end do 946xmoldip=0.0D0 947ymoldip=0.0D0 948zmoldip=0.0D0 949do i=1,ncenter 950 xmoldip=xmoldip+a(i)%x*charge(i) 951 ymoldip=ymoldip+a(i)%y*charge(i) 952 zmoldip=zmoldip+a(i)%z*charge(i) 953end do 954totdip=dsqrt(xmoldip**2+ymoldip**2+zmoldip**2) 955write(*,"(' Total dipole moment from atomic charges:',f12.6,' a.u.')") totdip 956write(*,"(' X/Y/Z of dipole from atomic charge:',3f12.6,' a.u.')") xmoldip,ymoldip,zmoldip 957totatmdip=dsqrt(sum(atmdipx)**2+sum(atmdipy)**2+sum(atmdipz)**2) 958write(*,"(' Total atomic dipole moment:',f12.6,' a.u.')") totatmdip 959write(*,"(' X/Y/Z of total atomic dipole:',3f12.6,' a.u.')") sum(atmdipx),sum(atmdipy),sum(atmdipz) 960corrdipx=xmoldip+sum(atmdipx) !Corresponding to actual molecular dipole moment derived from molecular density 961corrdipy=ymoldip+sum(atmdipy) 962corrdipz=zmoldip+sum(atmdipz) 963realdip=dsqrt(corrdipx**2+corrdipy**2+corrdipz**2) 964 965if (chgtype==5) call doADC(atmdipx,atmdipy,atmdipz,charge,realdip,5) 966if (chgtype==6) call doADC(atmdipx,atmdipy,atmdipz,charge,realdip,6) 967if (chgtype==7) call doCM5(charge) 968 969!Show fragment charge 970if (allocated(frag1)) write(*,"(/,' Fragment charge:',f12.6)") sum(charge(frag1)) 971 972write(*,*) 973call walltime(nwalltime2) 974write(*,"(' Calculation took up',i8,' seconds wall clock time')") nwalltime2-nwalltime1 975 976call path2filename(firstfilename,chgfilename) 977write(*,*) 978write(*,"(a)") " If output atoms with charges to "//trim(chgfilename)//".chg in current folder? (y/n)" 979read(*,*) selectyn 980if (selectyn=="y".or.selectyn=="Y") then 981 open(10,file=trim(chgfilename)//".chg",status="replace") 982 do i=1,ncenter 983 write(10,"(a4,4f12.6)") a(i)%name,a(i)%x*b2a,a(i)%y*b2a,a(i)%z*b2a,charge(i) 984 end do 985 close(10) 986 write(*,"(a)") " Result have been saved to "//trim(chgfilename)//".chg in current folder" 987 write(*,"(a)") " Columns ranging from 1 to 5 are name,X,Y,Z,charge respectively, unit is Angstrom" 988end if 989end subroutine 990 991 992!!------ Calculate atomic dipole moment corrected charge based on existing atomic charge (charge) and atomic dipole moments (dipx/y/z) 993!This routine is previously specific for ADCH, but can be extended to any other types of atomic charges 994!The "charge" is inputted Hirshfeld charge, finally it is replaced by ADC charge 995!chgtype 5= Becke with/without ADC, 6= ADCH 996subroutine doADC(dipx,dipy,dipz,charge,realdip,chgtype) 997use defvar 998use util 999implicit real*8 (a-h,o-z) 1000integer chgtype 1001real*8 gammamat(3,3),mat(3,3),avgr(3,1),avgrr(3,3),r(3,1),dip(3,1),tmp(1,1),eigval(3),eigvecmat(3,3) 1002real*8 w(ncenter),chargecorr(ncenter),charge(ncenter) 1003real*8 dipx(ncenter),dipy(ncenter),dipz(ncenter),realdip 1004 1005write(*,*) 1006write(*,*) "Now calculating atomic dipole moment corrected charge..." 1007write(*,*) 1008chargecorr=charge 1009 1010do i=1,ncenter 1011 if (ishowchgtrans==1) write(*,"('Atom: 'i4,a)") i,a(i)%name !ishowchgtrans==1 means output detail of charge transferation process during atomic dipole moment correction 1012 !Initialize variables 1013 totq=0.0D0 1014 tottmpdipx=0.0D0 1015 tottmpdipy=0.0D0 1016 tottmpdipz=0.0D0 1017 avgr=0.0D0 1018 avgrr=0.0D0 1019 dip(1,1)=dipx(i) 1020 dip(2,1)=dipy(i) 1021 dip(3,1)=dipz(i) 1022 1023 !Calculate weight of every atom 1024 do j=1,ncenter 1025 r(1,1)=a(j)%x-a(i)%x 1026 r(2,1)=a(j)%y-a(i)%y 1027 r(3,1)=a(j)%z-a(i)%z 1028 r2=r(1,1)**2+r(2,1)**2+r(3,1)**2 1029 distij=dsqrt(r2) 1030 1031 !Use modified Becke weight function with vdW radii criterion 1032 rmaxdist=vdwr(a(i)%index)+vdwr(a(j)%index) 1033 tr=distij/(rmaxdist/2.0D0)-1 !Transform variable so that it can in 0~rmaxdist range 1034 tr=1.5*(tr)-0.5*(tr)**3 1035 tr=1.5*(tr)-0.5*(tr)**3 1036 w(j)=0.5*(1-(1.5*tr-0.5*tr**3)) 1037 if (distij>rmaxdist) w(j)=0.0D0 1038 1039 avgr=avgr+w(j)*r 1040 avgrr=avgrr+w(j)*matmul(r,transpose(r)) 1041 end do 1042 1043 wtot=sum(w) 1044 avgr=avgr/wtot !Get <r> 1045 avgrr=avgrr/wtot !Get <rr+> 1046 gammamat=avgrr-matmul(avgr,transpose(avgr)) 1047 call Diagmat(gammamat,eigvecmat,eigval,500,1D-10) 1048 1049 rmaxv=maxval(eigval) 1050 mat=0.0D0 1051 tmpmin=1D-5 1052 mat(1,1)=1.0D0/(eigval(1)+tmpmin*(rmaxv+tmpmin)) 1053 mat(2,2)=1.0D0/(eigval(2)+tmpmin*(rmaxv+tmpmin)) 1054 mat(3,3)=1.0D0/(eigval(3)+tmpmin*(rmaxv+tmpmin)) 1055 1056 !Use transform matrix to transform r in old coordinate to r' in new coordinate, and transform P to P' 1057 !r=matmul(eigvecmat,r'), so r'=matmul(eigvecmat^(-1),r), because eigvecmat is unitary matrix, r'=matmul(transepose(eigvecmat),r) 1058 avgr=matmul(transpose(eigvecmat),avgr) 1059 dip=matmul(transpose(eigvecmat),dip) 1060 1061 !All values need have been calculated, now calculate final result 1062 do j=1,ncenter 1063 r(1,1)=a(j)%x-a(i)%x 1064 r(2,1)=a(j)%y-a(i)%y 1065 r(3,1)=a(j)%z-a(i)%z 1066 r=matmul(transpose(eigvecmat),r) ! Get r(i,j) vector in new coordinate 1067 tmp=w(j)/wtot*matmul(matmul(transpose(r-avgr),mat) ,dip) !delta q, namely the charge which atom i gives atom j 1068 chargecorr(j)=chargecorr(j)+tmp(1,1) !Charge after corrected 1069 if (ishowchgtrans==1) write(*,"(' Give atom ',i4,a4,f15.12,' Weight',2f15.12)") j,a(j)%name,tmp(1,1),w(j) 1070 totq=totq+tmp(1,1) 1071 tottmpdipx=tottmpdipx+(a(j)%x-a(i)%x)*tmp(1,1) 1072 tottmpdipy=tottmpdipy+(a(j)%y-a(i)%y)*tmp(1,1) 1073 tottmpdipz=tottmpdipz+(a(j)%z-a(i)%z)*tmp(1,1) 1074 end do 1075 if (ishowchgtrans==1) write(*,*) 1076end do 1077 1078write(*,*) " ======= Summary of atomic dipole moment corrected (ADC) charges =======" 1079do i=1,ncenter 1080 write(*,"(' Atom: ',i4,a,' Corrected charge:',f12.6,' Before:',f12.6)") i,a(i)%name,chargecorr(i),charge(i) 1081end do 1082write(*,"(' Summing up all corrected charges:',f12.7)") sum(chargecorr) 1083if (chgtype==5) write(*,"(a)") " Note: The values shown after ""Corrected charge"" are atomic dipole moment corrected Becke charges, the ones after ""Before"" are normal Becke charges" 1084if (chgtype==6) write(*,"(a)") " Note: The values shown after ""Corrected charge"" are ADCH charges, the ones after ""Before"" are Hirshfeld charges" 1085ADCdipx=sum(a%x*chargecorr) 1086ADCdipy=sum(a%y*chargecorr) 1087ADCdipz=sum(a%z*chargecorr) 1088ADCdip=sqrt(ADCdipx**2+ADCdipy**2+ADCdipz**2) 1089write(*,*) 1090write(*,"(' Total dipole from ADC charges (a.u.)',f11.7,' Error:',f11.7)") ADCdip,abs(ADCdip-realdip) 1091write(*,"(' X/Y/Z of dipole moment from the charge (a.u.)',3f11.7)") ADCdipx,ADCdipy,ADCdipz 1092charge=chargecorr !Overlay charge array, then return to Hirshfeld module and output result to .chg file 1093end subroutine 1094 1095 1096!!--------- Calculate CM5 charge based on Hirshfeld charge 1097subroutine doCM5(charge) 1098use defvar 1099implicit real*8 (a-h,o-z) 1100real*8 charge(ncenter),CMcharge(ncenter),radius(118),Dparm(118) 1101alpha=2.474D0 1102Dparm=0D0 1103Dparm(1)=0.0056D0 1104Dparm(2)=-0.1543D0 1105Dparm(4)=0.0333D0 1106Dparm(5)=-0.1030D0 1107Dparm(6)=-0.0446D0 1108Dparm(7)=-0.1072D0 1109Dparm(8)=-0.0802D0 1110Dparm(9)=-0.0629D0 1111Dparm(10)=-0.1088D0 1112Dparm(11)=0.0184D0 1113Dparm(13)=-0.0726D0 1114Dparm(14)=-0.0790D0 1115Dparm(15)=-0.0756D0 1116Dparm(16)=-0.0565D0 1117Dparm(17)=-0.0444D0 1118Dparm(18)=-0.0767D0 1119Dparm(19)=0.0130D0 1120Dparm(31)=-0.0512D0 1121Dparm(32)=-0.0557D0 1122Dparm(33)=-0.0533D0 1123Dparm(34)=-0.0399D0 1124Dparm(35)=-0.0313D0 1125Dparm(36)=-0.0541D0 1126Dparm(37)=0.0092D0 1127Dparm(49)=-0.0361D0 1128Dparm(50)=-0.0393D0 1129Dparm(51)=-0.0376D0 1130Dparm(52)=-0.0281D0 1131Dparm(53)=-0.0220D0 1132Dparm(54)=-0.0381D0 1133Dparm(55)=0.0065D0 1134Dparm(81)=-0.0255D0 1135Dparm(82)=-0.0277D0 1136Dparm(83)=-0.0265D0 1137Dparm(84)=-0.0198D0 1138Dparm(85)=-0.0155D0 1139Dparm(86)=-0.0269D0 1140Dparm(87)=0.0046D0 1141Dparm(113)=-0.0179D0 1142Dparm(114)=-0.0195D0 1143Dparm(115)=-0.0187D0 1144Dparm(116)=-0.0140D0 1145Dparm(117)=-0.0110D0 1146Dparm(118)=-0.0189D0 1147!As shown in CM5 paper, the covalent radii used in CM5 equation are tabulated in CRC book 91th, where they are obtained as follows: 1148!For Z=1~96, the radii are the average of CSD radii (For Fe, Mn, Co the low-spin is used) and Pyykko radii 1149!For Z=97~118, the radii are Pyykko radii 1150radius(1:96)=(covr(1:96)+covr_pyy(1:96))/2D0 1151radius(97:118)=covr_pyy(97:118) 1152radius=radius*b2a !Because the radii have already been converted to Bohr, so we convert them back to Angstrom 1153 1154if (ishowchgtrans==1) write(*,"(/,a)") " Details of CM5 charge correction:" 1155 1156do iatm=1,ncenter 1157 CMcorr=0 1158 iZ=a(iatm)%index 1159 do jatm=1,ncenter 1160 if (iatm==jatm) cycle 1161 jZ=a(jatm)%index 1162 Bval=exp( -alpha*(distmat(iatm,jatm)*b2a-radius(iZ)-radius(jZ)) ) 1163 if (iZ==1.and.jZ==6) then 1164 Tval=0.0502D0 1165 else if (iZ==6.and.jZ==1) then 1166 Tval=-0.0502D0 1167 else if (iZ==1.and.jZ==7) then 1168 Tval=0.1747D0 1169 else if (iZ==7.and.jZ==1) then 1170 Tval=-0.1747D0 1171 else if (iZ==1.and.jZ==8) then 1172 Tval=0.1671D0 1173 else if (iZ==8.and.jZ==1) then 1174 Tval=-0.1671D0 1175 else if (iZ==6.and.jZ==7) then 1176 Tval=0.0556D0 1177 else if (iZ==7.and.jZ==6) then 1178 Tval=-0.0556D0 1179 else if (iZ==6.and.jZ==8) then 1180 Tval=0.0234D0 1181 else if (iZ==8.and.jZ==6) then 1182 Tval=-0.0234D0 1183 else if (iZ==7.and.jZ==8) then 1184 Tval=-0.0346D0 1185 else if (iZ==8.and.jZ==7) then 1186 Tval=0.0346D0 1187 else 1188 Tval=Dparm(iZ)-Dparm(jZ) 1189 end if 1190 CMcorr=CMcorr+Tval*Bval 1191 if (ishowchgtrans==1) then 1192 write(*,"(i4,a,i4,a,' B_term:',f10.5,' T_term:',f10.5,' Corr. charge:',f10.5)") iatm,a(iatm)%name,jatm,a(jatm)%name,Bval,Tval,Tval*Bval 1193 end if 1194 end do 1195 CMcharge(iatm)=charge(iatm)+CMcorr 1196end do 1197write(*,*) 1198write(*,*) " ======= Summary of CM5 charges =======" 1199do i=1,ncenter 1200 write(*,"(' Atom: ',i4,a,' CM5 charge:',f12.6,' Hirshfeld charge:',f12.6)") i,a(i)%name,CMcharge(i),charge(i) 1201end do 1202write(*,"(' Summing up all CM5 charges:',f15.8)") sum(CMcharge) 1203CM5dipx=sum(a%x*CMcharge) 1204CM5dipy=sum(a%y*CMcharge) 1205CM5dipz=sum(a%z*CMcharge) 1206CM5dip=sqrt(CM5dipx**2+CM5dipy**2+CM5dipz**2) 1207write(*,*) 1208write(*,"(' Total dipole moment from CM5 charges',f12.7,' a.u.')") CM5dip 1209write(*,"(' X/Y/Z of dipole moment from CM5 charges',3f10.5, ' a.u.')") CM5dipx,CM5dipy,CM5dipz 1210charge=CMcharge 1211end subroutine 1212 1213 1214 1215!!------------ Calculate charges by fitting ESP, currently CHELPG grid and MK grid are used 1216!! TrEsp (transition charge from electrostatic potential) can also be calculated, see manual 1217!itype=1:CHELPG itype=2:MK 1218subroutine fitESP(itype) 1219use util 1220use defvar 1221use function 1222implicit real*8 (a-h,o-z) 1223character*200 addcenfile,extptfile 1224character selectyn,chgfilename*200 1225integer itype 1226integer :: nlayer=4 !Number of layers of points for MK 1227integer :: inuctype=1 1228real*8 :: espfitvdwr(0:nelesupp)=-1D0,sclvdwlayer(100)=(/1.4D0,1.6D0,1.8D0,2.0D0,(0D0,i=5,100)/) 1229real*8,allocatable :: ESPptval(:),ESPptx(:),ESPpty(:),ESPptz(:),Bmat(:),Amat(:,:),Amatinv(:,:),atmchg(:) 1230real*8,allocatable :: fitcenx(:),fitceny(:),fitcenz(:),fitcenvdwr(:),disptcen(:),origsphpt(:,:) 1231 1232! Read ESP and coordinates of fitting points from Gaussian Iop(6/33=2) output. Corresponding wavefunction file must be loaded to provide atom coordinates 1233! naddcen=0 1234! open(10,file="C:\gtest\benzene.out",status="old") 1235! call loclabel(10,"NAtoms") 1236! read(10,"(8x,i6)") nfitcen 1237! allocate(fitcenx(nfitcen),fitceny(nfitcen),fitcenz(nfitcen)) 1238! call loclabel(10,"Atomic Center 1 is at") 1239! do i=1,nfitcen 1240! read(10,"(32x,3f10.6)") fitcenx(i),fitceny(i),fitcenz(i) 1241! end do 1242! fitcenx=fitcenx/b2a 1243! fitceny=fitceny/b2a 1244! fitcenz=fitcenz/b2a 1245! call loclabel(10,"points will be used for",ifound) 1246! read(10,*) nESPpt 1247! write(*,"('Number of fitting points used:',i10)") nESPpt 1248! allocate(ESPptval(nESPpt),ESPptx(nESPpt),ESPpty(nESPpt),ESPptz(nESPpt)) 1249! matdim=nfitcen+1 1250! call loclabel(10,"ESP Fit Center",ifound) 1251! if (ifound==0) write(*,*) "Cannot locate ""ESP Fit Center"" field" 1252! do ipt=1,nESPpt 1253! read(10,"(32x,3f10.6)") ESPptx(ipt),ESPpty(ipt),ESPptz(ipt) 1254! end do 1255! ESPptx=ESPptx/b2a !Convert to Bohr unit 1256! ESPpty=ESPpty/b2a 1257! ESPptz=ESPptz/b2a 1258! call loclabel(10," Fit ",ifound,0) 1259! if (ifound==0) write(*,*) "Cannot locate "" Fit "" field" 1260! do ipt=1,nESPpt 1261! read(10,"(14x,f10.6)") ESPptval(ipt) 1262! end do 1263! ! goto 332 1264! goto 333 1265 1266fitspc=0.3D0/b2a !Spacing between grid for CHELPG 1267extdis=2.8D0/b2a !Extend 2.8 Angstrom to each side for CHELPG 1268densperarea=6D0*b2a**2 !Point density per Angstrom**2 for MK, in order to convert to Bohr**2, multiply by b2a**2 1269 1270iaddcen=0 !If give Additional center 1271iuseextpt=0 !If use external points 1272iskipespcalc=0 !If read ESP from external file directly rather than calculate here 1273 1274do while(.true.) 1275 write(*,*) 1276 write(*,*) "Note: All units in this module are in a.u." 1277 write(*,*) "-2 Load additional fitting centers from external file" 1278 if (iuseextpt==0) write(*,*) "-1 Use fitting points recorded in external file instead of generating them" 1279 write(*,*) "0 Return" 1280 write(*,*) "1 Start calculation!" 1281 if (iuseextpt==0) then 1282 if (itype==1) then 1283 write(*,"(' 2 Set grid spacing, current:',f7.3,' Bohr (',f7.3,' Angstrom)')") fitspc,fitspc*b2a 1284 write(*,"(' 3 Set box extension, current:',f7.3,' Bohr (',f7.3,' Angstrom)')") extdis,extdis*b2a 1285 else if (itype==2) then 1286 write(*,"(' 2 Set number of points per Angstrom^2, current:',f10.3)") densperarea/b2a**2 !Temporary convert to Angstrom**2 for convention 1287 write(*,"(' 3 Set number of layers per atom, current:',i4)") nlayer 1288 write(*,"(' 4 Set the value times van der Waals radius in each layer')") 1289 end if 1290 end if 1291 if (iuseextpt==0) then 1292 if (inuctype==1) write(*,*) "5 Switch if taking nuclear charge into account, current: Yes" 1293 if (inuctype==2) write(*,*) "5 Switch if taking nuclear charge into account, current: No" 1294 end if 1295 read(*,*) isel 1296 1297 if (isel==-2) then 1298 iaddcen=1 1299 write(*,*) "Input the name of the file recording coordinates of additional fitting centers" 1300 read(*,"(a)") addcenfile 1301 write(*,*) "Done!" 1302 else if (isel==-1) then 1303 iuseextpt=1 1304 write(*,*) "Input the name of the file recording coordinates of ESP fitting points" 1305 read(*,"(a)") extptfile 1306 write(*,*) "OK, the points recorded in this file will be used as fitting points" 1307 else if (isel==0) then 1308 Return 1309 else if (isel==1) then 1310 exit 1311 else if (isel==2) then 1312 write(*,*) "Input new value" 1313 if (itype==1) read(*,*) fitspc 1314 if (itype==2) then 1315 read(*,*) densperarea 1316 densperarea=densperarea*b2a**2 1317 end if 1318 else if (isel==3) then 1319 write(*,*) "Input new value" 1320 if (itype==1) read(*,*) extdis 1321 if (itype==2) read(*,*) nlayer 1322 else if (isel==4.and.itype==2) then 1323 write(*,*) "Current values:" 1324 do ilayer=1,nlayer 1325 write(*,"(' Layer',i4,' :',f8.4)") ilayer,sclvdwlayer(ilayer) 1326 end do 1327 write(*,*) 1328 do ilayer=1,nlayer 1329 write(*,"(a,i4,', e.g. 1.5')") " Input value for layer",ilayer 1330 read(*,*) sclvdwlayer(ilayer) 1331 end do 1332 else if (isel==5) then 1333 if (inuctype==1) then 1334 inuctype=2 1335 else if (inuctype==2) then 1336 inuctype=1 1337 end if 1338 end if 1339end do 1340 1341!Set vdW radius 1342if (itype==1) then !For CHELPG 1343 espfitvdwr(1:2)=1.45D0 !vdW radius copied from GetvdW routine (utilam), some of them are given in CHELPG original paper 1344 espfitvdwr(3:6)=1.5D0 1345 espfitvdwr(7:10)=1.7D0 1346 espfitvdwr(11:18)=2D0 1347 espfitvdwr(1:18)=espfitvdwr(1:18)/b2a !to Bohr 1348else if (itype==2) then !For MK, copied from GetvdW routine (utilam) 1349 espfitvdwr(1:17)=(/1.20d0,1.20d0,1.37d0,1.45d0,1.45d0,1.50d0,1.50d0,& 1350 1.40d0,1.35d0,1.30d0,1.57d0,1.36d0,1.24d0,1.17d0,1.80d0,1.75d0,1.70d0/) 1351 espfitvdwr(1:17)=espfitvdwr(1:17)/b2a 1352end if 1353write(*,*) "Atomic radii used:" 1354do ielem=1,nelesupp 1355 if (any(a%index==ielem).and.espfitvdwr(ielem)/=-1D0) write(*,"(' Element:',a,' vdW radius (Angstrom):',f6.3)") ind2name(ielem),espfitvdwr(ielem)*b2a 1356end do 1357 1358!Check sanity and complete vdW radius table for all involved elements 1359do iatm=1,ncenter 1360 if (espfitvdwr(a(iatm)%index)==-1D0) then 1361 write(*,"(' vdW radius used in fitting for element ',a,' is missing, input the radius (Bohr)')") ind2name(a(iatm)%index) 1362 write(*,"(a)") " Hint: If you don't know how to deal with the problem, simply input 3.4. (However, the radius of 3.4 Bohr may be not very appropriate for current element)" 1363 read(*,*) espfitvdwr(a(iatm)%index) 1364 end if 1365end do 1366 1367!Check total number of fitting centers 1368naddcen=0 1369if (iaddcen==1) then 1370 open(10,file=addcenfile,status="old") 1371 read(10,*) naddcen 1372end if 1373nfitcen=ncenter+naddcen 1374allocate(fitcenx(nfitcen),fitceny(nfitcen),fitcenz(nfitcen),fitcenvdwr(nfitcen),disptcen(nfitcen)) 1375 1376!Generate information of fitting centers 1377do iatm=1,ncenter 1378 fitcenx(iatm)=a(iatm)%x 1379 fitceny(iatm)=a(iatm)%y 1380 fitcenz(iatm)=a(iatm)%z 1381 fitcenvdwr(iatm)=espfitvdwr(a(iatm)%index) !vdW radius for each fitting center 1382end do 1383if (iaddcen==1) then 1384 do icen=ncenter+1,ncenter+naddcen 1385 read(10,*) fitcenx(icen),fitceny(icen),fitcenz(icen) 1386 fitcenvdwr(icen)=0D0 1387 end do 1388 close(10) 1389end if 1390 1391write(*,*) 1392if (iuseextpt==0) then !Count number and generate coordinates of fitting points 1393 if (itype==1) then !CHELPG 1394 xlow=minval(fitcenx)-extdis 1395 xhigh=maxval(fitcenx)+extdis 1396 ylow=minval(fitceny)-extdis 1397 yhigh=maxval(fitceny)+extdis 1398 zlow=minval(fitcenz)-extdis 1399 zhigh=maxval(fitcenz)+extdis 1400 xlen=xhigh-xlow 1401 ylen=yhigh-ylow 1402 zlen=zhigh-zlow 1403 nxfit=int(xlen/fitspc)+1 1404 nyfit=int(ylen/fitspc)+1 1405 nzfit=int(zlen/fitspc)+1 1406 nESPpt=0 1407 do ix=0,nxfit 1408 do iy=0,nyfit 1409 do iz=0,nzfit 1410 tmpx=xlow+ix*fitspc 1411 tmpy=ylow+iy*fitspc 1412 tmpz=zlow+iz*fitspc 1413 do icen=1,nfitcen 1414 disptcen(icen)=dsqrt( (fitcenx(icen)-tmpx)**2+(fitceny(icen)-tmpy)**2+(fitcenz(icen)-tmpz)**2 ) 1415 if (disptcen(icen)<=fitcenvdwr(icen)) exit 1416 if (icen==nfitcen.and.any(disptcen<=extdis)) nESPpt=nESPpt+1 1417 end do 1418 end do 1419 end do 1420 end do 1421 write(*,"(' Number of fitting points used:',i10)") nESPpt 1422 allocate(ESPptval(nESPpt),ESPptx(nESPpt),ESPpty(nESPpt),ESPptz(nESPpt)) 1423 iESPpt=0 1424 do ix=0,nxfit 1425 do iy=0,nyfit 1426 do iz=0,nzfit 1427 tmpx=xlow+ix*fitspc 1428 tmpy=ylow+iy*fitspc 1429 tmpz=zlow+iz*fitspc 1430 do icen=1,nfitcen 1431 disptcen(icen)=dsqrt( (fitcenx(icen)-tmpx)**2+(fitceny(icen)-tmpy)**2+(fitcenz(icen)-tmpz)**2 ) 1432 if (disptcen(icen)<=fitcenvdwr(icen)) exit 1433 if (icen==nfitcen.and.any(disptcen<=extdis)) then 1434 iESPpt=iESPpt+1 1435 ESPptx(iESPpt)=tmpx 1436 ESPpty(iESPpt)=tmpy 1437 ESPptz(iESPpt)=tmpz 1438 end if 1439 end do 1440 end do 1441 end do 1442 end do 1443 1444 else if (itype==2) then !MK, don't reserve spatial space for additional center 1445 cutinnerscl=minval(sclvdwlayer(1:nlayer)) 1446 write(*,"(' Note: If distance between a ESP point and any atom is smaller than',f6.3,' multiplied by corresponding vdW radius, then the point will be discarded')") cutinnerscl 1447 nESPpt=0 1448 maxsphpt=nint(4D0*pi*(maxval(fitcenvdwr)*maxval(sclvdwlayer))**2 *densperarea) !Find maximal possible number of points in unit sphere to allocate temporary origsphpt 1449 allocate(origsphpt(3,maxsphpt)) 1450 do icen=1,ncenter !Rather than nfitcen. Count how many possible ESP points in total 1451 do ilayer=1,nlayer 1452 numsphpt=nint(4D0*pi*(fitcenvdwr(icen)*sclvdwlayer(ilayer))**2 *densperarea) 1453 nESPpt=nESPpt+numsphpt 1454 end do 1455 end do 1456 allocate(ESPptval(nESPpt),ESPptx(nESPpt),ESPpty(nESPpt),ESPptz(nESPpt)) !Currently nESPpt is upper limit 1457 iESPpt=0 1458 do icen=1,ncenter 1459 do ilayer=1,nlayer 1460 radius=fitcenvdwr(icen)*sclvdwlayer(ilayer) 1461 numsphpt=nint(4D0*pi*radius**2 *densperarea) 1462 call unitspherept(origsphpt,numsphpt) !Input expected number of point in unit sphere, return actual number of points 1463 origsphpt(:,1:numsphpt)=origsphpt(:,1:numsphpt)*radius 1464 origsphpt(1,1:numsphpt)=origsphpt(1,1:numsphpt)+fitcenx(icen) !Move unit sphere to atomic center 1465 origsphpt(2,1:numsphpt)=origsphpt(2,1:numsphpt)+fitceny(icen) 1466 origsphpt(3,1:numsphpt)=origsphpt(3,1:numsphpt)+fitcenz(icen) 1467 do ipt=1,numsphpt 1468 tmpx=origsphpt(1,ipt) 1469 tmpy=origsphpt(2,ipt) 1470 tmpz=origsphpt(3,ipt) 1471 iok=1 1472 do icen2=1,ncenter 1473 if (icen2==icen) cycle 1474 disptcensq=(fitcenx(icen2)-tmpx)**2+(fitceny(icen2)-tmpy)**2+(fitcenz(icen2)-tmpz)**2 !distance between point and center 1475 if (disptcensq<(fitcenvdwr(icen2)*cutinnerscl)**2) then !Less than vdW RADIUS*cutinner of atom icen2, it should be ommitted 1476 iok=0 1477 exit 1478 end if 1479 end do 1480 if (iok==1) then 1481 iESPpt=iESPpt+1 1482 ESPptx(iESPpt)=tmpx 1483 ESPpty(iESPpt)=tmpy 1484 ESPptz(iESPpt)=tmpz 1485 end if 1486 end do 1487 end do 1488 end do 1489 nESPpt=iESPpt 1490 deallocate(origsphpt) 1491 write(*,"(' Number of fitting points used:',i10)") nESPpt 1492 end if 1493 1494else if (iuseextpt==1) then !Directly use external fitting points 1495 open(10,file=extptfile,status="old") 1496 read(10,*) nESPpt 1497 if (nESPpt<0) then 1498 iskipespcalc=1 !If the number of fitting points is negative, that means the fourth column records ESP value and needn't to be recalculated 1499 write(*,*) "ESP value of all fitting points are read from external file directly" 1500 end if 1501 nESPpt=abs(nESPpt) 1502 write(*,"(' Number of fitting points used:',i10)") nESPpt 1503 allocate(ESPptval(nESPpt),ESPptx(nESPpt),ESPpty(nESPpt),ESPptz(nESPpt)) 1504 do i=1,nESPpt 1505 if (iskipespcalc==0) read(10,*) ESPptx(i),ESPpty(i),ESPptz(i) 1506 if (iskipespcalc==1) read(10,*) ESPptx(i),ESPpty(i),ESPptz(i),ESPptval(i) 1507 end do 1508 close(10) 1509end if 1510 1511332 continue 1512!Generate ESP value of fitting points 1513if (iskipespcalc==0) then 1514 write(*,*) "Calculating ESP, please wait..." 1515 itmp=1 1516 do ipt=1,nESPpt 1517 if (ipt>=itmp*300) then 1518 write(*,"(' Finished:',i10,' /',i10)") ipt,nESPpt 1519 itmp=itmp+1 1520 end if 1521 if (inuctype==1) then 1522 ESPptval(ipt)=totesp((ESPptx(ipt)),(ESPpty(ipt)),(ESPptz(ipt))) 1523 else if (inuctype==2) then !Don't consider nuclear charge contribution 1524 ESPptval(ipt)=eleesp((ESPptx(ipt)),(ESPpty(ipt)),(ESPptz(ipt))) 1525 end if 1526 end do 1527 write(*,*) "Done!" 1528end if 1529 1530333 continue !We just read ESP points from Gaussian output directly, so jump to here 1531matdim=nfitcen+1 1532allocate(Bmat(matdim),Amat(matdim,matdim),Amatinv(matdim,matdim),atmchg(matdim)) 1533!See original paper of MK for detail of algorithem 1534!Forming Amat 1535Amat=0D0 1536do icen=1,nfitcen 1537 do jcen=icen,nfitcen 1538 do ipt=1,nESPpt 1539 dis1=dsqrt( (ESPptx(ipt)-fitcenx(icen))**2 + (ESPpty(ipt)-fitceny(icen))**2 + (ESPptz(ipt)-fitcenz(icen))**2 ) 1540 dis2=dsqrt( (ESPptx(ipt)-fitcenx(jcen))**2 + (ESPpty(ipt)-fitceny(jcen))**2 + (ESPptz(ipt)-fitcenz(jcen))**2 ) 1541 Amat(icen,jcen)=Amat(icen,jcen)+1D0/dis1/dis2 1542 end do 1543 end do 1544end do 1545Amat=Amat+transpose(Amat) 1546do i=1,nfitcen 1547 Amat(i,i)=Amat(i,i)/2D0 1548end do 1549Amat(matdim,:)=1D0 1550Amat(:,matdim)=1D0 1551Amat(matdim,matdim)=0D0 1552!Forming Bmat 1553Bmat=0D0 1554do icen=1,nfitcen 1555 do ipt=1,nESPpt 1556 dis=dsqrt( (ESPptx(ipt)-fitcenx(icen))**2 + (ESPpty(ipt)-fitceny(icen))**2 + (ESPptz(ipt)-fitcenz(icen))**2 ) 1557 Bmat(icen)=Bmat(icen)+ESPptval(ipt)/dis 1558 end do 1559end do 1560if (inuctype==1) then 1561 Bmat(matdim)=sum(a(:)%charge)-nelec !Net charge of the system 1562else if (inuctype==2) then 1563 Bmat(matdim)=-nelec 1564end if 1565Amatinv=invmat(Amat,matdim) 1566atmchg=matmul(Amatinv,Bmat) 1567 1568!Output summary 1569write(*,*) " Center X Y Z Charge" 1570do i=1,ncenter 1571 write(*,"(i6,a,3f12.6,f16.6)") i,ind2name(a(i)%index),fitcenx(i),fitceny(i),fitcenz(i),atmchg(i) 1572end do 1573do i=ncenter+1,ncenter+naddcen 1574 write(*,"(i6,2x,3f12.6,f16.6)") i,fitcenx(i),fitceny(i),fitcenz(i),atmchg(i) 1575end do 1576write(*,"(' Sum of charges:',f12.6)") sum(atmchg(1:nfitcen)) 1577!Calculate RMSE and RRMSE 1578RMSE=0D0 1579do ipt=1,nESPpt 1580 atmchgesp=0D0 1581 do icen=1,nfitcen 1582 dis=dsqrt( (ESPptx(ipt)-fitcenx(icen))**2 + (ESPpty(ipt)-fitceny(icen))**2 + (ESPptz(ipt)-fitcenz(icen))**2 ) 1583 atmchgesp=atmchgesp+atmchg(icen)/dis 1584 end do 1585 RMSE=RMSE+(ESPptval(ipt)-atmchgesp)**2 1586end do 1587RRMSE=dsqrt(RMSE/sum(ESPptval(1:nESPpt)**2)) 1588RMSE=dsqrt(RMSE/nESPpt) 1589write(*,"(' RMSE:',f12.6,' RRMSE:',f12.6)") RMSE,RRMSE 1590 1591!Show fragment charge 1592if (allocated(frag1)) write(*,"(/,' Fragment charge:',f12.6)") sum(atmchg(frag1)) 1593 1594write(*,*) 1595write(*,"(a)") " If output coordinates and ESP value of all fitting points to ESPfitpt.txt in current folder? (y/n)" 1596read(*,*) selectyn 1597if (selectyn=='y'.or.selectyn=="Y") then 1598 open(10,file="ESPfitpt.txt",status="replace") 1599 write(10,*) nESPpt 1600 do ipt=1,nESPpt 1601 write(10,"(3f12.6,f14.8)") ESPptx(ipt),ESPpty(ipt),ESPptz(ipt),ESPptval(ipt) 1602 end do 1603 write(*,*) "Data have been outputted to ESPfitpt.txt in current folder" 1604 write(*,"(a)") " All units are in a.u. The first line shows the number of fitting points, the first three columns are X,Y,Z coordinates, the last column corresponds to ESP value" 1605 close(10) 1606end if 1607 1608call path2filename(firstfilename,chgfilename) 1609write(*,"(a)") " If output atom coordinates with charges to "//trim(chgfilename)//".chg file in current folder? (y/n)" 1610read(*,*) selectyn 1611if (selectyn=='y'.or.selectyn=="Y") then 1612 open(10,file=trim(chgfilename)//".chg",status="replace") 1613 do i=1,ncenter 1614 write(10,"(a4,4f12.6)") a(i)%name,a(i)%x*b2a,a(i)%y*b2a,a(i)%z*b2a,atmchg(i) 1615 end do 1616 do i=ncenter+1,ncenter+naddcen 1617 write(10,"(a4,4f12.6)") " Bq",a(i)%x*b2a,a(i)%y*b2a,a(i)%z*b2a,atmchg(i) 1618 end do 1619 close(10) 1620 write(*,"(a)") " Result have been saved to "//trim(chgfilename)//".chg in current folder" 1621 write(*,"(a)") " Columns from 1 to 5 are name,X,Y,Z,charge respectively, unit is Angstrom" 1622end if 1623 1624!Output fitting points to pdb file for visualizing, only for debug purpose 1625! open(10,file="ESPfitpt.pdb",status="replace") 1626! do ipt=1,nESPpt 1627! write(10,"(a6,i5,1x,a4,1x,a3, 1x,a1,i4,4x,3f8.3,2f6.2,10x,a2)") & 1628! "HETATM",i,' '//"O "//' ',"MOL",'A',1,ESPptx(ipt),ESPpty(ipt),ESPptz(ipt),1.0,ESPptval(ipt)*au2kcal,"O " 1629! end do 1630! write(*,*) "Data have been outputted to ESPfitpt.pdb in current folder" 1631! close(10) 1632 1633end subroutine 1634 1635 1636!!!------------- Generate numpt points scattered evenly in an unit sphere 1637! Input argument numpt is the expected number of points, while the return value is actual number 1638! ptcrd store coordinates of the points 1639subroutine unitspherept(ptcrd,numpt) 1640implicit real*8 (a-h,o-z) 1641real*8 ptcrd(3,numpt) 1642integer numpt 1643pi=3.141592653589793D0 1644!The average number of equator points in all XY layes is numequ*2/pi, and there are numvert=numequ/2 layers 1645!Solve (numequ*2/pi)*numequ/2=numpt one can get numequ=sqrt(numpt*pi) 1646numequ=int(sqrt(numpt*pi)) !Maximal number of point in each XY layer 1647numvert=numequ/2 1648ipt=0 1649do ivert=0,numvert 1650 angz=dfloat(ivert)/numvert*pi 1651 scalexy=sin(angz) 1652 z=cos(angz) 1653 numxy=int(numequ*scalexy) 1654 if (numxy==0) numxy=1 1655 do ihori=1,numxy 1656 ipt=ipt+1 1657 if (ipt>numpt) then 1658 numpt=ipt-1 1659 return 1660 end if 1661 angxy=2D0*pi*ihori/numxy 1662 ptcrd(1,ipt)=cos(angxy)*scalexy 1663 ptcrd(2,ipt)=sin(angxy)*scalexy 1664 ptcrd(3,ipt)=z 1665 end do 1666end do 1667numpt=ipt 1668end subroutine 1669 1670 1671 1672 1673 1674 1675 1676 1677!!============================ Hirshfeld-I ============================!! 1678!!============================ Hirshfeld-I ============================!! 1679!!============================ Hirshfeld-I ============================!! 1680!!============================ Hirshfeld-I ============================!! 1681!!============================ Hirshfeld-I ============================!! 1682!Wrapper of Hirshfeld-I module to automatically set radpot and sphpot to proper values 1683!itype=1: Normal population analysis =2: Only used to generate proper atomic space (i.e. Don't do unnecessary things) 1684subroutine Hirshfeld_I_wrapper(itype) 1685use defvar 1686implicit real*8 (a-h,o-z) 1687radpotold=radpot 1688sphpotold=sphpot 1689if (iautointgrid==1) then 1690 radpot=30 1691 sphpot=170 1692 if (any(a%index>18)) radpot=40 1693 if (any(a%index>36)) radpot=50 1694 if (any(a%index>54)) radpot=60 1695end if 1696call Hirshfeld_I(itype) 1697if (iautointgrid==1) then 1698 radpot=radpotold 1699 sphpot=sphpotold 1700end if 1701end subroutine 1702 1703!!--------- Calculate Hirshfeld-I charge and yield final atomic radial density 1704!I've compared this module with hipart, this module is faster than hipart, and the accuracy under default setting is at least never lower than hipart 1705subroutine Hirshfeld_I(itype) 1706use defvar 1707use function 1708use util 1709implicit real*8 (a-h,o-z) 1710integer itype 1711type(content) gridatm(radpot*sphpot),gridatmorg(radpot*sphpot) 1712real*8 molrho(radpot*sphpot),promol(radpot*sphpot),tmpdens(radpot*sphpot),selfdens(radpot*sphpot),molrhoall(ncenter,radpot*sphpot) 1713real*8 charge(ncenter),lastcharge(ncenter) !Atomic charge of current iter. and last iter. 1714real*8 radrholow(200),radrhohigh(200) 1715character sep,c80tmp*80,chgfilename*200,selectyn 1716character*2 :: statname(-4:4)=(/ "-4","-3","-2","-1","_0","+1","+2","+3","+4" /) 1717integer :: maxcyc=50,ioutmedchg=0 1718real*8 :: crit=0.0002D0 1719real*8,external :: fdens_rad 1720!Used for mode 2. e.g. atmstatgrid(iatm,igrid,jatm,-1) means density of jatm with -1 charge state at igrid around iatm 1721real*8,allocatable :: atmstatgrid(:,:,:,:) 1722ntotpot=radpot*sphpot 1723 1724!Mode 1 use very low memory but expensive, because most data is computed every iteration 1725!Mode 2 use large memory but fast, because most data is only computed once at initial stage 1726!The result of the two modes differ with each other marginally, probably because in mode 1 radial density is related to max(npthigh,nptlow), which is not involved in mode 2 1727!In principle, result of mode 2 is slightly better 1728imode=2 1729 1730!Ignore jatm contribution to iatm centered grids if distance between iatm and jatm is larger than 1.5 times of sum of their vdwr 1731!This can decrease lots of time for large system, the lose of accuracy can be ignored (error is ~0.0001 per atom) 1732ignorefar=1 1733 1734write(*,*) 1735if (itype==1) write(*,*) " =============== Iterative Hirshfeld (Hirshfeld-I) ===============" 1736if (itype==2) write(*,*) " ============== Generate Hirshfeld-I atomic weights ==============" 1737do while(.true.) 1738 if (imode==1) write(*,*) "-2 Switch algorithm, current: Slow & low memory requirement" 1739 if (imode==2) write(*,*) "-2 Switch algorithm, current: Fast & large memory requirement" 1740 if (itype==1) then 1741 if (ioutmedchg==0) write(*,*) "-1 Switch if output intermediate results, current: No" 1742 if (ioutmedchg==1) write(*,*) "-1 Switch if output intermediate results, current: Yes" 1743 write(*,*) "0 Return" 1744 end if 1745 write(*,*) "1 Start calculation!" 1746 write(*,"(a,i4)") " 2 Set the maximum number of iterations, current:",maxcyc 1747 write(*,"(a,f10.6)") " 3 Set convergence criterion of atomic charges, current:",crit 1748 read(*,*) isel 1749 if (isel==-2) then 1750 if (imode==1) then 1751 imode=2 1752 else 1753 imode=1 1754 crit=0.001 !mode 1 is more time-consuming, use loose criterion 1755 end if 1756 else if (isel==-1) then 1757 if (ioutmedchg==1) then 1758 ioutmedchg=0 1759 else 1760 ioutmedchg=1 1761 end if 1762 else if (isel==0) then 1763 return 1764 else if (isel==1) then 1765 exit 1766 else if (isel==2) then 1767 write(*,*) "Input maximum number of iterations, e.g. 30" 1768 read(*,*) maxcyc 1769 else if (isel==3) then 1770 write(*,*) "Input convergence criterion of atomic charges, e.g. 0.001" 1771 read(*,*) crit 1772 end if 1773end do 1774 1775!Generate all needed .rad files 1776call genatmradfile 1777 1778!====== Start calculation ======! 1779call walltime(iwalltime1) 1780CALL CPU_TIME(time_begin) 1781 1782!Currently all atoms share the same radial points 1783nradpt=200 1784itmp=0 1785do irad=nradpt,1,-1 1786 radx=cos(irad*pi/(nradpt+1)) 1787 itmp=itmp+1 1788 atmradpos(itmp)=(1+radx)/(1-radx) 1789end do 1790 1791!Generate single center integration grid 1792call gen1cintgrid(gridatmorg,iradcut) 1793write(*,*) 1794write(*,"(' Radial grids:',i4,' Angular grids:',i5,' Total:',i7,' After pruning:',i7)") radpot,sphpot,radpot*sphpot,radpot*sphpot-iradcut*sphpot 1795 1796!Calculate molecular density 1797write(*,*) "Calculating density of actual molecule for all grids..." 1798nthreads=getNThreads() 1799!$OMP parallel do shared(molrhoall) private(iatm,ipt,gridatm) num_threads(nthreads) 1800do iatm=1,ncenter 1801 gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule 1802 gridatm%y=gridatmorg%y+a(iatm)%y 1803 gridatm%z=gridatmorg%z+a(iatm)%z 1804 do ipt=1+iradcut*sphpot,ntotpot 1805 molrhoall(iatm,ipt)=fdens(gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z) 1806 end do 1807end do 1808!$OMP end parallel do 1809 1810if (allocated(atmradnpt)) deallocate(atmradnpt) 1811if (allocated(atmradrho)) deallocate(atmradrho) 1812allocate(atmradnpt(ncenter),atmradrho(ncenter,200)) 1813sep='/' !Separation symbol of directory 1814if (isys==1) sep='\' 1815 1816!Calculate contribution of all atoms in every state to each atomic centered grids 1817if (imode==2) then 1818 allocate(atmstatgrid(ncenter,ntotpot,ncenter,-2:2)) 1819 atmstatgrid=0 1820 write(*,*) "Calculating atomic density contribution to grids..." 1821 do iatm=1,ncenter !The center of grids 1822 write(*,"(' Progress:',i5,' /',i5)") iatm,ncenter 1823 gridatm%value=gridatmorg%value !Weight in this grid point 1824 gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule 1825 gridatm%y=gridatmorg%y+a(iatm)%y 1826 gridatm%z=gridatmorg%z+a(iatm)%z 1827 do istat=-2,2 !Charge state 1828 do jatm=1,ncenter 1829 if (ignorefar==1) then 1830 atmdist=dsqrt( (a(iatm)%x-a(jatm)%x)**2+(a(iatm)%y-a(jatm)%y)**2+(a(iatm)%z-a(jatm)%z)**2 ) 1831 if (atmdist>(vdwr(iatm)+vdwr(jatm))*1.5D0) cycle 1832 end if 1833 if (a(jatm)%index==1.and.istat==1) cycle !H+ doesn't contains electron and cannot compute density 1834 c80tmp="atmrad"//sep//trim(a(jatm)%name)//statname(istat)//".rad" 1835 inquire(file=c80tmp,exist=alive) 1836 if (alive.eqv. .false.) cycle 1837 open(10,file=c80tmp,status="old") 1838 read(10,*) atmradnpt(jatm) 1839 do ipt=1,atmradnpt(jatm) 1840 read(10,*) rnouse,atmradrho(jatm,ipt) 1841 end do 1842 close(10) 1843 do ipt=1+iradcut*sphpot,ntotpot 1844 atmstatgrid(iatm,ipt,jatm,istat)=fdens_rad(jatm,gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z) 1845 end do 1846 end do 1847 end do 1848 end do 1849end if 1850 1851!Set atomic initial radial density as neutral state, which is loaded from corresponding .rad file 1852atmradrho=0 1853do iatm=1,ncenter 1854 open(10,file="atmrad"//sep//trim(a(iatm)%name)//"_0.rad",status="old") 1855 read(10,*) atmradnpt(iatm) 1856 do ipt=1,atmradnpt(iatm) 1857 read(10,*) rnouse,atmradrho(iatm,ipt) 1858 end do 1859 close(10) 1860end do 1861 1862write(*,*) 1863write(*,*) "Performing Hirshfeld-I iteration to refine atomic spaces..." 1864lastcharge=0 1865!Cycle each atom to calculate their charges 1866do icyc=1,maxcyc 1867 if (ioutmedchg==1) write(*,*) 1868 if (icyc==1) then 1869 write(*,"(' Cycle',i5)") icyc 1870 else 1871 write(*,"(' Cycle',i5,' Maximum change:',f10.6)") icyc,varmax 1872 end if 1873 1874 do iatm=1,ncenter 1875 gridatm%value=gridatmorg%value !Weight in this grid point 1876 gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule 1877 gridatm%y=gridatmorg%y+a(iatm)%y 1878 gridatm%z=gridatmorg%z+a(iatm)%z 1879 1880 !Molecular density 1881 molrho=molrhoall(iatm,:) 1882 1883 !Calculate promolecular and proatomic density 1884 promol=0D0 1885 do jatm=1,ncenter 1886 if (ignorefar==1) then 1887 atmdist=dsqrt( (a(iatm)%x-a(jatm)%x)**2+(a(iatm)%y-a(jatm)%y)**2+(a(iatm)%z-a(jatm)%z)**2 ) 1888 if (atmdist>(vdwr(iatm)+vdwr(jatm))*1.5D0) cycle 1889 end if 1890 if (imode==1) then 1891nthreads=getNThreads() 1892!$OMP parallel do shared(tmpdens) private(ipt) num_threads(nthreads) 1893 do ipt=1+iradcut*sphpot,ntotpot 1894 tmpdens(ipt)=fdens_rad(jatm,gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z) 1895 end do 1896!$OMP end parallel do 1897 else if (imode==2) then 1898 if (icyc==1) then 1899 tmpdens=atmstatgrid(iatm,:,jatm,0) 1900 else 1901 ichglow=floor(lastcharge(jatm)) 1902 ichghigh=ceiling(lastcharge(jatm)) 1903 tmpdens=(lastcharge(jatm)-ichglow)*atmstatgrid(iatm,:,jatm,ichghigh) + (ichghigh-lastcharge(jatm))*atmstatgrid(iatm,:,jatm,ichglow) 1904 end if 1905 end if 1906 promol=promol+tmpdens 1907 if (jatm==iatm) selfdens=tmpdens 1908 end do 1909 1910 !Calculate atomic charge 1911 electmp=0D0 1912 do ipt=1+iradcut*sphpot,ntotpot 1913 if (promol(ipt)/=0D0) electmp=electmp+selfdens(ipt)/promol(ipt)*molrho(ipt)*gridatm(ipt)%value 1914 end do 1915 if (nEDFelec==0) charge(iatm)=a(iatm)%charge-electmp 1916 if (nEDFelec>0) charge(iatm)=a(iatm)%index-electmp !Assume EDF information provides inner-core electrons for all atoms using ECP 1917 if (ioutmedchg==1) write(*,"(' Charge of atom',i5,'(',a2,')',': ',f12.6,' Delta:',f12.6)") & 1918 iatm,a(iatm)%name,charge(iatm),charge(iatm)-lastcharge(iatm) 1919 end do 1920 1921 !Check convergence 1922 varmax=maxval(abs(charge-lastcharge)) 1923 if (varmax<crit) then 1924 if (itype==1) then 1925 write(*,"(a,f10.6)") " All atomic charges have converged to criterion of",crit 1926 write(*,"(' Sum of all charges:',f14.8)") sum(charge) 1927 !Normalize to get rid of integration inaccuracy 1928 totnumelec=sum(a%charge-charge) 1929 facnorm=nelec/totnumelec 1930 do iatm=1,ncenter 1931 charge(iatm)=a(iatm)%charge-facnorm*(a(iatm)%charge-charge(iatm)) 1932 end do 1933 write(*,*) 1934 write(*,*) "Final atomic charges, after normalization to actual number of electrons" 1935 do iatm=1,ncenter 1936 write(*,"(' Atom',i5,'(',a2,')',': ',f12.6)") iatm,a(iatm)%name,charge(iatm) 1937 end do 1938 exit 1939 else 1940 write(*,*) "Hirshfeld-I atomic spaces converged successfully!" 1941 write(*,*) 1942 return 1943 end if 1944 else 1945 if (icyc==maxcyc) then 1946 write(*,"(/,' Convergence failed within',i4,' cycles!')") maxcyc 1947 exit 1948 end if 1949 end if 1950 1951 !Update atomic radial density by means of interpolation of adjacent charge state 1952 do iatm=1,ncenter 1953 !Read radial density of lower limit state 1954 ichglow=floor(charge(iatm)) 1955 radrholow=0 1956 c80tmp="atmrad"//sep//trim(a(iatm)%name)//statname(ichglow)//".rad" 1957 inquire(file=c80tmp,exist=alive) 1958 if (alive.eqv. .false.) then 1959 write(*,"(' Error: ',a,' was not prepared!')") trim(c80tmp) 1960 return 1961 end if 1962 open(10,file=c80tmp,status="old") 1963 read(10,*) nptlow 1964 do ipt=1,nptlow 1965 read(10,*) rnouse,radrholow(ipt) 1966 end do 1967 close(10) 1968 !Read radial density of upper limit state 1969 ichghigh=ceiling(charge(iatm)) 1970 radrhohigh=0 1971 c80tmp="atmrad"//sep//trim(a(iatm)%name)//statname(ichghigh)//".rad" 1972 inquire(file=c80tmp,exist=alive) 1973 if (alive.eqv. .false.) then 1974 write(*,"(' Error: ',a,' was not prepared!')") trim(c80tmp) 1975 return 1976 end if 1977 open(10,file=c80tmp,status="old") 1978 read(10,*) npthigh 1979 do ipt=1,npthigh 1980 read(10,*) rnouse,radrhohigh(ipt) 1981 end do 1982 close(10) 1983 !Update current radial density 1984 atmradrho(iatm,:)=(charge(iatm)-ichglow)*radrhohigh(:) + (ichghigh-charge(iatm))*radrholow(:) 1985 atmradnpt(iatm)=max(npthigh,nptlow) 1986 end do 1987 1988 lastcharge=charge 1989end do 1990 1991if (allocated(frag1)) write(*,"(/,' Fragment charge:',f12.6)") sum(charge(frag1)) 1992CALL CPU_TIME(time_end) 1993call walltime(iwalltime2) 1994write(*,"(' Calculation took up CPU time',f12.2,'s, wall clock time',i10,'s')") time_end-time_begin,iwalltime2-iwalltime1 1995 1996call path2filename(firstfilename,chgfilename) 1997write(*,*) 1998write(*,"(a)") " If output atoms with charges to "//trim(chgfilename)//".chg in current folder? (y/n)" 1999read(*,*) selectyn 2000if (selectyn=="y".or.selectyn=="Y") then 2001 open(10,file=trim(chgfilename)//".chg",status="replace") 2002 do i=1,ncenter 2003 write(10,"(a4,4f12.6)") a(i)%name,a(i)%x*b2a,a(i)%y*b2a,a(i)%z*b2a,charge(i) 2004 end do 2005 close(10) 2006 write(*,"(a)") " Result have been saved to "//trim(chgfilename)//".chg in current folder" 2007 write(*,"(a)") " Columns ranging from 1 to 5 are name,X,Y,Z,charge respectively, unit is Angstrom" 2008end if 2009end subroutine 2010 2011 2012!!------- Generate atomic radial density files at different states, used for such as Hirshfeld-I 2013!"atmrad" in current folder is used as working directory 2014!-2,-1,0,+1,+2 charge states of each element will be calculated to produce atomic .wfn file by Gaussian, predefined ground state multiplicity is used 2015!After that, radial density file (.rad) is generated for each state of each element 2016!If atomic wfn file is already existed, calculation will be skipped 2017!Radial distance values are the same as built-in atomic density, i.e. those in atmraddens.f90 2018subroutine genatmradfile 2019use defvar 2020use util 2021implicit real*8 (a-h,o-z) 2022character c80tmp*80,c200tmp*200,calclevel*80,radname*200,sep 2023character*2 :: statname(-3:3)=(/ "-3","-2","-1","_0","+1","+2","+3" /) 2024integer :: chgmulti(nelesupp,-3:3)=0 !Ground state multiplicity of each charge state of each element. If value=0, means undefined 2025 2026!Define chgmulti for elements for possible states 2027!H,Li,Na,K,Rb,Cs 2028chgmulti(1,0)=2 2029chgmulti(1,1)=1 2030chgmulti(1,-1)=1 2031chgmulti(3,:)=chgmulti(1,:) 2032chgmulti(11,:)=chgmulti(1,:) 2033chgmulti(19,:)=chgmulti(1,:) 2034chgmulti(37,:)=chgmulti(1,:) 2035chgmulti(55,:)=chgmulti(1,:) 2036!He,Ne,Ar,Kr,Xe,Rn 2037chgmulti(2,0)=1 2038chgmulti(2,1)=2 2039chgmulti(2,-1)=2 2040chgmulti(10,:)=chgmulti(2,:) 2041chgmulti(18,:)=chgmulti(2,:) 2042chgmulti(36,:)=chgmulti(2,:) 2043chgmulti(54,:)=chgmulti(2,:) 2044chgmulti(86,:)=chgmulti(2,:) 2045!Be,Mg,Ca,Sr,Ba 2046chgmulti(4,0)=1 2047chgmulti(4,1)=2 2048chgmulti(4,2)=1 2049chgmulti(4,-1)=2 2050chgmulti(12,:)=chgmulti(4,:) 2051chgmulti(20,:)=chgmulti(4,:) 2052chgmulti(38,:)=chgmulti(4,:) 2053chgmulti(56,:)=chgmulti(4,:) 2054!B,Al,Ga,In,Tl 2055chgmulti(5,0)=2 2056chgmulti(5,1)=1 2057chgmulti(5,2)=2 2058chgmulti(5,-1)=3 2059chgmulti(5,-2)=4 2060chgmulti(13,:)=chgmulti(5,:) 2061chgmulti(31,:)=chgmulti(5,:) 2062chgmulti(49,:)=chgmulti(5,:) 2063chgmulti(81,:)=chgmulti(5,:) 2064!C,Si,Ge,Sn,Pb 2065chgmulti(6,0)=3 2066chgmulti(6,1)=2 2067chgmulti(6,2)=1 2068chgmulti(6,-1)=4 2069chgmulti(6,-2)=3 2070chgmulti(14,:)=chgmulti(6,:) 2071chgmulti(32,:)=chgmulti(6,:) 2072chgmulti(50,:)=chgmulti(6,:) 2073chgmulti(82,:)=chgmulti(6,:) 2074!N,P,As,Sb,Bi 2075chgmulti(7,0)=4 2076chgmulti(7,1)=3 2077chgmulti(7,2)=2 2078chgmulti(7,-1)=3 2079chgmulti(7,-2)=2 2080chgmulti(15,:)=chgmulti(7,:) 2081chgmulti(33,:)=chgmulti(7,:) 2082chgmulti(51,:)=chgmulti(7,:) 2083chgmulti(83,:)=chgmulti(7,:) 2084!O,S,Se,Te,Po 2085chgmulti(8,0)=3 2086chgmulti(8,1)=4 2087chgmulti(8,2)=3 2088chgmulti(8,-1)=2 2089chgmulti(8,-2)=1 2090chgmulti(16,:)=chgmulti(8,:) 2091chgmulti(34,:)=chgmulti(8,:) 2092chgmulti(52,:)=chgmulti(8,:) 2093chgmulti(84,:)=chgmulti(8,:) 2094!F,Cl,Br,I,At 2095chgmulti(9,0)=2 2096chgmulti(9,1)=3 2097chgmulti(9,2)=4 2098chgmulti(9,-1)=1 2099chgmulti(17,:)=chgmulti(9,:) 2100chgmulti(35,:)=chgmulti(9,:) 2101chgmulti(53,:)=chgmulti(9,:) 2102chgmulti(85,:)=chgmulti(9,:) 2103!Spin multiplicity of transition metal for each state is determined by chemical intuition as well as a few single point energy data 2104!For simplicity, I assume that later elements in each row has identical configuration, of course this is incorrect but not too bad 2105!Sc (3d1,4s2) 2106chgmulti(21,0)=2 2107chgmulti(21,1)=3 2108chgmulti(21,2)=2 2109chgmulti(21,-1)=3 2110chgmulti(39,:)=chgmulti(21,:) !Y 2111chgmulti(57,:)=chgmulti(21,:) !La 2112!Ti (3d2,4s2) 2113chgmulti(22,0)=3 2114chgmulti(22,1)=4 2115chgmulti(22,2)=3 2116chgmulti(22,-1)=4 2117chgmulti(40,:)=chgmulti(22,:) !Zr 2118chgmulti(72,:)=chgmulti(22,:) !Hf 2119!V (3d3,4s2) 2120chgmulti(23,0)=4 2121chgmulti(23,1)=5 2122chgmulti(23,2)=4 2123chgmulti(23,-1)=5 2124chgmulti(41,:)=chgmulti(23,:) !Nb 2125chgmulti(73,:)=chgmulti(23,:) !Ta 2126!Cr (3d5,4s1) 2127chgmulti(24,0)=7 2128chgmulti(24,1)=6 2129chgmulti(24,2)=5 2130chgmulti(24,-1)=6 2131chgmulti(42,:)=chgmulti(24,:) !Mo 2132chgmulti(74,:)=chgmulti(24,:) !W 2133!Mn (3d5,4s2) 2134chgmulti(25,0)=6 2135chgmulti(25,1)=7 2136chgmulti(25,2)=6 2137chgmulti(25,-1)=5 2138chgmulti(43,:)=chgmulti(25,:) !Tc 2139chgmulti(75,:)=chgmulti(25,:) !Re 2140!Fe (3d6,4s2) 2141chgmulti(26,0)=5 2142chgmulti(26,1)=6 2143chgmulti(26,2)=7 2144chgmulti(26,-1)=4 2145chgmulti(44,:)=chgmulti(26,:) !Ru 2146chgmulti(76,:)=chgmulti(26,:) !Os 2147!Co (3d7,4s2) 2148chgmulti(27,0)=4 2149chgmulti(27,1)=5 2150chgmulti(27,2)=4 2151chgmulti(27,-1)=3 2152chgmulti(45,:)=chgmulti(27,:) !Rh 2153chgmulti(77,:)=chgmulti(27,:) !Ir 2154!Ni (3d8,4s2) 2155chgmulti(28,0)=3 2156chgmulti(28,1)=4 2157chgmulti(28,2)=3 2158chgmulti(28,-1)=2 2159chgmulti(46,:)=chgmulti(28,:) !Pd 2160chgmulti(78,:)=chgmulti(28,:) !Pt 2161!Cu (3d10,4s1) 2162chgmulti(29,0)=2 2163chgmulti(29,1)=1 2164chgmulti(29,2)=2 2165chgmulti(29,-1)=1 2166chgmulti(47,:)=chgmulti(29,:) !Ag 2167chgmulti(79,:)=chgmulti(29,:) !Au 2168!Zn (3d10,4s2) 2169chgmulti(30,0)=1 2170chgmulti(30,1)=2 2171chgmulti(30,2)=1 2172chgmulti(30,-1)=2 2173chgmulti(48,:)=chgmulti(30,:) !Cd 2174chgmulti(80,:)=chgmulti(30,:) !Hg 2175 2176sep='/' !Separation symbol of directory 2177if (isys==1) sep='\' 2178calclevel=" " 2179! calclevel="B3LYP/def2SVP" 2180 2181!Cycle each charge state of each atom. Each element is only calculated once. If the file is existing, don't recalculate again 2182do iatm=1,ncenter 2183 iele=a(iatm)%index 2184 do istat=-3,3 2185 if (chgmulti(iele,istat)==0) cycle !Undefined state 2186 radname="atmrad"//sep//trim(a(iatm)%name)//statname(istat)//".wfn" 2187 inquire(file=radname,exist=alive) 2188 if (alive) cycle 2189 2190 !Check Gaussian path 2191 inquire(file=gaupath,exist=alive) 2192 if (.not.alive) then 2193 write(*,*) "Couldn't find Gaussian path defined in ""gaupath"" variable in settings.ini" 2194 if (isys==1) write(*,*) "Input the path of Gaussian executable file, e.g. ""d:\study\g09w\g09.exe""" 2195 if (isys==2.or.isys==3) write(*,*) "Input the path of Gaussian executable file, e.g. ""/sob/g09/g09""" 2196 do while(.true.) 2197 read(*,"(a)") gaupath 2198 inquire(file=gaupath,exist=alive) 2199 if (alive) exit 2200 write(*,*) "Couldn't find Gaussian executable file, input again" 2201 end do 2202 end if 2203 2204 !Input calculation level 2205 if (calclevel==" ") then 2206 write(*,*) "Some atomic .wfn files are not found in ""atmrad"" folder in current directory" 2207 write(*,"(a)") " Now input the level for calculating these .wfn files, e.g. B3LYP/def2SVP" 2208 write(*,"(a)") " You can also add other keywords at the same time, e.g. M062X/6-311G(2df,2p) scf=xqc int=ultrafine" 2209 read(*,"(a)") calclevel 2210 end if 2211 2212 !Generate .gjf file 2213 inquire(file="./atmrad/.",exist=alive) 2214 if (alive.eqv. .false.) call system("mkdir atmrad") 2215 c200tmp="atmrad"//sep//trim(a(iatm)%name)//statname(istat)//".gjf" 2216 open(10,file=c200tmp,status="replace") 2217 write(10,"(a)") "# "//trim(calclevel)//" out=wfn" 2218 write(10,*) 2219 write(10,"(a)") trim(a(iatm)%name)//statname(istat) 2220 write(10,*) 2221 write(10,"(2i3)") istat,chgmulti(iele,istat) 2222 write(10,"(a)") a(iatm)%name 2223 write(10,*) 2224 c200tmp="atmrad"//sep//trim(a(iatm)%name)//statname(istat)//".wfn" 2225 write(10,"(a)") trim(c200tmp) 2226 write(10,*) 2227 write(10,*) 2228 close(10) 2229 2230 !Start calculation 2231 c80tmp="atmrad"//sep//trim(a(iatm)%name)//statname(istat) 2232 write(*,*) "Running: "//trim(Gaupath)//' "'//trim(c80tmp)//'.gjf" "'//trim(c80tmp)//'"' 2233 call system(trim(Gaupath)//' "'//trim(c80tmp)//'.gjf" "'//trim(c80tmp)//'"') 2234 2235 !Check if Gaussian task was successfully finished 2236 if (isys==1) then 2237 inquire(file=trim(c80tmp)//".out",exist=alive) 2238 else 2239 inquire(file=trim(c80tmp)//".log",exist=alive) 2240 end if 2241 if (alive) then 2242 if (isys==1) then 2243 open(10,file=trim(c80tmp)//".out",status="old") 2244 else 2245 open(10,file=trim(c80tmp)//".log",status="old") 2246 end if 2247 call loclabel(10,"Normal termination",igaunormal) 2248 close(10) 2249 if (igaunormal==0) then 2250 write(*,"(a)") " Gaussian running may be failed! Please manually check Gaussian input and output files in atmrad folder" 2251 write(*,*) "Press ENTER to continue" 2252 read(*,*) 2253 end if 2254 else 2255 write(*,"(a)") " Gaussian running may be failed! Please manually check Gaussian input and output files in atmrad folder" 2256 write(*,*) "Press ENTER to continue" 2257 read(*,*) 2258 end if 2259 end do 2260end do 2261 2262!All element wfn files have been generated, now calculate corresponding radial density file (.rad) 2263!Existing .rad file will not be re-calculated 2264write(*,*) 2265write(*,*) "Generating atomic radial density from atomic wfn file..." 2266do iatm=1,ncenter 2267 iele=a_org(iatm)%index 2268 do istat=-3,3 2269 if (chgmulti(iele,istat)==0) cycle !Undefined state 2270 c80tmp="atmrad"//sep//trim(a_org(iatm)%name)//statname(istat) 2271 inquire(file=trim(c80tmp)//".rad",exist=alive) 2272 if (alive) cycle 2273 inquire(file=trim(c80tmp)//".wfn",exist=alive) 2274 if (alive.eqv. .false.) then 2275 write(*,"(' Error: ',a,' was not found!')") trim(c80tmp)//".wfn" 2276 write(*,*) "If you want to skip, press ENTER directly" 2277 read(*,*) 2278 cycle 2279 end if 2280 write(*,"(' Converting ',a,' to ',a)") trim(c80tmp)//".wfn",trim(c80tmp)//".rad" 2281 call atmwfn2atmrad(trim(c80tmp)//".wfn",trim(c80tmp)//".rad") 2282 end do 2283end do 2284 2285!Recover to the firstly loaded file 2286call dealloall 2287call readinfile(firstfilename,1) 2288end subroutine 2289 2290 2291!!----- Generate atomic radial density from atomic wfn file 2292!The code is adapted from sphatmraddens 2293subroutine atmwfn2atmrad(infile,outfile) 2294use defvar 2295use function 2296implicit real*8 (a-h,o-z) 2297character(len=*) infile,outfile 2298real*8,allocatable :: potx(:),poty(:),potz(:),potw(:),radpos(:),sphavgval(:) 2299call dealloall 2300call readinfile(infile,1) 2301truncrho=1D-8 2302rlow=0D0 2303rhigh=12 2304nsphpt=974 2305nradpt=200 !Totally 200 radial points, but the number of point is truncated at truncrho (because the interpolation routine doesn't work well for very low value) 2306allocate(potx(nsphpt),poty(nsphpt),potz(nsphpt),potw(nsphpt),radpos(nradpt),sphavgval(nradpt)) 2307sphavgval=0 2308call Lebedevgen(nsphpt,potx,poty,potz,potw) 2309nthreads=getNThreads() 2310!$OMP PARALLEL DO SHARED(sphavgval,radpos) PRIVATE(irad,radx,radr,isph,rnowx,rnowy,rnowz) schedule(dynamic) NUM_THREADS(nthreads) 2311do irad=1,nradpt 2312 radx=cos(irad*pi/(nradpt+1)) 2313 radr=(1+radx)/(1-radx) !Becke transform 2314 radpos(irad)=radr 2315 do isph=1,nsphpt 2316 rnowx=potx(isph)*radr 2317 rnowy=poty(isph)*radr 2318 rnowz=potz(isph)*radr 2319 sphavgval(irad)=sphavgval(irad)+fdens(rnowx,rnowy,rnowz)*potw(isph) 2320 end do 2321end do 2322!$OMP END PARALLEL DO 2323open(10,file=outfile,status="replace") 2324write(10,*) count(sphavgval>truncrho) 2325do irad=nradpt,1,-1 2326 if (sphavgval(irad)>truncrho) write(10,"(f20.12,E18.10)") radpos(irad),sphavgval(irad) 2327end do 2328close(10) 2329end subroutine 2330 2331 2332!!---- Calculate density at a point for iatm based on loaded atomic radial density 2333real*8 function fdens_rad(iatm,x,y,z) 2334use defvar 2335use util 2336integer iatm,npt 2337real*8 x,y,z,r,rnouse 2338npt=atmradnpt(iatm) 2339r=dsqrt((a(iatm)%x-x)**2+(a(iatm)%y-y)**2+(a(iatm)%z-z)**2) 2340call lagintpol(atmradpos(1:npt),atmradrho(iatm,1:npt),npt,r,fdens_rad,rnouse,rnouse,1) 2341end function 2342 2343 2344 2345 2346!!---------------------------------------- 2347!!--------- Calculate EEM charge --------- 2348!!---------------------------------------- 2349subroutine EEM 2350use defvar 2351use util 2352integer,parameter :: maxBO=3 !Maximum of possible bond order 2353character c200tmp*200 2354real*8 EEMmat(ncenter+1,ncenter+1),EEMarr(ncenter+1),qarr(ncenter+1) 2355real*8 kappa,Aparm(nelesupp,maxBO),Bparm(nelesupp,maxBO) !If parameter is -1, means undefined parameter 2356real*8 :: chgnet=0 2357 2358if (ifiletype/=11) then 2359 write(*,"(a)") " Error: MDL Molfile (.mol) file must be used as input file, since it contains atomic connectivity information!" 2360 write(*,*) "Press Enter to return" 2361 read(*,*) 2362 return 2363end if 2364 2365iparmset=2 2366call genEEMparm(iparmset,kappa,Aparm,Bparm) 2367isel2=-10 2368 2369EEMcyc: do while(.true.) 2370 write(*,*) 2371 write(*,*) " ====== Electronegativity Equalization Method (EEM) ======" 2372 write(*,*) "-1 Return" 2373 write(*,*) "0 Start calculation" 2374 write(*,*) "1 Choose EEM parameters" 2375 write(*,"(a,f4.1)") " 2 Set net charge, current:",chgnet 2376 read(*,*) isel 2377 if (isel==-1) then 2378 return 2379 else if (isel==1) then 2380 do while(.true.) 2381 if (isel2/=-1) then 2382 write(*,*) 2383 write(*,*) "Present EEM parameters:" 2384 write(*,"(' kappa',f12.6)") kappa 2385 do iele=1,nelesupp 2386 do imulti=1,maxBO 2387 if (Aparm(iele,imulti)/=-1) write(*,"(1x,a,' Multiplicity:',i2,' A:',f9.5, ' B:',f9.5)") ind2name(iele),imulti,Aparm(iele,imulti),Bparm(iele,imulti) 2388 end do 2389 end do 2390 end if 2391 write(*,*) 2392 write(*,*) "-2 Return" 2393 write(*,*) "-1 Export present parameters to external file" 2394 write(*,*) "0 Load parameters from external file" 2395 write(*,*) "1 Use parameters fitted to HF/STO-3G Mulliken charge, IJMS, 8, 572 (2007)" 2396 write(*,*) "2 Use parameters fitted to B3LYP/6-31G* CHELPG charge, JCC, 30, 1174 (2009)" 2397 write(*,*) "3 Use parameters fitted to HF/6-31G* CHELPG charge, JCC, 30, 1174 (2009)" 2398 write(*,*) "4 Use parameters fitted to B3LYP/6-311G* NPA charge, J Cheminform, 8, 57(2016)" 2399 read(*,*) isel2 2400 if (isel2==-2) then 2401 exit 2402 else if (isel2==-1) then 2403 open(10,file="EEMparm.txt",status="replace") 2404 write(10,"(f12.6)") kappa 2405 do iele=1,nelesupp 2406 do imulti=1,maxBO 2407 if (Aparm(iele,imulti)/=-1) write(10,"(1x,a,i3,2f9.5)") ind2name(iele),imulti,Aparm(iele,imulti),Bparm(iele,imulti) 2408 end do 2409 end do 2410 close(10) 2411 write(*,*) "Parameters have been exported to EEMparm.txt in current folder" 2412 else if (isel2==0) then 2413 write(*,*) "Input path of parameter file, e.g. C:\aqours.txt" 2414 do while(.true.) 2415 read(*,"(a)") c200tmp 2416 inquire(file=c200tmp,exist=alive) 2417 if (alive) exit 2418 write(*,*) "Cannot find the file, input again" 2419 end do 2420 Aparm=-1 2421 Bparm=-1 2422 open(10,file=c200tmp,status="old") 2423 read(10,*) kappa 2424 nload=0 2425 do while(.true.) 2426 read(10,*,iostat=ierror) c200tmp,imulti,tmpA,tmpB 2427 if (ierror/=0) exit 2428 call lc2uc(c200tmp(1:1)) !Convert to upper case 2429 call uc2lc(c200tmp(2:2)) !Convert to lower case 2430 do iele=1,nelesupp 2431 if (c200tmp(1:2)==ind2name(iele)) exit 2432 end do 2433 Aparm(iele,imulti)=tmpA 2434 Bparm(iele,imulti)=tmpB 2435 nload=nload+1 2436 end do 2437 write(*,"(' Loaded',i5,' entries')") nload 2438 close(10) 2439 else 2440 call genEEMparm(isel2,kappa,Aparm,Bparm) 2441 end if 2442 end do 2443 else if (isel==2) then 2444 write(*,*) "Input net charge of the system, e.g. -1" 2445 read(*,*) chgnet 2446 else if (isel==0) then 2447 write(*,*) "Calculating..." 2448 write(*,*) 2449 !Construct EEM array 2450 EEMarr(ncenter+1)=chgnet 2451 do iatm=1,ncenter 2452 imulti=maxval(connmat(iatm,:)) 2453 if (imulti>maxBO) then 2454 write(*,"(' Error: Multiplicity of atom',i5,' exceeded upper limit (',i2,')')") iatm,maxBO 2455 cycle EEMcyc 2456 end if 2457 tmpval=Aparm(a(iatm)%index,imulti) 2458 if (tmpval==-1) then 2459 write(*,"(' Error: Parameter for atom',i5,'(',a,') is unavailable!')") iatm,a(iatm)%name 2460 cycle EEMcyc 2461 else 2462 EEMarr(iatm)=-tmpval 2463 end if 2464 end do 2465 2466 !Construct EEM matrix 2467 EEMmat=0 2468 EEMmat(ncenter+1,1:ncenter)=1 2469 EEMmat(1:ncenter,ncenter+1)=-1 2470 do i=1,ncenter 2471 imulti=maxval(connmat(i,:)) 2472 do j=1,ncenter 2473 if (i==j) then 2474 EEMmat(i,j)=Bparm(a(i)%index,imulti) 2475 else 2476 EEMmat(i,j)=kappa/(distmat(i,j)*b2a) 2477 end if 2478 end do 2479 end do 2480 2481 !Solve EEM equation 2482 qarr=matmul(invmat(EEMmat,ncenter+1),EEMarr) 2483 do iatm=1,ncenter 2484 write(*,"(' EEM charge of atom',i5,'(',a,'):',f12.6)") iatm,a(iatm)%name,qarr(iatm) 2485 end do 2486 write(*,"(' Electronegativity:',f12.6)") qarr(ncenter+1) 2487 end if 2488 2489end do EEMcyc 2490end subroutine 2491 2492!---- Generate EEM parameters 2493subroutine genEEMparm(iset,kappa,Aparm,Bparm) 2494use defvar 2495integer,parameter :: maxBO=3 !Maximum of bond order 2496real*8 kappa,Aparm(nelesupp,maxBO),Bparm(nelesupp,maxBO) 2497Aparm=-1 2498Bparm=-1 2499if (iset==1) then !Parameters fitted to Mulliken charge at HF/STO-3G, Int. J. Mol. Sci., 8, 572-582 (2007) 2500 write(*,"(a)") " Parameters have been set to those fitted to Mulliken charges at HF/STO-3G, see Int. J. Mol. Sci., 8, 572-582 (2007)" 2501 kappa=0.44D0 2502 Aparm(1,1)= 2.396D0 !H 2503 Bparm(1,1)= 0.959D0 2504 Aparm(6,1)= 2.459D0 !C,multi=1 2505 Bparm(6,1)= 0.611D0 2506 Aparm(7,1)= 2.597D0 !N,multi=1 2507 Bparm(7,1)= 0.790D0 2508 Aparm(8,1)= 2.625D0 !O,multi=1 2509 Bparm(8,1)= 0.858D0 2510 Aparm(16,1)= 2.407D0 !S,multi=1 2511 Bparm(16,1)= 0.491D0 2512 Aparm(6,2)= 2.464D0 !C,multi=2 2513 Bparm(6,2)= 0.565D0 2514 Aparm(7,2)= 2.554D0 !N,multi=2 2515 Bparm(7,2)= 0.611D0 2516 Aparm(8,2)= 2.580D0 !O,multi=2 2517 Bparm(8,2)= 0.691D0 2518else if (iset==2) then !Parameters fitted to CHELPG charges at B3LYP/6-31G*, J. Comput. Chem., 30, 1174 (2009) 2519 write(*,"(a)") " Parameters have been set to those fitted to CHELPG charges at B3LYP/6-31G*, see J. Comput. Chem., 30, 1174 (2009)" 2520 kappa=0.302D0 2521 Aparm(35,1)= 2.659D0 !Br,multi=1 2522 Bparm(35,1)= 1.802D0 2523 Aparm(6,1)= 2.482D0 !C,multi=1 2524 Bparm(6,1)= 0.464D0 2525 Aparm(17,1)= 2.519D0 !Cl,multi=1 2526 Bparm(17,1)= 1.450D0 2527 Aparm(9,1)= 3.577D0 !F,multi=1 2528 Bparm(9,1)= 3.419D0 2529 Aparm(1,1)= 2.385D0 !H,multi=1 2530 Bparm(1,1)= 0.737D0 2531 Aparm(7,1)= 2.595D0 !N,multi=1 2532 Bparm(7,1)= 0.468D0 2533 Aparm(8,1)= 2.825D0 !O,multi=1 2534 Bparm(8,1)= 0.844D0 2535 Aparm(16,1)= 2.452D0 !S,multi=1 2536 Bparm(16,1)= 0.362D0 2537 Aparm(30,1)= 2.298D0 !Zn,multi=1 2538 Bparm(30,1)= 0.420D0 2539 Aparm(6,2)= 2.464D0 !C,multi=2 2540 Bparm(6,2)= 0.392D0 2541 Aparm(7,2)= 2.556D0 !N,multi=2 2542 Bparm(7,2)= 0.377D0 2543 Aparm(8,2)= 2.789D0 !O,multi=2 2544 Bparm(8,2)= 0.834D0 2545else if (iset==3) then !Parameters fitted to CHELPG charges at HF/6-31G*, J. Comput. Chem., 30, 1174 (2009) 2546 write(*,"(a)") " Parameters have been set to those fitted to CHELPG charges at HF/6-31G*, see J. Comput. Chem., 30, 1174 (2009)" 2547 kappa=0.227D0 2548 Aparm(35,1)= 2.615D0 !Br,multi=1 2549 Bparm(35,1)= 1.436D0 2550 Aparm(6,1)= 2.481D0 !C,multi=1 2551 Bparm(6,1)= 0.373D0 2552 Aparm(17,1)= 2.517D0 !Cl,multi=1 2553 Bparm(17,1)= 1.043D0 2554 Aparm(9,1)= 3.991D0 !F,multi=1 2555 Bparm(9,1)= 3.594D0 2556 Aparm(1,1)= 2.357D0 !H,multi=1 2557 Bparm(1,1)= 0.688D0 2558 Aparm(7,1)= 2.585D0 !N,multi=1 2559 Bparm(7,1)= 0.329D0 2560 Aparm(8,1)= 2.870D0 !O,multi=1 2561 Bparm(8,1)= 0.717D0 2562 Aparm(16,1)= 2.450D0 !S,multi=1 2563 Bparm(16,1)= 0.269D0 2564 Aparm(30,1)= 2.185D0 !Zn,multi=1 2565 Bparm(30,1)= 0.375D0 2566 Aparm(6,2)= 2.475D0 !C,multi=2 2567 Bparm(6,2)= 0.292D0 2568 Aparm(7,2)= 2.556D0 !N,multi=2 2569 Bparm(7,2)= 0.288D0 2570 Aparm(8,2)= 2.757D0 !O,multi=2 2571 Bparm(8,2)= 0.621D0 2572else if (iset==4) then !Parameters fitted to NPA charges at B3LYP/6-311G*, see J. Cheminform., 8, 57 (2016) 2573!The data were taken from "13321_2016_171_MOESM5_ESM.rar Additional file 5" of supplmental material 2574!13321_2016_171_MOESM5_ESM\neemp\ideal_q1\set3_DE_RMSD_B3LYP_6311Gd_NPA_cross_ideal\output_set3_DE_RMSD_B3LYP_6311Gd_NPA_cross_ideal_all 2575 write(*,"(a)") " Parameters have been set to those fitted to NPA charges at B3LYP/6-311G*, they were extracted from SI of J. Cheminform., 8, 57 (2016)" 2576 kappa=0.4024D0 2577 Aparm(1,1)= 2.4598D0 !H,multi=1 2578 Bparm(1,1)= 0.9120D0 2579 Aparm(6,1)= 2.5957D0 !C,multi=1 2580 Bparm(6,1)= 0.5083D0 2581 Aparm(6,2)= 2.6029D0 !C,multi=2 2582 Bparm(6,2)= 0.5021D0 2583 Aparm(6,3)= 2.5326D0 !C,multi=3 2584 Bparm(6,3)= 0.5932D0 2585 Aparm(7,1)= 2.7802D0 !N,multi=1 2586 Bparm(7,1)= 0.7060D0 2587 Aparm(7,2)= 2.7141D0 !N,multi=2 2588 Bparm(7,2)= 0.5366D0 2589 Aparm(7,3)= 2.6391D0 !N,multi=3 2590 Bparm(7,3)= 0.5171D0 2591 Aparm(8,1)= 2.9496D0 !O,multi=1 2592 Bparm(8,1)= 0.8264D0 2593 Aparm(8,2)= 2.8595D0 !O,multi=2 2594 Bparm(8,2)= 0.6589D0 2595 Aparm(9,1)= 2.9165D0 !F,multi=1 2596 Bparm(9,1)= 0.8427D0 2597 Aparm(15,2)= 2.1712D0 !P,multi=2 2598 Bparm(15,2)= 0.4802D0 2599 Aparm(16,1)= 2.5234D0 !S,multi=1 2600 Bparm(16,1)= 0.3726D0 2601 Aparm(16,2)= 2.5334D0 !S,multi=2 2602 Bparm(16,2)= 0.3519D0 2603 Aparm(17,1)= 2.5625D0 !Cl,multi=1 2604 Bparm(17,1)= 0.9863D0 2605 Aparm(35,1)= 2.4772D0 !Br,multi=1 2606 Bparm(35,1)= 1.2131D0 2607end if 2608end subroutine 2609