1!------------- Modify & Check wavefunction 2subroutine modwfn 3use defvar 4use util 5implicit real*8 (a-h,o-z) 6character seltmpc*10,orbtype*5,selectyn,c1000tmp*1000,c2000tmp*2000 7real*8 eigval(nbasis),eigvec(nbasis,nbasis),tmpmat(nbasis,nbasis) 8integer orbarr(nmo) 9integer,allocatable :: exclfragatm(:),tmparrint(:) 10 11do while(.true.) 12 write(*,*) " ============ Modify & Check wavefunction ============ " 13 write(*,"(' Number of GTFs:',i6,', Orb:',i6,', Atoms:',i5,', A/B elec:',2f8.3)") nprims,nmo,ncenter,naelec,nbelec 14 if (ifragcontri/=1) write(*,*) "-4 Exclude contribution of some atoms to real space functions" 15 if (ifragcontri/=1) write(*,*) "-3 Only retain contribution of some atoms to real space functions" 16 write(*,*) "-1 Return" 17 write(*,*) "0 Save the modified wavefunction to a new .wfn file" 18 if (allocated(CObasa)) then 19 write(*,*) "1 List all primitive function 2 List all basis function" 20 else 21 write(*,*) "1 List all primitive function" 22 end if 23 write(*,*) "3 List all orbitals 4 Print detail information of an orbital" 24 if (allocated(CObasa)) write(*,*) "5 Print coefficient matrix in basis function" 25 if (allocated(CObasa)) write(*,*) "6 Print density matrix in basis function" 26 if (allocated(CObasa)) write(*,*) "7 Print various kinds of integral matrix between basis functions" 27 write(*,*) "11 Swap some information of two primitive functions" 28 write(*,*) "21 Set center of a primitive 22 Set type of a primitive" 29 write(*,*) "23 Set exponent of a primitive 24 Set coefficient of a primitive" 30 if (allocated(CObasa)) then 31 write(*,*) "25 Set coefficients of GTFs/basis functions that satisfied certain conditions" 32 else 33 write(*,*) "25 Set coefficients of GTFs that satisfied certain conditions" 34 end if 35 write(*,*) "26 Set occupation number of some orbitals" 36 write(*,*) "27 Set type of some orbitals" 37 write(*,*) "31 Translate the system 32 Translate and duplicate the system" 38 write(*,*) "33 Rotate wavefunction, namely X->Y, Y->Z, Z->X" 39 if (imodwfn==0) write(*,*) "34 Set occupation number of inner orbitals to zero" !If occupation has been modified, don't do this to complicate things 40 if (allocated(MOsym)) write(*,*) "35 Keep or discard orbital contributions according to irreducible rep." 41 write(*,*) "36 Invert phase of some orbitals" 42 read(*,*) iselect 43 44 write(*,*) 45 if (iselect==-1) then 46 if (allocated(CObasa).and.imodwfn==1) then 47 write(*,*) "Updating density matrix..." 48 call genP 49 write(*,*) "Density matrix has been updated" 50 end if 51 exit 52 53 else if (iselect==-3.or.iselect==-4) then 54 deallocate(fragatm) !fragatm has been defined previously by default, fragatm contains all atoms 55 if (iselect==-3) then 56 ! "fragatm" is convertion relationship from fragment to the whole, 57 ! e.g. fragatm(4) is the actual atom index corresponding the 4th atom in fragment list 58 write(*,"(a)") " Input atomic indices to define the fragment, e.g. 1,3-6,8,10-11 means atoms 1,3,4,5,6,8,10,11 will constitute the fragment" 59 read(*,"(a)") c2000tmp 60 call str2arr(c2000tmp,nfragatmnum) 61 allocate(fragatm(nfragatmnum)) 62 call str2arr(c2000tmp,nfragatmnum,fragatm) 63 call sorti4(fragatm,"val") 64 else if (iselect==-4) then 65 write(*,*) "How many atoms will be excluded?" 66 write(*,*) "e.g. 1,3-6,8,10-11 means the atoms 1,3,4,5,6,8,10,11 will be excluded" 67 read(*,"(a)") c2000tmp 68 call str2arr(c2000tmp,nexclatm) 69 nfragatmnum=ncenter-nexclatm 70 allocate(fragatm(nfragatmnum),exclfragatm(nexclatm)) 71 call str2arr(c2000tmp,nexclatm,exclfragatm) 72 j=0 73 do i=1,ncenter 74 if (all(exclfragatm/=i)) then 75 j=j+1 76 fragatm(j)=i 77 end if 78 end do 79 end if 80 j=0 81 do i=1,nprims 82 if (any(fragatm==b(i)%center)) then 83 j=j+1 !Move function in the fragment to head of list 84 CO(:,j)=CO(:,i) 85 b(j)=b(i) 86 end if 87 end do 88 ifragcontri=1 !Fragment has been defined by users 89 write(*,"(' Done,',i8,' GTFs have been discarded,',i8,' GTFs reserved')") nprims-j,j 90 nprims=j !Cut list at j, all functions after j seem non exist 91 if (iselect==-4) deallocate(exclfragatm) 92 93 !Modification of wavefunction has finished, now reduce size of b, CO... to current nprims and nmo to avoid potential problems 94 if (allocated(b)) then !Only for input file contains wavefunctions 95 call resizebynmo(nmo,nprims) !Reduce size of CO, MOene, MOocc, MOtype 96 allocate(b_tmp(nprims)) 97 b_tmp(:)=b(1:nprims) 98 deallocate(b) 99 allocate(b(nprims)) 100 b=b_tmp 101 deallocate(b_tmp) 102 end if 103 else if (iselect==0) then 104 call outwfn("new.wfn",1,1,10) 105 write(*,*) "New .wfn file has been outputted to new.wfn in current folder" 106 107 else if (iselect==1) then 108 do i=1,nprims 109 write(*,"(i6,' Center:',i5,'(',a2,')',' Type: ',a,' Exponent:',D16.7)") i,b(i)%center,a(b(i)%center)%name,GTFtype2name(b(i)%functype),b(i)%exp 110 end do 111 112 else if (iselect==2) then 113 do i=1,nbasis 114 write(*,"(' Basis:',i5,' Shell:',i5,' Center:',i5,'(',a2,') Type:',a)")& 115 i,basshell(i),bascen(i),a(bascen(i))%name,GTFtype2name(bastype(i)) 116 end do 117 118 else if (iselect==3) then 119 do i=1,nmo 120 if (MOtype(i)==0) orbtype="A+B" 121 if (MOtype(i)==1) orbtype="A" 122 if (MOtype(i)==2) orbtype="B" 123 if (allocated(MOsym)) then 124 write(*,"(' Orb:',i6,' Ene(au/eV):',f13.6,f13.4,' Occ:',f9.6,' Type: ',a,' (',a,')')") i,MOene(i),MOene(i)*au2eV,MOocc(i),trim(orbtype),trim(MOsym(i)) 125 else 126 write(*,"(' Orb:',i6,' Ene(au/eV):',f13.6,f13.4,' Occ:',f9.6,' Type: ',a)") i,MOene(i),MOene(i)*au2eV,MOocc(i),trim(orbtype) 127 end if 128 end do 129 else if (iselect==4) then 130 write(*,*) "Input the orbital index, e.g. 12" 131 read(*,*) i 132 if (i<1.or.i>nmo) then 133 write(*,"(' Invalid orbital index, should within range of',i5,' and ',i5)") 1,nmo 134 else 135 write(*,"(' Occupation number is ',f12.7,' Energy is',f12.6,' Hartree')") MOocc(i),MOene(i) 136 if (MOtype(i)==0) write(*,*) "This is a close-shell orbital" 137 if (MOtype(i)==1) write(*,*) "This is an alpha orbital" 138 if (MOtype(i)==2) write(*,*) "This is a beta orbital" 139 write(*,*) 140 do j=1,nprims 141 write(*,"(' GTF:',i6,' Cen:',i5,'(',a2,')',' Type: ',a,' Coeff:',1PD15.8,' Exp: ',1PD13.7)") & 142 j,b(j)%center,a(b(j)%center)%name,GTFtype2name(b(j)%functype),co(i,j),b(j)%exp 143 end do 144 write(*,"(a,/)") " Note: The ""coeff."" are expansion coefficients of orbitals with respect to GTFs, including normalization constant" 145 if (allocated(b)) then 146 do j=1,nbasis 147 if (MOtype(i)==0.or.MOtype(i)==1) covalue=CObasa(j,i) 148 if (MOtype(i)==2) covalue=CObasb(j,i-nbasis) 149 write(*,"(' Basis func:',i6,' Cen:',i5,'(',a2,')',' Shell:',i5,' Type: ',a,' Coeff:',f12.8)") & 150 j,bascen(j),a(bascen(j))%name,basshell(j),GTFtype2name(bastype(j)),covalue 151 end do 152 end if 153 write(*,"(a,/)") " Note: The ""coeff."" are expansion coefficients of orbitals with respect to basis functions, which are normalized functions" 154 end if 155 156 else if (iselect==5) then 157 write(*,*) "0 Return" 158 write(*,*) "1 Print on screen" 159 write(*,*) "2 Print to Cmat.txt in current folder" 160 read(*,*) iseltmp 161 if (iseltmp==1.or.iseltmp==2) then 162 if (iseltmp==1) ides=6 163 if (iseltmp==2) then 164 ides=10 165 open(ides,file="Cmat.txt",status="replace") 166 end if 167 write(ides,*) "Note: (i,j) element means coefficient of ith basis function in jth orbital" 168 if (wfntype==0.or.wfntype==2.or.wfntype==3) then 169 call showmatgau(cobasa,"Coefficient matrix",0,fileid=ides) 170 else if (wfntype==1.or.wfntype==4) then 171 call showmatgau(cobasa,"Alpha coefficient matrix",0,fileid=ides) 172 call showmatgau(cobasb,"Beta coefficient matrix",0,fileid=ides) 173 end if 174 if (iseltmp==2) then 175 write(*,*) "Done! The matrix has been outputted to Cmat.txt in current folder" 176 close(ides) 177 end if 178 end if 179 180 else if (iselect==6) then 181 write(*,*) "0 Return" 182 write(*,*) "1 Print on screen" 183 write(*,*) "2 Print to Pmat.txt in current folder" 184 read(*,*) iseltmp 185 if (iseltmp==1.or.iseltmp==2) then 186 if (iseltmp==1) ides=6 187 if (iseltmp==2) then 188 ides=10 189 open(ides,file="Pmat.txt",status="replace") 190 end if 191 call showmatgau(Ptot,"Total density matrix",1,fileid=ides) 192 sumt=0 193 do i=1,nbasis 194 sumt=sumt+Ptot(i,i) 195 end do 196 write(ides,"(' The trace of the density matrix:',f12.6)") sumt 197 write(ides,"(' The trace of the density matrix multiplied by overlap matrix:',f12.6)") sum(Ptot*Sbas) 198 if (wfntype==1.or.wfntype==2.or.wfntype==4) then 199 suma=0 200 sumb=0 201 do i=1,nbasis 202 suma=suma+Palpha(i,i) 203 sumb=sumb+Pbeta(i,i) 204 end do 205 write(ides,*) 206 call showmatgau(Palpha-Pbeta,"Spin density matrix",1,fileid=ides) 207 write(ides,*) 208 call showmatgau(Palpha,"Alpha density matrix",1,fileid=ides) 209 write(ides,*) 210 call showmatgau(Pbeta,"Beta density matrix",1,fileid=ides) 211 write(ides,*) 212 write(ides,"(' The trace of the alpha and beta density matrix:',2f12.6)") suma,sumb 213 write(ides,"(' The trace of the density matrix multiplied by overlap matrix:',/,' Alpha:',f12.6,' Beta:',f12.6)") sum(Palpha*Sbas),sum(Pbeta*Sbas) 214 end if 215 if (iseltmp==2) then 216 write(*,*) "Done! The matrix has been outputted to Pmat.txt in current folder" 217 close(ides) 218 end if 219 end if 220 221 !Print various kinds of integral matrix between basis functions 222 else if (iselect==7) then 223 write(*,*) "Print which kind of integral matrix?" 224 write(*,*) "1 Overlap integral" 225 write(*,*) "2 Electric dipole moment integral" 226 write(*,*) "3 Magnetic dipole moment integral" 227 write(*,*) "4 Velocity integral" 228 write(*,*) "5 Kinetic energy integral" 229 read(*,*) imattype 230 write(*,*) "Select destination for output" 231 write(*,*) "1 Print on screen" 232 write(*,*) "2 Print to intmat.txt in current folder" 233 read(*,*) iout 234 if (iout==1) ides=6 235 if (iout==2) then 236 ides=10 237 open(ides,file="intmat.txt",status="replace") 238 end if 239 if (imattype==1) then 240 call showmatgau(Sbas,"Overlap matrix",1,fileid=ides) 241 call diagsymat(Sbas,eigvec,eigval,ierror) 242 write(ides,*) 243 write(ides,*) "Eigenvalues:" 244 write(ides,"(6f12.8)") eigval 245 else if (imattype==2) then 246 if (allocated(Dbas)) then 247 write(ides,*) 248 call showmatgau(Dbas(1,:,:),"Electric dipole moment matrix (X component)",1,fileid=ides) 249 write(ides,*) 250 call showmatgau(Dbas(2,:,:),"Electric dipole moment matrix (Y component)",1,fileid=ides) 251 write(ides,*) 252 call showmatgau(Dbas(3,:,:),"Electric dipole moment matrix (Z component)",1,fileid=ides) 253 else 254 write(*,"(a)") " Error: Electric dipole moment integral matrix has not been calculated, please set ""igenDbas"" & 255 in settings.ini to 1, so that this matrix can be generated when loading input file" 256 write(*,*) "Press Enter to skip" 257 read(*,*) 258 cycle 259 end if 260 else if (imattype==3) then 261 if (allocated(Magbas)) then 262 write(ides,*) 263 call showmatgau(Magbas(1,:,:),"Magnetic dipole moment matrix (X component)",0,fileid=ides) 264 write(ides,*) 265 call showmatgau(Magbas(2,:,:),"Magnetic dipole moment matrix (Y component)",0,fileid=ides) 266 write(ides,*) 267 call showmatgau(Magbas(3,:,:),"Magnetic dipole moment matrix (Z component)",0,fileid=ides) 268 else 269 write(*,"(a)") " Error: Magnetic dipole moment integral matrix has not been calculated, please set ""igenMagbas"" & 270 in settings.ini to 1, so that this matrix can be generated when loading input file" 271 write(*,*) "Press Enter to skip" 272 read(*,*) 273 cycle 274 end if 275 else if (imattype==4.or.imattype==5) then 276 if (isphergau==1) then !If you want, you can generate the matrix and perform Cartesian->spherical transformation at file loading stage for Velbas 277 write(*,"(a)") " Error: Spherical-harmonic type of basis functions are found. This function only works when all basis functions are Cartesian type!" 278 write(*,*) "Press Enter to skip" 279 read(*,*) 280 cycle 281 end if 282 if (imattype==4) then 283 if (.not.allocated(Velbas)) then 284 write(*,*) "Calculating velocity matrix..." 285 allocate(Velbas(3,nbasis,nbasis)) 286 call genvelbas 287 end if 288 write(ides,*) 289 call showmatgau(Velbas(1,:,:),"Velocity matrix (X component)",0,fileid=ides) 290 write(ides,*) 291 call showmatgau(Velbas(2,:,:),"Velocity matrix (Y component)",0,fileid=ides) 292 write(ides,*) 293 call showmatgau(Velbas(3,:,:),"Velocity matrix (Z component)",0,fileid=ides) 294 else if (imattype==5) then 295 if (.not.allocated(Tbas)) then 296 write(*,*) "Calculating Kinetic energy matrix..." 297 allocate(Tbas(nbasis,nbasis)) 298 call genTbas 299 end if 300 write(ides,*) 301 call showmatgau(Tbas(:,:),"Kinetic energy matrix",1,fileid=ides) 302 end if 303 end if 304 305 if (iout==2) then 306 write(*,*) "Done! The matrix has been outputted to intmat.txt in current folder" 307 close(ides) 308 end if 309 310 else if (iselect==11) then 311 write(*,*) "Swap information of which two GTFs? Input their indices e.g. 18,21" 312 read(*,*) i,j 313 write(*,*) "Swap which information for the two GTFs?" 314 write(*,*) "1 Swap all propertie" 315 write(*,*) "2 Swap center" 316 write(*,*) "3 Swap function type" 317 write(*,*) "4 Swap exponent" 318 write(*,*) "5 Swap orbital expansion coefficient" 319 read(*,*) iswapcontent 320 if (iswapcontent==1) call swap(i,j,"all") 321 if (iswapcontent==2) call swap(i,j,"cen") 322 if (iswapcontent==3) call swap(i,j,"typ") 323 if (iswapcontent==4) call swap(i,j,"exp") 324 if (iswapcontent==5) call swap(i,j,"MO ") 325 write(*,*) "Swapping finished!" 326 327 else if (iselect==21) then 328 write(*,*) "Input the index of primitive function" 329 read(*,*) i 330 write(*,*) "Input the center" 331 read(*,*) j 332 if (j<=ncenter.and.j>0) then 333 b(i)%center=j 334 else 335 write(*,"('Error: The value should >0 and <=',i7)") ncenter 336 end if 337 338 else if (iselect==22) then 339 write(*,*) "Input the index of primitive function" 340 read(*,*) i 341 write(*,*) "Input the type" 342 write(*,*) "Valid input: S/X/Y/Z/XX/YY/ZZ/XY/XZ/YZ/XXX/YYY/ZZZ/XXY/XXZ/YYZ/XYY/XZZ/YZZ/XYZ" 343 write(*,*) "ZZZZ/YZZZ/YYZZ/YYYZ/YYYY/XZZZ/XYZZ/XYYZ/XYYY/XXZZ/XXYZ/XXYY/XXXZ/XXXY/XXXX" 344 write(*,"(a)") " ZZZZZ/YZZZZ/YYZZZ/YYYZZ/YYYYZ/YYYYY/XZZZZ/XYZZZ/XYYZZ/XYYYZ/XYYYY/XXZZZ/XXYZZ/XXYYZ/XXYYY/XXXZZ/XXXYZ/XXXYY/XXXXZ/XXXXY/XXXXX" 345 read(*,*) seltmpc 346 do j=1,size(GTFtype2name) 347 if (seltmpc==GTFtype2name(j)) then 348 b(i)%functype=j 349 exit 350 end if 351 if (j==20) write(*,*) "Error: Couldn't recognize this type" 352 end do 353 354 else if (iselect==23) then 355 write(*,*) "Input the index of primitive function" 356 read(*,*) i 357 write(*,*) "Input the exponent" 358 read(*,*) rexp 359 b(i)%exp=rexp 360 361 else if (iselect==24) then 362 write(*,*) "Input the index of primitive function" 363 read(*,*) iprm 364 write(*,*) "Input the orbital index, e.g. 12" 365 read(*,*) imonum 366 if (iprm<=nprims.and.iprm>0.and.imonum<=nmo.and.imonum>0) then 367 write(*,*) "Input the coefficient" 368 read(*,*) rcoeff 369 CO(imonum,iprm)=rcoeff 370 else 371 write(*,"('Error: The index of function or orbital exceed valid range')") 372 end if 373 374 else if (iselect==25) then 375 isetmode=1 376 if (allocated(CObasa)) then 377 write(*,*) "1 Set coefficients of GTFs" 378 write(*,*) "2 Set coefficients of basis functions" 379 read(*,*) isetmode 380 end if 381 if (isetmode==1) then 382 write(*,*) "The following your inputs are conditions for filtering GTFs" 383 write(*,*) "Rule of range input: 3,17 means from 3 to 17, 6,6 means only 6, 0,0 means all" 384 write(*,*) 385 write(*,*) "Input the range of index of GTFs" 386 read(*,*) ind1,ind2 387 write(*,*) "Input the range of atoms" 388 read(*,*) iatm1,iatm2 389 write(*,*) "Input the type of GTFs (one of S,X,Y,Z,XX,XY... ALL means all types)" 390 write(*,*) "You can also input S,P,D,F,G,H to select GTF according to angular moment" 391 read(*,*) seltmpc 392 write(*,*) "Input the range of orbitals" 393 read(*,*) imo1,imo2 394 write(*,*) "Input the expansion coefficient you want to set, e.g. 0.5" 395 read(*,*) coval 396 if (ind1==0) ind1=1 397 if (ind2==0) ind2=nprims 398 if (iatm1==0) iatm1=1 399 if (iatm2==0) iatm2=ncenter 400 if (imo1==0) imo1=1 401 if (imo2==0) imo2=nmo 402 nsel=0 403 do iGTF=ind1,ind2 404 if (b(iGTF)%center>=iatm1.and.b(iGTF)%center<=iatm2) then 405 if (seltmpc=="ALL".or.seltmpc=="all".or.GTFtype2name(b(iGTF)%functype)==trim(seltmpc).or.type2ang(b(iGTF)%functype)==trim(seltmpc)) then 406 CO(imo1:imo2,iGTF)=coval 407 nsel=nsel+1 408 end if 409 end if 410 end do 411 write(*,"(' Coefficient of',i8,' GTFs are set')") nsel 412 else if (isetmode==2) then 413 write(*,*) "The following your inputs are conditions for filtering basis functions" 414 write(*,*) "Rule of range input: 3,17 means from 3 to 17, 6,6 means only 6, 0,0 means all" 415 write(*,*) "You can also input S,P,D,F,G,H to select according to angular moment" 416 write(*,*) 417 write(*,*) "Input the range of index of basis functions" 418 read(*,*) ind1,ind2 419 write(*,*) "Input the range of atoms" 420 read(*,*) iatm1,iatm2 421 write(*,"(a)") " Input the type of basis functions (one of S,X,Y,Z,XX... ALL means all types)" 422 read(*,*) seltmpc 423 write(*,*) "Input the range of orbitals" 424 read(*,*) imo1,imo2 425 ispinsel=1 426 if (wfntype==1.or.wfntype==4) then 427 write(*,*) "For which type of orbitals? 0=Both 1=Alpha 2=Beta" 428 read(*,*) ispinsel 429 end if 430 write(*,*) "Input the expansion coefficient you want to set, e.g. 0.5" 431 read(*,*) coval 432 if (ind1==0) ind1=1 433 if (ind2==0) ind2=nbasis 434 if (iatm1==0) iatm1=1 435 if (iatm2==0) iatm2=ncenter 436 if (imo1==0) imo1=1 437 if (imo2==0) imo2=nbasis 438 nsel=0 439 do ibas=ind1,ind2 440 if (bascen(ibas)>=iatm1.and.bascen(ibas)<=iatm2) then 441 if (seltmpc=="ALL".or.seltmpc=="all".or.GTFtype2name(bastype(ibas))==trim(seltmpc).or.type2ang(bastype(ibas))==trim(seltmpc)) then 442 if (ispinsel==0.or.ispinsel==1) CObasa(ibas,imo1:imo2)=coval 443 if (ispinsel==0.or.ispinsel==2) CObasb(ibas,imo1:imo2)=coval 444 nsel=nsel+1 445 end if 446 end if 447 end do 448 write(*,"(' Coefficient of',i8,' basis functions are set')") nsel 449 imodwfn=1 450 end if 451 write(*,*) "Done!" 452 453 else if (iselect==26) then 454 do while(.true.) 455 write(*,*) "Select the orbitals for which the occupation numbers are needed to be changed" 456 write(*,*) "e.g. 2,4,13-16,20 means select orbital 2,4,13,14,15,16,20" 457 write(*,*) "Input 0 can select all orbitals, input q or 00 can return" 458 read(*,"(a)") c1000tmp 459 if (c1000tmp(1:1)=='q'.or.c1000tmp(1:2)=='00') exit 460 if (c1000tmp(1:1)=='0') then 461 numorbsel=nmo 462 do i=1,nmo 463 orbarr(i)=i 464 end do 465 else 466 call str2arr(c1000tmp,numorbsel,orbarr) 467 if ( any(orbarr(1:numorbsel)<1).or.any(orbarr(1:numorbsel)>nmo) ) then 468 write(*,*) "Error: One or more orbital indices exceeded valid range!" 469 cycle 470 end if 471 end if 472 write(*,*) "Set to which value? e.g. 1.2" 473 write(*,*) "Note:" 474 write(*,"(a)") " You can also input for example ""+1.1"" ""-1.1"" ""*1.1"" ""/1.1"" to add, minus, multiply and divide the occupation numbers by 1.1" 475 write(*,"(a)") " To recover the initial occupation numbers, input ""i""" 476 write(*,"(a)") " To generate occupation state for calculating odd electron density, input ""odd""" 477 read(*,"(a)") c1000tmp 478 if (index(c1000tmp,"odd")/=0) then 479 do iorb=1,numorbsel 480 MOocc(orbarr(iorb))=min(2-MOocc(orbarr(iorb)),MOocc(orbarr(iorb))) 481 end do 482 write(*,*) "Done!" 483 else if (index(c1000tmp,"i")/=0) then 484 MOocc(orbarr(1:numorbsel))=MOocc_org(orbarr(1:numorbsel)) 485 write(*,*) "The occupation numbers have been recovered" 486 else if (c1000tmp(1:1)=='+'.or.c1000tmp(1:1)=='-'.or.c1000tmp(1:1)=='*'.or.c1000tmp(1:1)=='/') then 487 read(c1000tmp(2:),*) tmpval 488 if (c1000tmp(1:1)=='+') MOocc(orbarr(1:numorbsel))=MOocc(orbarr(1:numorbsel))+tmpval 489 if (c1000tmp(1:1)=='-') MOocc(orbarr(1:numorbsel))=MOocc(orbarr(1:numorbsel))-tmpval 490 if (c1000tmp(1:1)=='*') MOocc(orbarr(1:numorbsel))=MOocc(orbarr(1:numorbsel))*tmpval 491 if (c1000tmp(1:1)=='/') MOocc(orbarr(1:numorbsel))=MOocc(orbarr(1:numorbsel))/tmpval 492 else 493 read(c1000tmp,*) tmpval 494 MOocc(orbarr(1:numorbsel))=tmpval 495 write(*,*) "Done!" 496 end if 497 call updatenelec !Update the number of electrons 498 if (occtmp==0D0) write(*,"(/,a)") " Note: When saving present wavefunction to new .wfn file, the orbitals with zero occupation number will be discarded" 499 imodwfn=1 500 if (any(MOocc/=int(MOocc))) then 501 if (wfntype==0) then 502 wfntype=3 !RHF-> Restricted post-HF wavefunction 503 write(*,*) "Note: Now the wavefunction is recognized as restricted post-HF wavefunction" 504 else if (wfntype==1.or.wfntype==2) then !UHF/ROHF-> Unrestricted post-HF wavefunction 505 wfntype=4 506 write(*,*) "Note: Now the wavefunction is recognized as unrestricted post-HF wavefunction" 507 end if 508 end if 509 write(*,*) 510 end do 511 512 else if (iselect==27) then 513 write(*,*) "Set type for which range of orbitals? e.g. 14,17" 514 write(*,*) "Hint: Input 0,0 can select all orbitals" 515 read(*,*) iorb1,iorb2 516 if (iorb1==0.and.iorb2==0) then 517 iorb1=1 518 iorb2=nmo 519 end if 520 write(*,*) "Set to which type? 0=Alpha&Beta 1=Alpha 2=Beta" 521 read(*,*) isettype 522 MOtype(iorb1:iorb2)=isettype 523 !Recount alpha and beta electrons 524 call updatenelec 525 write(*,*) "Done!" 526 !Update wavefunction type 527 if (all(MOtype==0)) then 528 if (all(MOocc==nint(MOocc))) then !All A+B orbital & integer occupation 529 wfntype=0 530 write(*,"('Note: Now the wavefunction is recognized as restricted close-shell single-determinant wavefunction')") 531 else !All A+B orbital & partial occupation 532 wfntype=3 533 write(*,"('Note: Now the wavefunction is recognized as close-shell post-HF wavefunction')") 534 end if 535 write(*,*) 536 else 537 if (any(MOocc/=nint(MOocc)).and.all(MOtype/=0)) then !Either A or B, and partial occupation 538 wfntype=4 539 write(*,"('Note: Now the wavefunction is recognized as open-shell post-HF wavefunction')") 540 else if (all(MOocc==nint(MOocc)).and.all(MOtype/=0)) then !Integer occupation and either A or B 541 wfntype=1 542 write(*,"('Note: Now the wavefunction is recognized as unrestricted single-determinant wavefunction')") 543 else if (all(MOocc==nint(MOocc)).and.any(MOtype==0).and.all(MOtype/=2)) then !Integer occupation and at least one orbital is A, and B is unexisted 544 wfntype=2 545 write(*,"('Note: Now the wavefunction is recognized as restricted open-shell wavefunction')") 546 else 547 write(*,"('Warning: The type of present wavefunction cannot be identified! You need to reset orbital types')") 548 end if 549 write(*,*) 550 end if 551 imodwfn=1 552 553 else if (iselect==31) then 554 write(*,*) "Input X,Y,Z of translation vector (e.g. 3.2,1.0,0)" 555 read(*,*) pbctransx,pbctransy,pbctransz 556 write(*,*) "You inputted coordinates are in which unit? 1:Bohr 2:Angstrom" 557 read(*,*) iunit 558 if (iunit==2) then 559 pbctransx=pbctransx/b2a 560 pbctransy=pbctransy/b2a 561 pbctransz=pbctransz/b2a 562 end if 563 do i=1,ncenter 564 a(i)%x=a(i)%x+pbctransx 565 a(i)%y=a(i)%y+pbctransy 566 a(i)%z=a(i)%z+pbctransz 567 end do 568 imodwfn=1 569 570 else if (iselect==32) then 571 write(*,*) "Input X,Y,Z of translation vector (e.g. 3.2,1.0,0)" 572 read(*,*) pbctransx,pbctransy,pbctransz 573 write(*,*) "You inputted coordinates are in which unit? 1:Bohr 2:Angstrom" 574 read(*,*) iunit 575 if (iunit==2) then 576 pbctransx=pbctransx/b2a 577 pbctransy=pbctransy/b2a 578 pbctransz=pbctransz/b2a 579 end if 580 write(*,*) "Duplicate system how many times?" 581 read(*,*) numdup 582 !_tmp is for backing up current information 583 allocate(a_tmp(ncenter)) 584 allocate(b_tmp(nprims)) 585 allocate(CO_tmp(nmo,nprims)) 586 a_tmp=a 587 b_tmp=b 588 CO_tmp=CO 589 deallocate(a,b,CO) 590 nprims_tmp=nprims 591 ncenter_tmp=ncenter 592 nprims=nprims*(numdup+1) 593 ncenter=ncenter*(numdup+1) 594 nelec=nelec*(numdup+1) 595 naelec=naelec*(numdup+1) 596 nbelec=nbelec*(numdup+1) 597 allocate(a(ncenter)) 598 allocate(b(nprims)) 599 allocate(CO(nmo,nprims)) 600 do idup=0,numdup 601 a(ncenter_tmp*idup+1:ncenter_tmp*(idup+1))=a_tmp(1:center_tmp) 602 a(ncenter_tmp*idup+1:ncenter_tmp*(idup+1))%x=a_tmp(1:center_tmp)%x+pbctransx*idup 603 a(ncenter_tmp*idup+1:ncenter_tmp*(idup+1))%y=a_tmp(1:center_tmp)%y+pbctransy*idup 604 a(ncenter_tmp*idup+1:ncenter_tmp*(idup+1))%z=a_tmp(1:center_tmp)%z+pbctransz*idup 605 b(nprims_tmp*idup+1:nprims_tmp*(idup+1))=b_tmp(1:nprims_tmp) 606 b(nprims_tmp*idup+1:nprims_tmp*(idup+1))%center=b_tmp(1:nprims_tmp)%center+ncenter_tmp*idup 607 CO(:,nprims_tmp*idup+1:nprims_tmp*(idup+1))=CO_tmp(:,1:nprims_tmp) !Notice that the orbitals do not satisify normalization condition any more, and the orbital occupation number will be artifical 608 end do 609 deallocate(a_tmp,b_tmp,CO_tmp) 610 imodwfn=1 611 call gendistmat !The number of atoms have changed, so we must update distance matrix 612 613 else if (iselect==33) then 614 write(*,*) "Rotate which orbital? (Input 0 to rotate all orbitals)" 615 read(*,*) iorb 616 if (iorb/=0) then 617 call orbcoeffrotate(iorb) 618 else if (iorb==0) then 619 do imo=1,nmo 620 call orbcoeffrotate(imo) 621 end do 622 write(*,*) "Also rotate atomic coordinates? (y/n)" 623 read(*,*) selectyn 624 if (selectyn=='y'.or.selectyn=='Y') then 625 do iatm=1,ncenter 626 tmpval=a(iatm)%x 627 a(iatm)%x=a(iatm)%z 628 a(iatm)%z=a(iatm)%y 629 a(iatm)%y=tmpval 630 end do 631 end if 632 end if 633 write(*,*) "Done!" 634 635 else if (iselect==34) then 636 innerel=0 637 do i=1,ncenter 638 if (int(a(i)%charge)/=a(i)%index) then 639 write(*,"(' Note: Atom',i5,' is not taken into account since it utilizes pseudopotential')") i 640 cycle 641 end if 642 if (a(i)%index>2.and.a(i)%index<=10) innerel=innerel+2 643 if (a(i)%index>10.and.a(i)%index<=18) innerel=innerel+10 644 if (a(i)%index>18.and.a(i)%index<=36) innerel=innerel+18 645 if (a(i)%index>36.and.a(i)%index<=54) innerel=innerel+36 646 if (a(i)%index>54.and.a(i)%index<=86) innerel=innerel+54 647 if (a(i)%index>86) innerel=innerel+86 648 end do 649 nelec=nelec-innerel 650 naelec=naelec-innerel/2 651 nbelec=nbelec-innerel/2 652 if (wfntype==1.or.wfntype==4) then !UHF and U-post-HF wfn 653 MOocc(1:innerel/2)=0D0 654 do j=1,nmo !Where the first beta orbital appear now 655 if (motype(j)==2) exit 656 end do 657 MOocc(1:j+innerel/2-1)=0D0 658 write(*,"(' The occupation of',i7,' lowest energy orbitals have been set to zero')") innerel 659 else if (wfntype==0.or.wfntype==2.or.wfntype==3) then !restricted(=0) or RO(=2) or post-R(=3) wavefunction 660 MOocc(1:innerel/2)=0D0 661 write(*,"(' The occupation of',i7,' lowest energy orbitals have been set to zero')") innerel/2 662 end if 663 if (wfntype==3.or.wfntype==4) write(*,"(' Warning: Discarding inner orbitals for post-HF wavefunction will lead to unexpected result!')") 664 imodwfn=1 665 666 else if (iselect==35) then 667 call SelMO_IRREP 668 669 else if (iselect==36) then 670 write(*,*) "Input index of the orbitals, e.g. 2,3,7-10" 671 read(*,"(a)") c1000tmp 672 call str2arr(c1000tmp,ntmp,orbarr) 673 do idxtmp=1,ntmp 674 idx=orbarr(idxtmp) 675 CO(idx,:)=-CO(idx,:) 676 if (allocated(CObasa)) then 677 if (idx<=nbasis) then 678 CObasa(:,idx)=-CObasa(:,idx) 679 else 680 CObasb(:,idx-nbasis)=-CObasb(:,idx-nbasis) 681 end if 682 end if 683 end do 684 write(*,*) "Done!" 685 imodwfn=1 686 end if 687 write(*,*) 688end do 689end subroutine 690 691 692 693 694!!---------------- Select MOs according to irreducible representation 695subroutine SelMO_IRREP 696use defvar 697use util 698character symlab(nmo)*4,c2000tmp*2000,symstat(nmo)*9 !Allocate the array lengths as upper limit 699integer tmparr(nmo),symNorb(nmo) !Allocate the array lengths as upper limit 700if (wfntype/=0.and.wfntype/=1) then 701 write(*,"(a)") " Error: This function only works for RHF or UHF wavefunctions (or the DFT counterparts)" 702 return 703end if 704nsym=0 705do imo=1,nmo 706 if (MOocc_org(imo)==0D0) cycle 707 if (all(symlab(1:nsym)/=MOsym(imo))) then 708 nsym=nsym+1 709 symlab(nsym)=MOsym(imo) 710 end if 711end do 712symNorb=0 713do imo=1,nmo 714 if (MOocc_org(imo)==0D0) cycle 715 do isym=1,nsym 716 if (MOsym(imo)==symlab(isym)) symNorb(isym)=symNorb(isym)+1 717 end do 718end do 719symstat="Normal" 720MOocc=MOocc_org 721if (imodwfn==1) write(*,*) "Note: Original occupation status has been recovered" 722write(*,*) "Note: Only the orbitals that originally occupied are taken into account here" 723do while(.true.) 724 write(*,*) 725 write(*,*) "Information of various irreducible representations:" 726 do isym=1,nsym 727 write(*,"(i5,' Sym: ',a,' N_orb:',i5,' Status: ',a)") isym,symlab(isym),symNorb(isym),symstat(isym) 728 end do 729 write(*,*) 730 write(*,*) "0 Save and return" 731 write(*,*) "1 Discard specific irreducible representations" 732 write(*,*) "2 Recover original status" 733 write(*,*) "3 Reverse status" 734 read(*,*) isel 735 736 if (isel==0) then 737 call updatenelec 738 imodwfn=1 739 write(*,*) "The current orbital occupation status has been saved" 740 write(*,*) "Updating density matrix..." 741 call genP 742 write(*,*) "Density matrix has been updated" 743 exit 744 else if (isel==2) then 745 MOocc=MOocc_org 746 symstat="Normal" 747 else if (isel==1.or.isel==3) then 748 if (isel==1) then 749 write(*,*) "Input the index of the irreducible representations to be discarded, e.g. 1,3-5" 750 read(*,"(a)") c2000tmp 751 call str2arr(c2000tmp,nsymsel,tmparr) 752 do isym=1,nsymsel 753 symstat(tmparr(isym))="Discarded" 754 end do 755 else if (isel==3) then 756 do isym=1,nsym 757 if (symstat(isym)=="Normal") then 758 symstat(isym)="Discarded" 759 else 760 symstat(isym)="Normal" 761 end if 762 end do 763 end if 764 do imo=1,nmo 765 if (MOocc_org(imo)==0D0) cycle 766 do isym=1,nsym 767 if (MOsym(imo)==symlab(isym)) then 768 if (symstat(isym)=="Normal") MOocc(imo)=MOocc_org(imo) 769 if (symstat(isym)=="Discarded") MOocc(imo)=0D0 770 exit 771 end if 772 end do 773 end do 774 end if 775 write(*,*) "Done!" 776end do 777end subroutine 778 779 780!!---------- Update the number of electrons 781subroutine updatenelec 782use defvar 783integer imo 784nelec=0 785naelec=0 786nbelec=0 787do imo=1,nmo 788 if (MOtype(imo)==0) then 789 naelec=naelec+MOocc(imo)/2D0 790 nbelec=nbelec+MOocc(imo)/2D0 791 else if (MOtype(imo)==1) then 792 naelec=naelec+MOocc(imo) 793 else if (MOtype(imo)==2) then 794 nbelec=nbelec+MOocc(imo) 795 end if 796end do 797nelec=naelec+nbelec 798end subroutine 799 800 801!!!-------- Check if present wavefunction is sanity, i.e. all orbital satisfies normalization condition 802subroutine wfnsanity 803use defvar 804implicit real*8 (a-h,o-z) 805real*8 GTFSmat(nprims*(nprims+1)/2) 806call genGTFSmat(GTFSmat,nprims*(nprims+1)/2) 807rmaxdev=0 808do imo=1,nmo 809 tmp=0 810 do iGTF=1,nprims 811 do jGTF=iGTF+1,nprims 812 tmp=tmp+2*co(imo,iGTF)*co(imo,jGTF)*GTFSmat(jGTF*(jGTF-1)/2+iGTF) 813 end do 814 tmp=tmp+co(imo,iGTF)**2*GTFSmat(iGTF*(iGTF-1)/2+iGTF) 815 end do 816 write(*,"(' Orbital',i7,', Occ:',f8.4,' Value:',f16.10)") imo,MOocc(imo),tmp 817 if (tmp>rmaxdev) rmaxdev=tmp 818end do 819write(*,"(' Maximum deviation to 1:',f16.10)") rmaxdev-1 820end subroutine 821 822 823 824!!---------- Return normalization coefficient for specific type of cartesian type GTF, see Levine 5ed p487 825!The meaning of itype is defined in GTFtype2name 826real*8 function normgau(itype,exp) 827use defvar 828use util 829implicit real*8 (a-h,o-z) 830ix=type2ix(itype) 831iy=type2iy(itype) 832iz=type2iz(itype) 833normgau=(2*exp/pi)**0.75D0*dsqrt( (8*exp)**(ix+iy+iz)*ft(ix)*ft(iy)*ft(iz)/(ft(2*ix)*ft(2*iy)*ft(2*iz)) ) 834end function 835!!!---------- Return normalization coefficient for specific type of spherical harmonic GTF 836!See http://en.wikipedia.org/wiki/Gaussian_orbital for the formula 837!!!This is useless for resolving the contraction coefficent problem of Molden input file of ORCA , I don't know why 838!Lval is the angular moment, 0/1/2/3/4=s/p/d/f/g 839! real*8 function rnormgau_sph(Lval,exp) 840! use defvar 841! use util 842! integer Lval 843! real*8 exp 844! rnormgau_sph=exp**(Lval/2D0+0.75D0) / (dsqrt(pi)*2**(-0.25D0-Lval/2D0)) / dsqrt(gamma_ps(Lval+1)) 845! end function 846 847subroutine genn1n2nf(Lval,n1,n2,nf) 848integer Lval,n1,n2,nf 849if (Lval==0) then 850 n1=3 851 n2=3 852 nf=1 853else if (Lval==1) then 854 n1=7 855 n2=5 856 nf=1 857else if (Lval==2) then 858 n1=11 859 n2=7 860 nf=9 861else if (Lval==3) then 862 n1=15 863 n2=9 864 nf=225 865else if (Lval==4) then 866 n1=19 867 n2=11 868 nf=11025 869end if 870end subroutine 871!!---- Used to produce normalization factor to counteract the contraction coefficient problem of the Molden input file generated by ORCA 872!I don't know where the formula comes from 873!Lval=0/1/2/3/4 corresponds to s/p/d/f/g 874real*8 function rnormgau_ORCA(exp,Lval) 875real*8 exp,pi 876integer Lval,n1,n2,nf 877pi=acos(-1D0) 878call genn1n2nf(Lval,n1,n2,nf) 879rnormgau_ORCA=dsqrt(dsqrt(2**n1*exp**n2/(pi**3*nf*nf))) 880end function 881!----- Renormalization of Gauss basis functions for Molden input file 882subroutine renormmoldengau(nlen,Lval,exp,con) 883implicit real*8 (a-h,o-z) 884integer nlen,Lval 885real*8 exp(nlen),con(nlen),ctmp(nlen) 886pi=acos(-1D0) 887call genn1n2nf(Lval,n1,n2,nf) 888fc=(2D0**n1)/(pi**3*nf) 889do i=1,nlen 890 prmnormfac=sqrt(sqrt(fc*(exp(i)**n2))) 891 ctmp(i)=con(i)*prmnormfac 892end do 893facnorm=0D0 894do i=1,nlen 895 do j=1,i 896 expavg=(exp(i)+exp(j))/2D0 897 facadd=ctmp(i)*ctmp(j)/sqrt(fc*(expavg**n2)) 898 if (i/=j) facadd=facadd*2 899 facnorm=facnorm+facadd 900 end do 901end do 902if (facnorm>1D-10) facnorm=1/sqrt(facnorm) 903con=con*facnorm 904end subroutine 905 906 907!!----- Use Lowdin orthogonalization method to transform density matrix (Xmat) and coefficient to ortho-normalized basis, meanwhile update Sbas to identity matrix 908!see Szabo p143 909subroutine symmortho 910use defvar 911use util 912implicit real*8 (a-h,o-z) 913real*8 Umat(nbasis,nbasis),svalvec(nbasis),Xmat(nbasis,nbasis) !workvec(3*nbasis-1) 914! call DSYEV('V','U',nbasis,sbas,nbasis,svalvec,workvec,3*nbasis-1,ierror) !lapack, the resultant sbas is eigenvector matrix 915! call diagmat(Sbas,Umat,svalvec) !My slow diagonalization routine 916call diagsymat(Sbas,Umat,svalvec,ierror) 917if (ierror/=0) write(*,*) "Error: Diagonalization of overlap matrix failed!" 918!Now Sbas is diagonalized matrix 919forall (i=1:nbasis) Sbas(i,i)=dsqrt(svalvec(i)) 920Xmat=matmul(matmul(Umat,Sbas),transpose(Umat)) !Then Xmat is S^0.5 921Ptot=matmul(matmul(Xmat,Ptot),Xmat) 922if (allocated(Palpha)) then 923 Palpha=matmul(matmul(Xmat,Palpha),Xmat) 924 Pbeta=Ptot-Palpha 925end if 926Cobasa=matmul(Xmat,Cobasa) 927if (allocated(Cobasb)) Cobasb=matmul(Xmat,Cobasb) 928forall(i=1:nbasis) Sbas(i,i)=1D0 !Reconstruct overlap matrix in ortho basis function 929end subroutine 930!!----- Like symmortho, but input overlap matrix, and only output Lowdin orthogonalization transformation matrix Xmat 931! innbas=input number of basis. Smatin is input overlap matrix, which will be modified 932! inv=1 get S^-0.5, other get S^0.5 933subroutine symmorthomat(innbas,Smatin,Xmat,inv) 934use defvar 935use util 936integer inv 937real*8 Umat(innbas,innbas),svalvec(innbas),Smatin(nbasis,nbasis),Smat(innbas,innbas),Xmat(innbas,innbas) 938Smat=Smatin 939call diagsymat(Smat,Umat,svalvec,ierror) 940forall (i=1:nbasis) Smat(i,i)=dsqrt(svalvec(i)) 941if (inv==1) forall (i=1:nbasis) Smat(i,i)=1D0/Smat(i,i) 942if (ierror/=0) write(*,*) "Error: Diagonalization of overlap matrix is fail!" 943Xmat=matmul(matmul(Umat,Smat),transpose(Umat)) 944end subroutine 945 946 947!!!------------------------- Generate distance matrix 948subroutine gendistmat 949use defvar 950implicit real*8 (a-h,o-z) 951if (allocated(distmat)) deallocate(distmat) 952allocate(distmat(ncenter,ncenter)) 953distmat=0.0D0 954do i=1,ncenter 955 do j=i+1,ncenter 956 distmat(i,j)=dsqrt((a(i)%x-a(j)%x)**2+(a(i)%y-a(j)%y)**2+(a(i)%z-a(j)%z)**2) 957 end do 958end do 959distmat=distmat+transpose(distmat) 960end subroutine 961 962 963!!!------------------------- Swap two gaussian primitive functions 964subroutine swap(i,j,swaptype) 965use defvar 966integer n,i,j 967character*3 swaptype 968type(primtype) tempb !For exchanging basis functions' order 969if (swaptype=="all") then 970 tempb=b(i) 971 b(i)=b(j) 972 b(j)=tempb 973else if (swaptype=="cen") then 974 tempb%center=b(i)%center 975 b(i)%center=b(j)%center 976 b(j)%center=tempb%center 977else if (swaptype=="typ") then 978 tempb%functype=b(i)%functype 979 b(i)%functype=b(j)%functype 980 b(j)%functype=tempb%functype 981else if (swaptype=="exp") then 982 tempb%exp=b(i)%exp 983 b(i)%exp=b(j)%exp 984 b(j)%exp=tempb%exp 985end if 986if (swaptype=="all".or.swaptype=="MO ") then 987 do n=1,nmo 988 temp=co(n,i) 989 co(n,i)=co(n,j) 990 co(n,j)=temp 991 end do 992end if 993end subroutine 994 995 996!!!---- Rotate(exchange) GTF and basis function coefficients within all shell in different direction of specific orbital 997! use this three times, namely XYZ->ZXY->YZX->XYZ, the coefficient recovered. 998! In detail, for examples, for d-type will lead to such coefficient exchange: XX to YY, YY to ZZ, ZZ to XX, XY to YZ, XZ to XY, YZ to XZ 999! exchange only involve the GTFs/basis func. in the same shell 1000subroutine orbcoeffrotate(orb) !orb=Rotate which orbital 1001use defvar 1002implicit real*8 (a-h,o-z) 1003integer orb 1004real*8 COorborg(nprims) !For backing up origin CO 1005real*8 CObasa_tmp(nbasis) !For backing up origin CObasa 1006real*8 CObasb_tmp(nbasis) !For backing up origin CObasb 1007COorborg(:)=CO(orb,:) 1008do i=1,nprims 1009 ixtmp=type2iz(b(i)%functype) 1010 iytmp=type2ix(b(i)%functype) 1011 iztmp=type2iy(b(i)%functype) 1012 do j=1,nprims 1013 if (type2ix(b(j)%functype)==ixtmp.and.type2iy(b(j)%functype)==iytmp.and.& 1014 type2iz(b(j)%functype)==iztmp.and.b(j)%exp==b(i)%exp.and.b(j)%center==b(i)%center) CO(orb,j)=COorborg(i) 1015 end do 1016end do 1017if (allocated(CObasa)) then 1018 CObasa_tmp=CObasa(:,orb) 1019 if (wfntype==1.or.wfntype==4) CObasb_tmp=CObasb(:,orb) 1020 do iatm=1,ncenter 1021 do ibas=basstart(iatm),basend(iatm) 1022 ityp=bastype(ibas) 1023 ixtmp=type2iz(ityp) 1024 iytmp=type2ix(ityp) 1025 iztmp=type2iy(ityp) 1026 do jbas=basstart(iatm),basend(iatm) 1027 jtyp=bastype(jbas) 1028 if (type2ix(jtyp)==ixtmp.and.type2iy(jtyp)==iytmp.and.type2iz(jtyp)==iztmp.and.& 1029 basshell(ibas)==basshell(jbas)) then 1030 CObasa(jbas,orb)=CObasa_tmp(ibas) 1031 if (wfntype==1.or.wfntype==4) CObasb(jbas,orb)=CObasb_tmp(ibas) 1032 end if 1033 end do 1034 end do 1035 end do 1036end if 1037end subroutine 1038 1039 1040!!!------ Define property/origin/spacing/grid number and then save to a 3D matrix, infomode=1 means silent 1041!! iorb is used to choose the orbital for whose wavefunction will be calculated. This can be an arbitrary value if functype/=4 1042subroutine savecubmat(functype,infomode,iorb) 1043use defvar 1044use util 1045use function 1046implicit none 1047integer walltime1,walltime2 1048integer :: i,j,k,ii,infomode,functype,calcfunc,ifinish,iorb !Calculate which orbital wavefunction for fmo routine 1049real*8 t,time_begin,time_end,time_endtmp,tmpx,tmpy,tmpz,tmprho,xarr(nx),yarr(ny),zarr(nz) 1050character c80tmp*80 1051iorbsel=iorb 1052calcfunc=functype 1053if (functype==12) calcfunc=8 !If calculate total ESP, first calculate nuclear ESP, and finally call espcub 1054if (functype==100.and.iuserfunc==14) calcfunc=0 !If calculate electronic ESP, first skip grid calculation, and finally call espcub 1055if (infomode==0.and.functype/=12) then 1056 if (expcutoff<0) write(*,"(' Note: All exponential functions exp(x) with x<',f8.3,' will be ignored ')") expcutoff 1057end if 1058 1059ii=10 1060ifinish=0 1061!Cut the minimal noise at the end of the coordinate, otherwise the originally symmetry points may become unsymmetry 1062do k=1,nz 1063 write(c80tmp,"(D20.13)") orgz+(k-1)*dz 1064 read(c80tmp,*) zarr(k) 1065end do 1066do j=1,ny 1067 write(c80tmp,"(D20.13)") orgy+(j-1)*dy 1068 read(c80tmp,*) yarr(j) 1069end do 1070do i=1,nx 1071 write(c80tmp,"(D20.13)") orgx+(i-1)*dx 1072 read(c80tmp,*) xarr(i) 1073end do 1074 1075call walltime(walltime1) 1076CALL CPU_TIME(time_begin) 1077nthreads=getNThreads() 1078!$OMP PARALLEL DO SHARED(cubmat,ifinish) PRIVATE(i,j,k,tmpx,tmpy,tmpz,tmprho) schedule(dynamic) NUM_THREADS(nthreads) 1079do k=1,nz 1080 tmpz=zarr(k) 1081 do j=1,ny 1082 tmpy=yarr(j) 1083 do i=1,nx 1084 tmpx=xarr(i) 1085 if (calcfunc==1513) then !Only involved by funcvsfunc routine, when RDG and signlambda2rho is combined 1086 call signlambda2rho_RDG(tmpx,tmpy,tmpz,cubmat(i,j,k),cubmattmp(i,j,k),tmprho) 1087 else if (calcfunc==1614) then !Only involved by funcvsfunc routine, when RDG and signlambda2rho is combined 1088 call signlambda2rho_RDG_prodens(tmpx,tmpy,tmpz,cubmat(i,j,k),cubmattmp(i,j,k)) 1089 else 1090 cubmat(i,j,k)=calcfuncall(calcfunc,tmpx,tmpy,tmpz) 1091 end if 1092 end do 1093 end do 1094 if (infomode==0.and.functype/=12) then 1095 if ( nthreads ==1) then 1096 CALL CPU_TIME(time_endtmp) 1097 t=dfloat(k)/nz*100 !Completed percent 1098 if (k==1) then !Show approximate time at start 1099 t=(time_endtmp-time_begin)*(nz-1) !How many time to work out remain work 1100 write(*,"(' Calculation will take up CPU time about',f9.2,' seconds (',f8.2,' minutes)')") t,t/60D0 1101 else if (t>=ii) then 1102 ii=ii+10 1103 write(*,"(f5.1,'% completed,',f8.1,' seconds remain')") t,(time_endtmp-time_begin)/k*(nz-k) 1104 end if 1105 else 1106 ifinish=ifinish+1 1107 write(*,"(' Finished:',i5,' /',i5)") ifinish,nz 1108 end if 1109 end if 1110end do 1111!$OMP END PARALLEL DO 1112CALL CPU_TIME(time_end) 1113call walltime(walltime2) 1114if (functype==12.or.(functype==100.and.iuserfunc==14)) then 1115 call espcub !Calculate electronic ESP and accumulate to existing cubmat 1116else 1117 if (infomode==0) write(*,"(' Calculation took up CPU time',f12.2,'s, wall clock time',i10,'s')") time_end-time_begin,walltime2-walltime1 1118end if 1119end subroutine 1120 1121 1122!!!------ Output molecular formula 1123subroutine showformula 1124use defvar 1125implicit real*8 (a-h,o-z) 1126character*6 tmp 1127write(*,"(' Formula: ')",advance="no") 1128do i=0,nelesupp 1129 n=0 1130 do iatm=1,ncenter 1131 if (a(iatm)%index==i) n=n+1 1132 end do 1133 write(tmp,"(i6)") n 1134 if (n/=0) write(*,"(a,a,' ')",advance="no") trim(ind2name(i)),trim(adjustl(tmp)) 1135end do 1136write(*,*) 1137end subroutine 1138 1139 1140!!!----------- Resize number of orbitals of CO, MOene, MOtype, MOocc to "newnmo", also resize number of GTFs of CO to "newnprims" 1141subroutine resizebynmo(newnmo,newnprims) 1142use defvar 1143implicit real*8 (a-h,o-z) 1144real*8,allocatable :: CO_bk(:,:),MOene_bk(:),MOocc_bk(:) 1145integer,allocatable :: MOtype_bk(:) 1146integer newnmo,oldnmo,newnprims,oldnprims 1147oldnmo=size(CO,1) 1148oldnprims=size(CO,2) 1149allocate(CO_bk(oldnmo,oldnprims),MOene_bk(oldnmo),MOocc_bk(oldnmo),MOtype_bk(oldnmo)) 1150CO_bk=CO 1151MOene_bk=MOene 1152MOocc_bk=MOocc 1153MOtype_bk=MOtype 1154deallocate(CO,MOene,MOocc,MOtype) 1155allocate(CO(newnmo,newnprims),MOene(newnmo),MOocc(newnmo),MOtype(newnmo)) 1156if (newnmo>=oldnmo) then !Enlarge array size, don't forget to fill the gap afterwards 1157 if (newnprims>=oldnprims) CO(1:oldnmo,1:oldnprims)=CO_bk(:,:) 1158 if (newnprims<oldnprims) CO(1:oldnmo,:)=CO_bk(:,1:newnprims) 1159 MOene(1:oldnmo)=MOene_bk(:) 1160 MOocc(1:oldnmo)=MOocc_bk(:) 1161 MOtype(1:oldnmo)=MOtype_bk(:) 1162else if (newnmo<oldnmo) then !Reduce array size 1163 if (newnprims>=oldnprims) CO(:,1:oldnprims)=CO_bk(1:newnmo,:) 1164 if (newnprims<oldnprims) CO(:,:)=CO_bk(1:newnmo,1:newnprims) 1165 MOene(:)=MOene_bk(1:newnmo) 1166 MOocc(:)=MOocc_bk(1:newnmo) 1167 MOtype(:)=MOtype_bk(1:newnmo) 1168end if 1169deallocate(CO_bk,MOene_bk,MOocc_bk,MOtype_bk) 1170end subroutine 1171 1172 1173!!!------------ Generate gjf of atoms in molecule, and invoke Gaussian to get .wfn, then input them into custom list 1174subroutine setpromol 1175use defvar 1176use util 1177implicit real*8 (a-h,o-z) 1178integer :: i,j,k,itype=0 1179character*2 typename(100),nametmp 1180character*80 basisset,tmpdir,c80tmp 1181character*80 outwfnname 1182logical alivegauout,alivewfntmp,aliveatomwfn 1183 1184!The only difference between c80tmp and tmpdir is that the latter has \ or / separator at the end 1185if (iwfntmptype==1) then 1186 if (isys==1) tmpdir="wfntmp\" 1187 if (isys==2.or.isys==3) tmpdir="wfntmp/" 1188 c80tmp="wfntmp" 1189 inquire(file='./wfntmp/.',exist=alivewfntmp) 1190 if ((isys==1).and.(alivewfntmp.eqv..true.)) then !delete old wfntmp folder 1191 write(*,*) "Running: rmdir /S /Q wfntmp" 1192 call system("rmdir /S /Q wfntmp") 1193 else if ((isys==2.or.isys==3).and.alivewfntmp.eqv..true.) then 1194 write(*,*) "Running: rm -rf wfntmp" 1195 call system("rm -rf wfntmp") 1196 end if 1197else if (iwfntmptype==2) then 1198 do i=1,9999 !Find a proper name of temporary folder 1199 write(c80tmp,"('wfntmp',i4.4)") i 1200 inquire(file='./c80tmp/.',exist=alivewfntmp) 1201 if (alivewfntmp.eqv..false.) exit 1202 end do 1203 if (isys==1) write(tmpdir,"('wfntmp',i4.4,'\')") i 1204 if (isys==2.or.isys==3) write(tmpdir,"('wfntmp',i4.4,'/')") i 1205end if 1206write(*,*) "Running: mkdir "//trim(c80tmp) !Build new temporary folder 1207call system("mkdir "//trim(c80tmp)) 1208inquire(file="./atomwfn/.",exist=aliveatomwfn) 1209if (isys==1.and.aliveatomwfn.eqv..true.) then 1210 write(*,*) "Running: copy atomwfn\*.wfn "//trim(tmpdir) 1211 call system("copy atomwfn\*.wfn "//trim(tmpdir)) 1212else if ((isys==2.or.isys==3).and.aliveatomwfn.eqv..true.) then 1213 write(*,*) "Running: cp atomwfn/*.wfn "//trim(tmpdir) 1214 call system("cp atomwfn/*.wfn "//trim(tmpdir)) 1215end if 1216 1217noatmwfn=0 !Check if the atomic wfn file have pre-stored in atomwfn folder, if not, invoke gaussian to calc it 1218do i=1,nfragatmnum 1219 if (isys==1) inquire(file="atomwfn\"//a(fragatm(i))%name//".wfn",exist=alive) 1220 if (isys==2.or.isys==3) inquire(file="atomwfn/"//a(fragatm(i))%name//".wfn",exist=alive) 1221 if (.not.alive) then 1222 noatmwfn=1 1223 exit 1224 end if 1225end do 1226 1227if (noatmwfn==0) then 1228 write(*,"(a)") " All atom .wfn files needed have already presented in ""atomwfn"" folder, we will not calculate them" 1229else if (noatmwfn==1) then !Some or all atomic wfn don't exist, calc them 1230 !Select calculation level 1231 write(*,"(a)") " Note: Some or all atom .wfn files needed are not present in ""atomwfn"" folder, they must be calculated now. See Section 3.7.3 of the manual for detail." 1232 write(*,"(a)") " Now please input the level for calculating atom wfn files, theoretical method is optional." 1233 write(*,"(a)") " For example: 6-31G* or B3LYP/cc-pVDZ You can also add other keywords at the same time, e.g. M052X/6-311G(2df,2p) scf=(vshift=500) int=ultrafine" 1234 read(*,"(a)") basisset !Note: 6d 10f is not required for generating wfn files, since the work has been done in L607 internally 1235 !Check Gaussian path 1236 inquire(file=gaupath,exist=alive) 1237 if (.not.alive) then 1238 write(*,*) "Couldn't find Gaussian path defined in ""gaupath"" variable in settings.ini" 1239 if (isys==1) write(*,*) "Input the path of Gaussian executable file, e.g. ""d:\study\g03w\g03.exe""" 1240 if (isys==2.or.isys==3) write(*,*) "Input the path of Gaussian executable file, e.g. ""/sob/g09/g09""" 1241 do while(.true.) 1242 read(*,"(a)") gaupath 1243 inquire(file=gaupath,exist=alive) 1244 if (alive) exit 1245 write(*,*) "Couldn't find Gaussian executable file, input again" 1246 end do 1247 end if 1248end if 1249 1250!Generate .gjf file for all elements, regardless if their wfn file have already presented, meanwhile count the total number of elements 1251itype=0 1252do i=1,nfragatmnum 1253 inquire(file=trim(tmpdir)//a(fragatm(i))%name//".gjf",exist=alive) 1254 if (.not.alive) then 1255 itype=itype+1 !How many different types 1256 typename(itype)=a(fragatm(i))%name 1257 1258 if (a_org(fragatm(i))%index>36) then 1259 inquire(file=trim(tmpdir)//a(fragatm(i))%name//".wfn",exist=alive) 1260 if (.not.alive) then !The wfn file of the heavy element hasn't been provided in "atomwfn" and hence cannot be found in "wfntmp" here 1261 write(*,"(a,a,a)") " Error: Multiwfn can't invoke Gaussian to generate wavefunction file and sphericalize density for ",a(fragatm(i))%name,", since its & 1262 index is larger than 36! You should provide corresponding atom .wfn files in ""atomwfn"" folder manually" 1263 write(*,*) "Press ENTER to continue" 1264 read(*,*) 1265 return 1266 end if 1267 end if 1268 1269 open(14,file=trim(tmpdir)//a(fragatm(i))%name//".gjf",status="replace") 1270 !If user inputted including "/" e.g. B3LYP/6-31g*, will replace default theoretical method 1271 if (index(basisset,'/')==0) then 1272 if (a(fragatm(i))%index<=20.or.a(fragatm(i))%index>=31) then 1273 write(14,"(a,/)") "#T out=wfn ROHF/"//trim(basisset) !Main group elements. If not use scf=sp, in g09, RO calculations for IIIA elements are to converge 1274 write(14,"(a,/)") "Temporary file for promolecule, ROHF"//trim(basisset) 1275 else 1276 write(14,"(a,/)") "#T out=wfn UB3LYP/"//trim(basisset) !Transition metals 1277 write(14,"(a,/)") "Temporary file for promolecule, UB3LYP"//trim(basisset) 1278 end if 1279 else 1280 if (a(fragatm(i))%index<=20.or.a(fragatm(i))%index>=31) then 1281 write(14,"(a,/)") "#T out=wfn RO"//trim(basisset) !Main group elements 1282 write(14,"(a,/)") "Temporary file for promolecule, RO"//trim(basisset) 1283 else 1284 write(14,"(a,/)") "#T out=wfn U"//trim(basisset) !Transition metals (RO may leads to convergence problem) 1285 write(14,"(a,/)") "Temporary file for promolecule, U"//trim(basisset) 1286 end if 1287 end if 1288 1289 !Currently support up to the fourth row 1290 if (a(fragatm(i))%name=="H ".or.a(fragatm(i))%name=="Li".or.a(fragatm(i))%name=="Na".or.a(fragatm(i))%name=="K") write(14,*) "0 2" 1291 if (a(fragatm(i))%name=="Be".or.a(fragatm(i))%name=="Mg".or.a(fragatm(i))%name=="Ca") write(14,*) "0 1" 1292 if (a(fragatm(i))%name=="B ".or.a(fragatm(i))%name=="Al".or.a(fragatm(i))%name=="Ga") write(14,*) "0 2" 1293 if (a(fragatm(i))%name=="C ".or.a(fragatm(i))%name=="Si".or.a(fragatm(i))%name=="Ge") then 1294 if (SpherIVgroup==0) write(14,*) "0 5" 1295 if (SpherIVgroup==1) write(14,*) "0 3" 1296 end if 1297 if (a(fragatm(i))%name=="N ".or.a(fragatm(i))%name=="P ".or.a(fragatm(i))%name=="As") write(14,*) "0 4" 1298 if (a(fragatm(i))%name=="O ".or.a(fragatm(i))%name=="S ".or.a(fragatm(i))%name=="Se") write(14,*) "0 3" 1299 if (a(fragatm(i))%name=="F ".or.a(fragatm(i))%name=="Cl".or.a(fragatm(i))%name=="Br") write(14,*) "0 2" 1300 if (a(fragatm(i))%name=="He".or.a(fragatm(i))%name=="Ne".or.a(fragatm(i))%name=="Ar".or.a(fragatm(i))%name=="Kr") write(14,*) "0 1" 1301 if (a(fragatm(i))%name=="Sc") write(14,*) "0 2" !3d1 4s2 1302 if (a(fragatm(i))%name=="Ti") write(14,*) "0 3" !3d2 4s2 1303 if (a(fragatm(i))%name=="V ") write(14,*) "0 4" !3d3 4s2 1304 if (a(fragatm(i))%name=="Cr") write(14,*) "0 7" !3d5 4s1, needn't sphericalization 1305 if (a(fragatm(i))%name=="Mn") write(14,*) "0 6" !3d5 4s2, needn't sphericalization 1306 if (a(fragatm(i))%name=="Fe") write(14,*) "0 5" !3d6 4s2 1307 if (a(fragatm(i))%name=="Co") write(14,*) "0 4" !3d7 4s2 1308 if (a(fragatm(i))%name=="Ni") write(14,*) "0 3" !3d8 4s2 1309 if (a(fragatm(i))%name=="Cu") write(14,*) "0 2" !3d10 4s1, needn't sphericalization 1310 if (a(fragatm(i))%name=="Zn") write(14,*) "0 1" !3d10 4s2, needn't sphericalization 1311 write(14,*) a(fragatm(i))%name,0.0,0.0,0.0 1312 write(14,*) 1313 write(14,*) trim(tmpdir)//a(fragatm(i))%name//".wfn" !The output path of wfn file 1314 write(14,*) 1315 write(14,*) 1316 close(14) 1317 end if 1318end do 1319 1320if (noatmwfn==0) then 1321 if (isys==1) call system("del "//trim(tmpdir)//"*.gjf /Q") !The .gjf generated have valueless now, delete them for avoiding user's misunderstanding 1322 if (isys==2.or.isys==3) call system("rm "//trim(tmpdir)//"*.gjf -f") 1323else if (noatmwfn==1) then !Some wfn needs to be genereated by Gaussian and sphericalized here 1324 do i=1,nfragatmnum 1325 nametmp=a_org(fragatm(i))%name 1326 inquire(file=trim(tmpdir)//nametmp//".wfn",exist=alive) 1327 if (alive) cycle !If the .wfn file had copied from atomwfn folder, needn't recalculate 1328 1329 write(*,*) "Running:" 1330 write(*,*) trim(Gaupath)//' "'//trim(tmpdir)//nametmp//'.gjf" "'//trim(tmpdir)//nametmp//'"' 1331 call system(trim(Gaupath)//' "'//trim(tmpdir)//nametmp//'.gjf" "'//trim(tmpdir)//nametmp//'"') 1332 !Check if Gaussian task was successfully finished 1333 if (isys==1) inquire(file=trim(tmpdir)//trim(nametmp)//".out",exist=alivegauout) 1334 if (isys==2.or.isys==3) inquire(file=trim(tmpdir)//trim(nametmp)//".log",exist=alivegauout) 1335 if (alivegauout) then 1336 if (isys==1) open(10,file=trim(tmpdir)//trim(nametmp)//".out",status="old") 1337 if (isys==2.or.isys==3) open(10,file=trim(tmpdir)//trim(nametmp)//".log",status="old") 1338 call loclabel(10,"Normal termination",igaunormal) 1339 close(10) 1340 if (igaunormal==0) then 1341 write(*,"(a)") " Gaussian running may be failed! Please manually check Gaussian input and output files in wfntmp folder. Press ENTER to continue" 1342 read(*,*) 1343 else if (igaunormal==1) then 1344 write(*,*) "Finished successfully!" 1345 end if 1346 else 1347 write(*,"(a)") " Gaussian running may be failed! Please manually check Gaussian input and output files in wfntmp folder" 1348 read(*,*) 1349 end if 1350 1351 !Load and sphericalize electron density for the just generated wfn, and then save 1352 if (ispheratm==1.and.igaunormal==1) then 1353 call dealloall 1354 call readwfn(trim(tmpdir)//nametmp//".wfn",1) 1355 !Main group, restrict open-shell 1356 if (nametmp=="H ".or.nametmp=="Li".or.nametmp=="Na".or.nametmp=="K") MOocc(nmo)=1.0D0 1357 if (nametmp=="B ".or.nametmp=="Al".or.nametmp=="Ga") then 1358 nmo=nmo+2 1359 call resizebynmo(nmo,nprims) !Enlarge nmo by 2, but don't interfere nprims 1360 MOene(nmo-1:nmo)=MOene(nmo-2) 1361 MOtype(nmo-2:nmo)=1 !actually no use, because we only use atomic wfn. files to get total density 1362 MOocc(nmo-2:nmo)=1D0/3D0 1363 call orbcoeffrotate(nmo-2) !XYZ->ZXY, note: nmo-2 is original single occupied orbital 1364 CO(nmo-1,:)=CO(nmo-2,:) 1365 call orbcoeffrotate(nmo-2) !ZXY->YZX 1366 CO(nmo,:)=CO(nmo-2,:) 1367 call orbcoeffrotate(nmo-2) !YZX->XYZ, namely recovered 1368 !Now nmo-2,nmo-1,nmo correspond XYZ,ZXY,YZX 1369 end if 1370 if (nametmp=="C ".or.nametmp=="Si".or.nametmp=="Ge") then 1371 if (SpherIVgroup==0) then 1372 MOocc(nmo-3:nmo)=1.0D0 1373 else if (SpherIVgroup==1) then 1374 nmo=nmo+1 1375 call resizebynmo(nmo,nprims) 1376 MOene(nmo)=MOene(nmo-2) !MOene(nmo-1) is degenerate to MOene(nmo-2) 1377 MOtype(nmo-2:nmo)=1 1378 MOocc(nmo-2:nmo)=2D0/3D0 1379 !Rotate and copy the first occupied p orbital (nmo-5) 1380 call orbcoeffrotate(nmo-2) !XYZ->ZXY 1381 CO(nmo-1,:)=CO(nmo-2,:) !Overlap the already occupied orbital 1382 call orbcoeffrotate(nmo-2) !ZXY->YZX 1383 CO(nmo,:)=CO(nmo-2,:) 1384 call orbcoeffrotate(nmo-2) !YZX->XYZ, namely recovered 1385 end if 1386 end if 1387 if (nametmp=="N ".or.nametmp=="P ".or.nametmp=="As") MOocc(nmo-2:nmo)=1.0D0 1388 if (nametmp=="O ".or.nametmp=="S ".or.nametmp=="Se") MOocc(nmo-2:nmo)=4D0/3D0 1389 if (nametmp=="F ".or.nametmp=="Cl".or.nametmp=="Br") MOocc(nmo-2:nmo)=5D0/3D0 1390 !Transition metals, unrestrict open-shell, find boundary of alpha and beta first 1391 do ibound=2,nmo 1392 if (MOene(ibound)<MOene(ibound-1)) exit !from ii is beta orbitals 1393 end do 1394 !For Sc, Ti and V, rotate and duplicate d orbitals in each diection to get *near* spherical density, as for III main group 1395 !Note: Don't use Hartree-Fock, because correct energy sequence couldn't be reproduced, so can't be sphericalized correctly! 1396 if (nametmp=="Sc".or.nametmp=="Ti".or.nametmp=="V ") then !3d1 4s2, 3d2 4s2, 3d3 4s2 1397 if (nametmp=="Sc") then 1398 ibeg=1 !alpha 4s orbital, because this s orbital shows very strong unequlitity 1399 iend=2 1400 else if (nametmp=="Ti") then 1401 ibeg=2 1402 iend=3 1403 else if (nametmp=="V") then 1404 ibeg=2 1405 iend=4 1406 end if 1407 ienlarge=(iend-ibeg+1)*2 1408 call resizebynmo(nmo+ienlarge,nprims) 1409 ipass=0 1410 do iavgorb=ibeg,iend 1411 call orbcoeffrotate(ibound-iavgorb) !rotate this orbital 1412 CO(nmo+1+ipass,:)=CO(ibound-iavgorb,:) !Duplicate this orbital 1413 call orbcoeffrotate(ibound-iavgorb) 1414 CO(nmo+2+ipass,:)=CO(ibound-iavgorb,:) 1415 call orbcoeffrotate(ibound-iavgorb) !recover 1416 MOocc(ibound-iavgorb)=1D0/3D0 1417 MOene(nmo+1+ipass:nmo+2+ipass)=MOene(ibound-iavgorb) 1418 ipass=ipass+2 !next time skip nmo+1 and nmo+2 1419 end do 1420 MOocc(nmo+1:nmo+ienlarge)=1D0/3D0 1421 MOtype(nmo+1:nmo+ienlarge)=1 1422 nmo=nmo+ienlarge 1423 else if (nametmp=="Fe") then !3d6 4s2 1424 MOocc(nmo-1)=0D0 !delete the only d-beta orbital, the "nmo"th orbital is 4s-beta 1425 MOocc(ibound-6:ibound-2)=1.2D0 !Scatter one electron in beta-d orbital to alpha orbitals evenly. MOocc(ibound) is 4s orbital 1426 else if (nametmp=="Co") then !3d7 4s2 1427 MOocc(nmo-2:nmo-1)=0D0 1428 MOocc(ibound-6:ibound-2)=1.4D0 1429 else if (nametmp=="Ni") then !3d8 4s2 1430 MOocc(nmo-3:nmo-1)=0D0 1431 MOocc(ibound-6:ibound-2)=1.6D0 1432 end if 1433 call outwfn(trim(tmpdir)//nametmp//".wfn",0,0,10) 1434 end if 1435 end do 1436end if 1437write(*,*) 1438 1439!Setup custom operation list 1440ncustommap=nfragatmnum 1441if (allocated(custommapname)) deallocate(custommapname) 1442if (allocated(customop)) deallocate(customop) 1443allocate(custommapname(ncustommap)) 1444allocate(customop(ncustommap)) 1445customop='-' 1446if (ipromol==1) customop='+' !Do promolecular map 1447 1448!Generate atomic wfn file from element wfn file, meanwhile take them into custom operation list 1449do i=1,itype !Scan each atomtype in current system 1450 call dealloall 1451 call readwfn(trim(tmpdir)//typename(i)//".wfn",1) 1452 do j=1,nfragatmnum 1453 if (a_org(fragatm(j))%name==typename(i)) then !Find atoms attributed to current element 1454 a(1)%x=a_org(fragatm(j))%x !Modify the atomic .wfn, then output to new .wfn 1455 a(1)%y=a_org(fragatm(j))%y 1456 a(1)%z=a_org(fragatm(j))%z 1457 write(outwfnname,"(a2,i4,a4)") typename(i),fragatm(j),".wfn" 1458 call outwfn(trim(tmpdir)//outwfnname,0,0,10) 1459 custommapname(j)=trim(tmpdir)//outwfnname !Sequence is identical to atom in fragment 1460 end if 1461 end do 1462end do 1463 1464call dealloall 1465call readinfile(firstfilename,1) 1466! nmo=nmo_org !Original .wfn may be maniplated by user, retrieve them 1467! b=b_org 1468! CO=CO_org 1469! nprims=nprims_org 1470end subroutine 1471 1472 1473 1474!!!------------------ Generate density matrix, can be used when basis function information is available 1475subroutine genP 1476use defvar 1477implicit real*8 (a-h,o-z) 1478if (allocated(Ptot)) deallocate(Ptot) 1479if (allocated(Palpha)) deallocate(Palpha) 1480if (allocated(Pbeta)) deallocate(Pbeta) 1481allocate(Ptot(nbasis,nbasis)) 1482Ptot=0 1483if (wfntype==1.or.wfntype==2.or.wfntype==4) then !open-shell 1484 allocate(Palpha(nbasis,nbasis)) 1485 allocate(Pbeta(nbasis,nbasis)) 1486 Palpha=0D0 1487 Pbeta=0D0 1488end if 1489 1490!For SCF wavefunction, if the wavefunction has not been modified (imodwfn==0), use fast way to construct it 1491!However, if the wavefunction has been modified, the case may be complicated, for example, there is a hole orbital. In these cases 1492!We use general way (as used for post-HF) to construct density matrix 1493 1494if (wfntype==0.and.imodwfn==0) then !RHF 1495 Ptot=2*matmul(CObasa(:,1:nint(naelec)),transpose(CObasa(:,1:nint(naelec)))) 1496else if (wfntype==1.and.imodwfn==0) then !UHF 1497 Palpha=matmul(CObasa(:,1:nint(naelec)),transpose(CObasa(:,1:nint(naelec)))) 1498 Pbeta=matmul(CObasb(:,1:nint(nbelec)),transpose(CObasb(:,1:nint(nbelec)))) 1499 Ptot=Palpha+Pbeta 1500else if (wfntype==2.and.imodwfn==0) then !ROHF 1501 Palpha=matmul(CObasa(:,1:nint(naelec)),transpose(CObasa(:,1:nint(naelec)))) 1502 Pbeta=matmul(CObasa(:,1:nint(nbelec)),transpose(CObasa(:,1:nint(nbelec)))) 1503 Ptot=Palpha+Pbeta 1504else if (wfntype==3.or.((wfntype==0.or.wfntype==2).and.imodwfn==1)) then !Restricted post-HF 1505 do imo=1,nmo 1506 if (MOocc(imo)==0D0) cycle 1507 Ptot=Ptot+MOocc(imo)*matmul(CObasa(:,imo:imo),transpose(CObasa(:,imo:imo))) 1508 end do 1509else if (wfntype==4.or.(wfntype==1.and.imodwfn==1)) then !Unrestricted post-HF 1510 do imo=1,nbasis 1511 if (MOocc(imo)==0D0) cycle 1512 Palpha=Palpha+MOocc(imo)*matmul(CObasa(:,imo:imo),transpose(CObasa(:,imo:imo))) 1513 end do 1514 do imo=1,nbasis 1515 if (MOocc(imo+nbasis)==0D0) cycle 1516 Pbeta=Pbeta+MOocc(imo+nbasis)*matmul(CObasb(:,imo:imo),transpose(CObasb(:,imo:imo))) 1517 end do 1518 Ptot=Palpha+Pbeta 1519end if 1520 1521end subroutine 1522 1523 1524 1525!!!------------------ Evaluate overlap integral for two unnormalized GTFs, a warpper of doSintactual for simplicity 1526real*8 function doSint(iGTF,jGTF) 1527integer iGTF,jGTF 1528real*8,external :: doSintactual 1529doSint=doSintactual(iGTF,jGTF,0,0,0,0,0,0) 1530end function 1531!!!------------------ Evaluate overlap integral for two unnormalized GTFs 1532!~p arguments are the shifts of GTF index, used by doKint, doveloint but not by doSint 1533real*8 function doSintactual(iGTF,jGTF,ix1p,iy1p,iz1p,ix2p,iy2p,iz2p) 1534use util 1535use defvar 1536implicit real*8(a-h,o-z) 1537integer iGTF,jGTF,ix1p,iy1p,iz1p,ix2p,iy2p,iz2p 1538x1=a(b(iGTF)%center)%x 1539y1=a(b(iGTF)%center)%y 1540z1=a(b(iGTF)%center)%z 1541x2=a(b(jGTF)%center)%x 1542y2=a(b(jGTF)%center)%y 1543z2=a(b(jGTF)%center)%z 1544ee1=b(iGTF)%exp 1545ee2=b(jGTF)%exp 1546ep=ee1+ee2 1547sqrtep=dsqrt(ep) 1548px=(ee1*x1+ee2*x2)/ep 1549py=(ee1*y1+ee2*y2)/ep 1550pz=(ee1*z1+ee2*z2)/ep 1551expterm=dexp( -ee1*ee2*((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)/ep ) 1552ix1=type2ix(b(iGTF)%functype)+ix1p 1553iy1=type2iy(b(iGTF)%functype)+iy1p 1554iz1=type2iz(b(iGTF)%functype)+iz1p 1555ix2=type2ix(b(jGTF)%functype)+ix2p 1556iy2=type2iy(b(jGTF)%functype)+iy2p 1557iz2=type2iz(b(jGTF)%functype)+iz2p 1558!chen book,P103 1559numx=ceiling( (ix1+ix2+1)/2D0 ) !Need to calculate n points 1560sx=0.0D0 1561do i=1,numx 1562 tmp=Rhm(numx,i)/sqrtep+px 1563 term1=(tmp-x1)**ix1 1564 term2=(tmp-x2)**ix2 1565 sx=sx+Whm(numx,i)*term1*term2 1566end do 1567sx=sx/sqrtep 1568 1569numy=ceiling( (iy1+iy2+1)/2D0 ) 1570sy=0.0D0 1571do i=1,numy 1572 tmp=Rhm(numy,i)/sqrtep+py 1573 term1=(tmp-y1)**iy1 1574 term2=(tmp-y2)**iy2 1575 sy=sy+Whm(numy,i)*term1*term2 1576end do 1577sy=sy/sqrtep 1578 1579numz=ceiling( (iz1+iz2+1)/2D0 ) 1580sz=0.0D0 1581do i=1,numz 1582 tmp=Rhm(numz,i)/sqrtep+pz 1583 term1=(tmp-z1)**iz1 1584 term2=(tmp-z2)**iz2 1585 sz=sz+Whm(numz,i)*term1*term2 1586end do 1587sz=sz/sqrtep 1588 1589doSintactual=sx*sy*sz*expterm 1590end function 1591!!!------------------ Generate overlap matrix between all GTFs 1592!nsize should be nprims*(nprims+1)/2 1593subroutine genGTFSmat(GTFSmat,nsize) 1594use defvar 1595implicit real*8 (a-h,o-z) 1596integer nsize 1597real*8 GTFSmat(nsize) 1598nthreads=getNThreads() 1599!$OMP PARALLEL DO SHARED(GTFSmat) PRIVATE(ides,iGTF,jGTF) schedule(dynamic) NUM_THREADS(nthreads) 1600do iGTF=1,nprims 1601 do jGTF=iGTF,nprims 1602 ides=jGTF*(jGTF-1)/2+iGTF 1603 GTFSmat(ides)=doSint(iGTF,jGTF) 1604 end do 1605end do 1606!$OMP END PARALLEL DO 1607end subroutine 1608!!!------------------ Generate overlap matrix between all basis functions 1609!Sbas should be allocated first. The resultant matrix is for Cartesian basis functions, may be converted to spherical-harmonic later 1610subroutine genSbas 1611use defvar 1612implicit real*8 (a-h,o-z) 1613Sbas=0D0 1614nthreads=getNThreads() 1615!$OMP PARALLEL DO SHARED(Sbas) PRIVATE(i,ii,j,jj) schedule(dynamic) NUM_THREADS(nthreads) 1616do i=1,nbasis 1617 do j=i,nbasis 1618 do ii=primstart(i),primend(i) 1619 do jj=primstart(j),primend(j) 1620 Sbas(i,j)=Sbas(i,j)+primconnorm(ii)*primconnorm(jj)*doSint(ii,jj) 1621 end do 1622 end do 1623 Sbas(j,i)=Sbas(i,j) 1624 end do 1625end do 1626!$OMP END PARALLEL DO 1627end subroutine 1628 1629 1630 1631!!!-------- Evaluate dipole moment integral for two unnormalized GTFs. The negative charge of electron has been considered! 1632!~p arguments are the shifts of GTF index as doSintactual 1633!xint/yint/zint correspond to dipole moment integral in X/Y/Z 1634subroutine dodipoleint(iGTF,jGTF,ix1p,iy1p,iz1p,ix2p,iy2p,iz2p,xint,yint,zint) 1635use util 1636use defvar 1637implicit real*8(a-h,o-z) 1638real*8 xint,yint,zint 1639integer iGTF,jGTF,ix1p,iy1p,iz1p,ix2p,iy2p,iz2p 1640x1=a(b(iGTF)%center)%x 1641y1=a(b(iGTF)%center)%y 1642z1=a(b(iGTF)%center)%z 1643x2=a(b(jGTF)%center)%x 1644y2=a(b(jGTF)%center)%y 1645z2=a(b(jGTF)%center)%z 1646ee1=b(iGTF)%exp 1647ee2=b(jGTF)%exp 1648ep=ee1+ee2 1649sqrtep=dsqrt(ep) 1650px=(ee1*x1+ee2*x2)/ep 1651py=(ee1*y1+ee2*y2)/ep 1652pz=(ee1*z1+ee2*z2)/ep 1653expterm=dexp( -ee1*ee2*((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)/ep ) 1654ix1=type2ix(b(iGTF)%functype)+ix1p 1655iy1=type2iy(b(iGTF)%functype)+iy1p 1656iz1=type2iz(b(iGTF)%functype)+iz1p 1657ix2=type2ix(b(jGTF)%functype)+ix2p 1658iy2=type2iy(b(jGTF)%functype)+iy2p 1659iz2=type2iz(b(jGTF)%functype)+iz2p 1660!First, calculate sx,sy,sz as usual as doSint 1661numx=ceiling( (ix1+ix2+1)/2D0 ) !Need to calculate n points 1662sx=0.0D0 1663do i=1,numx 1664 tmp=Rhm(numx,i)/sqrtep+px 1665 term1=(tmp-x1)**ix1 1666 term2=(tmp-x2)**ix2 1667 sx=sx+Whm(numx,i)*term1*term2 1668end do 1669sx=sx/sqrtep 1670numy=ceiling( (iy1+iy2+1)/2D0 ) 1671sy=0.0D0 1672do i=1,numy 1673 tmp=Rhm(numy,i)/sqrtep+py 1674 term1=(tmp-y1)**iy1 1675 term2=(tmp-y2)**iy2 1676 sy=sy+Whm(numy,i)*term1*term2 1677end do 1678sy=sy/sqrtep 1679numz=ceiling( (iz1+iz2+1)/2D0 ) 1680sz=0.0D0 1681do i=1,numz 1682 tmp=Rhm(numz,i)/sqrtep+pz 1683 term1=(tmp-z1)**iz1 1684 term2=(tmp-z2)**iz2 1685 sz=sz+Whm(numz,i)*term1*term2 1686end do 1687sz=sz/sqrtep 1688!Second, calculate overlap integral in X,Y,Z directions but with X,Y,Z coordinate variables (relative to the original point of the whole system) to produce sxx,syy,szz 1689numx=ceiling( (ix1+ix2+2)/2D0 ) !Because X variable is introduced, ix1+ix2+2 is used instead of ix1+ix2+1 1690sxx=0.0D0 1691do i=1,numx 1692 tmp=Rhm(numx,i)/sqrtep+px 1693 term1=(tmp-x1)**ix1 1694 term2=(tmp-x2)**ix2 1695 sxx=sxx+Whm(numx,i)*term1*term2*tmp 1696end do 1697sxx=sxx/sqrtep 1698numy=ceiling( (iy1+iy2+2)/2D0 ) 1699syy=0.0D0 1700do i=1,numy 1701 tmp=Rhm(numy,i)/sqrtep+py 1702 term1=(tmp-y1)**iy1 1703 term2=(tmp-y2)**iy2 1704 syy=syy+Whm(numy,i)*term1*term2*tmp 1705end do 1706syy=syy/sqrtep 1707numz=ceiling( (iz1+iz2+2)/2D0 ) 1708szz=0.0D0 1709do i=1,numz 1710 tmp=Rhm(numz,i)/sqrtep+pz 1711 term1=(tmp-z1)**iz1 1712 term2=(tmp-z2)**iz2 1713 szz=szz+Whm(numz,i)*term1*term2*tmp 1714end do 1715szz=szz/sqrtep 1716 1717xint=-sxx*sy*sz*expterm 1718yint=-sx*syy*sz*expterm 1719zint=-sx*sy*szz*expterm 1720end subroutine 1721!------ A warpper of dodipoleint, used to directly get a single component of dipole moment integral. icomp=1/2/3 corresponds to X/Y/Z component 1722real*8 function dipintcomp(icomp,iGTF,jGTF,ix1p,iy1p,iz1p,ix2p,iy2p,iz2p) 1723integer icomp,iGTF,jGTF,ix1p,iy1p,iz1p,ix2p,iy2p,iz2p 1724real*8 xcomp,ycomp,zcomp 1725call dodipoleint(iGTF,jGTF,ix1p,iy1p,iz1p,ix2p,iy2p,iz2p,xcomp,ycomp,zcomp) 1726if (icomp==1) dipintcomp=xcomp 1727if (icomp==2) dipintcomp=ycomp 1728if (icomp==3) dipintcomp=zcomp 1729end function 1730!!!--------------- Generate dipole moment integral matrix between all GTFs 1731!nsize should be nprims*(nprims+1)/2 1732subroutine genGTFDmat(GTFdipmat,nsize) 1733use defvar 1734implicit real*8 (a-h,o-z) 1735integer nsize 1736real*8 GTFdipmat(3,nsize) 1737nthreads=getNThreads() 1738!$OMP PARALLEL DO SHARED(GTFdipmat) PRIVATE(ides,iGTF,jGTF,xdiptmp,ydiptmp,zdiptmp) schedule(dynamic) NUM_THREADS(nthreads) 1739do iGTF=1,nprims 1740 do jGTF=iGTF,nprims 1741 ides=jGTF*(jGTF-1)/2+iGTF 1742 call dodipoleint(iGTF,jGTF,0,0,0,0,0,0,xdiptmp,ydiptmp,zdiptmp) 1743 GTFdipmat(1,ides)=xdiptmp 1744 GTFdipmat(2,ides)=ydiptmp 1745 GTFdipmat(3,ides)=zdiptmp 1746 end do 1747end do 1748!$OMP END PARALLEL DO 1749end subroutine 1750!!!------------------ Generate dipole moment integral matrix between all basis functions 1751!Dbas should be allocated first. The resultant matrix is for Cartesian basis functions, may be converted to spherical-harmonic later 1752subroutine genDbas 1753use defvar 1754implicit real*8 (a-h,o-z) 1755Dbas=0D0 1756nthreads=getNThreads() 1757!$OMP PARALLEL DO SHARED(Dbas) PRIVATE(i,ii,j,jj,xdiptmp,ydiptmp,zdiptmp) schedule(dynamic) NUM_THREADS(nthreads) 1758do i=1,nbasis 1759 do j=i,nbasis 1760 do ii=primstart(i),primend(i) 1761 do jj=primstart(j),primend(j) 1762 call dodipoleint(ii,jj,0,0,0,0,0,0,xdiptmp,ydiptmp,zdiptmp) 1763 Dbas(1,i,j)=Dbas(1,i,j)+primconnorm(ii)*primconnorm(jj)*xdiptmp 1764 Dbas(2,i,j)=Dbas(2,i,j)+primconnorm(ii)*primconnorm(jj)*ydiptmp 1765 Dbas(3,i,j)=Dbas(3,i,j)+primconnorm(ii)*primconnorm(jj)*zdiptmp 1766 end do 1767 end do 1768 Dbas(:,j,i)=Dbas(:,i,j) 1769 end do 1770end do 1771!$OMP END PARALLEL DO 1772end subroutine 1773 1774 1775 1776!!!------- Evaluate magnetic integral for two unnormalized GTFs 1777!The imaginary sign i is ignored. Note that the negative sign of magnetic operator is not occurred here 1778!Consult doVelint for the method for evaluation of <a|d/dx|b>, and TCA,6,341 for the formula of magnetic integral 1779subroutine doMagint(iGTF,jGTF,xcomp,ycomp,zcomp) 1780use defvar 1781implicit real*8(a-h,o-z) 1782integer iGTF,jGTF 1783real*8 xcomp,ycomp,zcomp,term(4) 1784ee1=b(iGTF)%exp 1785ee2=b(jGTF)%exp 1786ix1=type2ix(b(iGTF)%functype) 1787iy1=type2iy(b(iGTF)%functype) 1788iz1=type2iz(b(iGTF)%functype) 1789ix2=type2ix(b(jGTF)%functype) 1790iy2=type2iy(b(jGTF)%functype) 1791iz2=type2iz(b(jGTF)%functype) 1792term=0 1793!X component, <a|y*d/dz-z*d/dy|b>. Since <a|d/dz|b>=iz2*<a|b-1z>-2*ee2*<a|b+1z>, the term such as <a|y*d/dz|b> can be evaluated in terms of dipole integrals iz2*<a|y|b-1z>-2*ee2*<a|y|b+1z> 1794if(iz2>0) term(1)= iz2*dipintcomp(2,iGTF,jGTF,0,0,0,0,0,-1) !viz. iz2*<a|y|b-1z> 1795 term(2)=-2*ee2*dipintcomp(2,iGTF,jGTF,0,0,0,0,0, 1) 1796if(iy2>0) term(3)= -iy2*dipintcomp(3,iGTF,jGTF,0,0,0,0,-1,0) 1797 term(4)= 2*ee2*dipintcomp(3,iGTF,jGTF,0,0,0,0, 1,0) 1798xcomp=-sum(term) !Note that the result of dipintcomp has a negative sign due to the negative charge of electron, so here revise the sign 1799term=0 1800!Y component, <a|z*d/dx-x*d/dz|b> 1801if(ix2>0) term(1)= ix2*dipintcomp(3,iGTF,jGTF,0,0,0,-1,0,0) 1802 term(2)=-2*ee2*dipintcomp(3,iGTF,jGTF,0,0,0, 1,0,0) 1803if(iz2>0) term(3)= -iz2*dipintcomp(1,iGTF,jGTF,0,0,0,0,0,-1) 1804 term(4)= 2*ee2*dipintcomp(1,iGTF,jGTF,0,0,0,0,0, 1) 1805ycomp=-sum(term) 1806term=0 1807!Z component, <a|x*d/dy-y*d/dx|b> 1808if(iy2>0) term(1)= iy2*dipintcomp(1,iGTF,jGTF,0,0,0,0,-1,0) 1809 term(2)=-2*ee2*dipintcomp(1,iGTF,jGTF,0,0,0,0, 1,0) 1810if(ix2>0) term(3)= -ix2*dipintcomp(2,iGTF,jGTF,0,0,0,-1,0,0) 1811 term(4)= 2*ee2*dipintcomp(2,iGTF,jGTF,0,0,0, 1,0,0) 1812zcomp=-sum(term) 1813end subroutine 1814!!!--------------- Generate magnetic dipole moment integral matrix between all GTFs 1815!nsize should be nprims*(nprims+1)/2 1816!Beware that when using this result, (j,i) element should be set to negative value of (i,j) due to the Hermitean character of this operator! 1817subroutine genGTFMmat(GTFdipmat,nsize) 1818use defvar 1819implicit real*8 (a-h,o-z) 1820integer nsize 1821real*8 GTFdipmat(3,nsize) 1822GTFdipmat=0D0 1823nthreads=getNThreads() 1824!$OMP PARALLEL DO SHARED(GTFdipmat) PRIVATE(ides,iGTF,jGTF,xdiptmp,ydiptmp,zdiptmp) schedule(dynamic) NUM_THREADS(nthreads) 1825do iGTF=1,nprims 1826 do jGTF=iGTF+1,nprims !For iGTF=jGTF, the value must exactly zero, so don't calculate 1827 ides=jGTF*(jGTF-1)/2+iGTF 1828 call doMagint(iGTF,jGTF,xdiptmp,ydiptmp,zdiptmp) 1829 GTFdipmat(1,ides)=xdiptmp 1830 GTFdipmat(2,ides)=ydiptmp 1831 GTFdipmat(3,ides)=zdiptmp 1832 end do 1833end do 1834!$OMP END PARALLEL DO 1835end subroutine 1836!!!------------------ Generate magnetic integral matrix between all basis functions 1837!Magbas should be allocated first. The resultant matrix is for Cartesian basis functions, may be converted to spherical-harmonic later 1838!Notice that the diagonal element of magnetic integral matrix is zero, and (i,j)=-(j,i) due to Hermitean character 1839subroutine genMagbas 1840use defvar 1841implicit real*8 (a-h,o-z) 1842Magbas=0D0 1843nthreads=getNThreads() 1844!$OMP PARALLEL DO SHARED(Magbas) PRIVATE(i,ii,j,jj,xcomp,ycomp,zcomp) schedule(dynamic) NUM_THREADS(nthreads) 1845do i=1,nbasis 1846 do j=i+1,nbasis 1847 do ii=primstart(i),primend(i) 1848 do jj=primstart(j),primend(j) 1849 call doMagint(ii,jj,xcomp,ycomp,zcomp) 1850 Magbas(1,i,j)=Magbas(1,i,j)+primconnorm(ii)*primconnorm(jj)*xcomp 1851 Magbas(2,i,j)=Magbas(2,i,j)+primconnorm(ii)*primconnorm(jj)*ycomp 1852 Magbas(3,i,j)=Magbas(3,i,j)+primconnorm(ii)*primconnorm(jj)*zcomp 1853 end do 1854 end do 1855 Magbas(:,j,i)=-Magbas(:,i,j) 1856 end do 1857end do 1858!$OMP END PARALLEL DO 1859end subroutine 1860 1861 1862 1863!!!------- Evaluate velocity integral for two unnormalized GTFs 1864!There are three components. e.g. X direction: i<a|d/dx|b>. The imaginary sign i is ignored. Note that the negative sign of momentum operator is not occurred here 1865!One can consult p97 of Chen's book for the derivative of GTF. Namely <a|d/dx|b>=ix2*<a|b-1x>-2*ee2*<a|b+1x> 1866subroutine doVelint(iGTF,jGTF,xcomp,ycomp,zcomp) 1867use defvar 1868implicit real*8(a-h,o-z) 1869integer iGTF,jGTF 1870real*8 xcomp,ycomp,zcomp 1871ee1=b(iGTF)%exp 1872ee2=b(jGTF)%exp 1873ix1=type2ix(b(iGTF)%functype) 1874iy1=type2iy(b(iGTF)%functype) 1875iz1=type2iz(b(iGTF)%functype) 1876ix2=type2ix(b(jGTF)%functype) 1877iy2=type2iy(b(jGTF)%functype) 1878iz2=type2iz(b(jGTF)%functype) 1879term1=0 1880if(ix2>0) term1= ix2*doSintactual(iGTF,jGTF,0,0,0,-1,0,0) 1881 term2=-2*ee2*doSintactual(iGTF,jGTF,0,0,0, 1,0,0) 1882xcomp=term1+term2 1883term1=0 1884if(iy2>0) term1= iy2*doSintactual(iGTF,jGTF,0,0,0,0,-1,0) 1885 term2=-2*ee2*doSintactual(iGTF,jGTF,0,0,0,0, 1,0) 1886ycomp=term1+term2 1887term1=0 1888if(iz2>0) term1= iz2*doSintactual(iGTF,jGTF,0,0,0,0,0,-1) 1889 term2=-2*ee2*doSintactual(iGTF,jGTF,0,0,0,0,0, 1) 1890zcomp=term1+term2 1891end subroutine 1892!!!--------------- Generate velocity integral matrix between all GTFs 1893!nsize should be nprims*(nprims+1)/2 1894!Beware that when using this result, (j,i) element should be set to negative value of (i,j) due to the Hermitean character of this operator! 1895subroutine genGTFVelmat(GTFVelmat,nsize) 1896use defvar 1897implicit real*8 (a-h,o-z) 1898integer nsize 1899real*8 GTFVelmat(3,nsize) 1900nthreads=getNThreads() 1901!$OMP PARALLEL DO SHARED(GTFVelmat) PRIVATE(ides,iGTF,jGTF,xtmp,ytmp,ztmp) schedule(dynamic) NUM_THREADS(nthreads) 1902do iGTF=1,nprims 1903 do jGTF=iGTF,nprims 1904 ides=jGTF*(jGTF-1)/2+iGTF 1905 call doVelint(iGTF,jGTF,xtmp,ytmp,ztmp) 1906 GTFVelmat(1,ides)=xtmp 1907 GTFVelmat(2,ides)=ytmp 1908 GTFVelmat(3,ides)=ztmp 1909 end do 1910end do 1911!$OMP END PARALLEL DO 1912end subroutine 1913!!!------------------ Generate velocity integral matrix between all basis functions 1914!Velbas should be allocated first. The resultant matrix is for Cartesian basis functions, may be converted to spherical-harmonic later 1915!Notice that the diagonal element of velocity integral matrix is zero, and (i,j)=-(j,i) due to Hermitean character 1916subroutine genVelbas 1917use defvar 1918implicit real*8 (a-h,o-z) 1919Velbas=0D0 1920nthreads=getNThreads() 1921!$OMP PARALLEL DO SHARED(Velbas) PRIVATE(i,ii,j,jj,xcomp,ycomp,zcomp) schedule(dynamic) NUM_THREADS(nthreads) 1922do i=1,nbasis 1923 do j=i+1,nbasis 1924 do ii=primstart(i),primend(i) 1925 do jj=primstart(j),primend(j) 1926 call doVelint(ii,jj,xcomp,ycomp,zcomp) 1927 Velbas(1,i,j)=Velbas(1,i,j)+primconnorm(ii)*primconnorm(jj)*xcomp 1928 Velbas(2,i,j)=Velbas(2,i,j)+primconnorm(ii)*primconnorm(jj)*ycomp 1929 Velbas(3,i,j)=Velbas(3,i,j)+primconnorm(ii)*primconnorm(jj)*zcomp 1930 end do 1931 end do 1932 Velbas(:,j,i)=-Velbas(:,i,j) 1933 end do 1934end do 1935!$OMP END PARALLEL DO 1936end subroutine 1937 1938 1939 1940!!!------------------- Evaluate kinetic integral (i.e. -(1/2)der2 )for two unnormalized GTFs, see Chen's book, p104 1941real*8 function doTint(iGTF,jGTF) 1942use defvar 1943implicit real*8(a-h,o-z) 1944integer iGTF,jGTF 1945real*8 term(4) 1946ee1=b(iGTF)%exp 1947ee2=b(jGTF)%exp 1948ix1=type2ix(b(iGTF)%functype) 1949iy1=type2iy(b(iGTF)%functype) 1950iz1=type2iz(b(iGTF)%functype) 1951ix2=type2ix(b(jGTF)%functype) 1952iy2=type2iy(b(jGTF)%functype) 1953iz2=type2iz(b(jGTF)%functype) 1954term=0 1955if(ix1>0.and.ix2>0) term(1)= ix1*ix2*doSintactual(iGTF,jGTF,-1,0,0,-1,0,0) 1956if(ix1>0) term(2)=-2*ee2*ix1*doSintactual(iGTF,jGTF,-1,0,0, 1,0,0) 1957if(ix2>0) term(3)=-2*ee1*ix2*doSintactual(iGTF,jGTF, 1,0,0,-1,0,0) 1958 term(4)= 4*ee1*ee2*doSintactual(iGTF,jGTF, 1,0,0, 1,0,0) 1959Tx=sum(term) 1960term=0 1961if(iy1>0.and.iy2>0) term(1)= iy1*iy2*doSintactual(iGTF,jGTF,0,-1,0,0,-1,0) 1962if(iy1>0) term(2)=-2*ee2*iy1*doSintactual(iGTF,jGTF,0,-1,0,0, 1,0) 1963if(iy2>0) term(3)=-2*ee1*iy2*doSintactual(iGTF,jGTF,0, 1,0,0,-1,0) 1964 term(4)= 4*ee1*ee2*doSintactual(iGTF,jGTF,0, 1,0,0, 1,0) 1965Ty=sum(term) 1966term=0 1967if(iz1>0.and.iz2>0) term(1)= iz1*iz2*doSintactual(iGTF,jGTF,0,0,-1,0,0,-1) 1968if(iz1>0) term(2)=-2*ee2*iz1*doSintactual(iGTF,jGTF,0,0,-1,0,0, 1) 1969if(iz2>0) term(3)=-2*ee1*iz2*doSintactual(iGTF,jGTF,0,0, 1,0,0,-1) 1970 term(4)= 4*ee1*ee2*doSintactual(iGTF,jGTF,0,0, 1,0,0, 1) 1971Tz=sum(term) 1972doTint=(Tx+Ty+Tz)/2 1973end function 1974!!!------------------ Generate kinetic energy matrix between all GTFs 1975!nsize should be nprims*(nprims+1)/2 1976subroutine genGTFTmat(GTFTmat,nsize) 1977use defvar 1978implicit real*8 (a-h,o-z) 1979integer nsize 1980real*8 GTFTmat(nsize) 1981nthreads=getNThreads() 1982!$OMP PARALLEL DO SHARED(GTFTmat) PRIVATE(ides,iGTF,jGTF) schedule(dynamic) NUM_THREADS(nthreads) 1983do iGTF=1,nprims 1984 do jGTF=iGTF,nprims 1985 ides=jGTF*(jGTF-1)/2+iGTF 1986 GTFTmat(ides)=doTint(iGTF,jGTF) 1987 end do 1988end do 1989!$OMP END PARALLEL DO 1990end subroutine 1991!!!------------------ Generate kinetic energy matrix between all basis functions 1992!Tbas should be allocated first. The resultant matrix is for Cartesian basis functions, may be converted to spherical-harmonic later 1993subroutine genTbas 1994use defvar 1995implicit real*8 (a-h,o-z) 1996Tbas=0D0 1997nthreads=getNThreads() 1998!$OMP PARALLEL DO SHARED(Tbas) PRIVATE(i,ii,j,jj) schedule(dynamic) NUM_THREADS(nthreads) 1999do i=1,nbasis 2000 do j=i,nbasis 2001 do ii=primstart(i),primend(i) 2002 do jj=primstart(j),primend(j) 2003 Tbas(i,j)=Tbas(i,j)+primconnorm(ii)*primconnorm(jj)*doTint(ii,jj) 2004 end do 2005 end do 2006 Tbas(j,i)=Tbas(i,j) 2007 end do 2008end do 2009!$OMP END PARALLEL DO 2010end subroutine 2011 2012 2013 2014!!!------ Show system one-electron properties based on density matrix and integral matrix between basis functions 2015!The results are correct only when Cartesian basis functions are used 2016subroutine sys1eprop 2017use defvar 2018if (.not.allocated(Sbas)) allocate(Sbas(nbasis,nbasis)) 2019call genSbas 2020write(*,"(' Total number of electrons:',f16.8)") sum(Ptot*Sbas) 2021if (.not.allocated(Tbas)) allocate(Tbas(nbasis,nbasis)) 2022call genTbas 2023write(*,"(' Kinetic energy:',f18.9,' a.u.')") sum(Ptot*Tbas) 2024if (.not.allocated(Dbas)) allocate(Dbas(3,nbasis,nbasis)) 2025call genDbas 2026write(*,"(' Electric dipole moment in X/Y/Z:',3f13.7,' a.u.')") sum(Ptot*Dbas(1,:,:)),sum(Ptot*Dbas(2,:,:)),sum(Ptot*Dbas(3,:,:)) 2027if (.not.allocated(Magbas)) allocate(Magbas(3,nbasis,nbasis)) 2028call genMagbas 2029write(*,"(' Magnetic dipole moment in X/Y/Z:',3f13.7,' a.u.')") sum(Ptot*Magbas(1,:,:)),sum(Ptot*Magbas(2,:,:)),sum(Ptot*Magbas(3,:,:)) 2030if (.not.allocated(Velbas)) allocate(Velbas(3,nbasis,nbasis)) 2031call genVelbas 2032write(*,"(' Linear momentum in X/Y/Z: ',3f13.7,' a.u.')") sum(Ptot*Velbas(1,:,:)),sum(Ptot*Velbas(2,:,:)),sum(Ptot*Velbas(3,:,:)) 2033end subroutine 2034 2035 2036 2037!!!------------- Show all properties at a point 2038!ifuncsel: controls the gradient and Hessian for which function 2039!ifileid: Output to which file destination, of course 6=screen 2040subroutine showptprop(inx,iny,inz,ifuncsel,ifileid) 2041use util 2042use defvar 2043use function 2044implicit real*8(a-h,o-z) 2045real*8 inx,iny,inz 2046real*8 elehess(3,3),eigvecmat(3,3),eigval(3),elegrad(3),funchess(3,3),funcgrad(3),tmparr(3,1) 2047integer ifuncsel,ifileid 2048if (allocated(b)) then !If loaded file contains wavefuntion information 2049 call gencalchessmat(2,1,inx,iny,inz,elerho,elegrad,elehess) !Generate electron density, gradient and hessian 2050 write(ifileid,"(' Density of all electrons:',E18.10)") elerho 2051 if (ipolarpara==0) then 2052 tmpval=fspindens(inx,iny,inz,'s') 2053 write(ifileid,"(' Density of Alpha electrons:',E18.10)") (elerho+tmpval)/2D0 2054 write(ifileid,"(' Density of Beta electrons:',E18.10)") (elerho-tmpval)/2D0 2055 write(ifileid,"(' Spin density of electrons:',E18.10)") tmpval 2056 else if (ipolarpara==1) then 2057 write(ifileid,"(' Spin polarization parameter function:',E18.10)") fspindens(inx,iny,inz,'s') 2058 end if 2059 valG=lagkin(inx,iny,inz,0) 2060 valGx=lagkin(inx,iny,inz,1) 2061 valGy=lagkin(inx,iny,inz,2) 2062 valGz=lagkin(inx,iny,inz,3) 2063 write(ifileid,"(' Lagrangian kinetic energy G(r):',E18.10)") valG 2064 write(ifileid,"(' G(r) in X,Y,Z:',3E18.10)") valGx,valGy,valGz 2065 valK=Hamkin(inx,iny,inz,0) 2066! valKx=Hamkin(inx,iny,inz,1) 2067! valKy=Hamkin(inx,iny,inz,2) 2068! valKz=Hamkin(inx,iny,inz,3) 2069 write(ifileid,"(' Hamiltonian kinetic energy K(r):',E18.10)") valK 2070! write(ifileid,"(' K(r) in X,Y,Z:',3E18.10)") valKx,valKy,valKz 2071 write(ifileid,"(' Potential energy density V(r):',E18.10)") -valK-valG !When without EDF, also equals to flapl(inx,iny,inz,'t')/4.0D0-2*valG 2072 write(ifileid,"(' Energy density E(r) or H(r):',E18.10)") -valK 2073 write(ifileid,"(' Laplacian of electron density:',E18.10)") laplfac*(elehess(1,1)+elehess(2,2)+elehess(3,3)) 2074 write(ifileid,"(' Electron localization function (ELF):',E18.10)") ELF_LOL(inx,iny,inz,"ELF") 2075 write(ifileid,"(' Localized orbital locator (LOL):',E18.10)") ELF_LOL(inx,iny,inz,"LOL") 2076 write(ifileid,"(' Local information entropy:',E18.10)") infoentro(1,inx,iny,inz) 2077 write(ifileid,"(' Reduced density gradient (RDG):',E18.10)") fgrad(inx,iny,inz,'r') 2078 write(ifileid,"(' Reduced density gradient with promolecular approximation:',E18.10)") RDGprodens(inx,iny,inz) 2079 write(ifileid,"(' Sign(lambda2)*rho:',E18.10)") signlambda2rho(inx,iny,inz) 2080 write(ifileid,"(' Sign(lambda2)*rho with promolecular approximation:',E18.10)") signlambda2rho_prodens(inx,iny,inz) 2081 if (pairfunctype==1) write(ifileid,"(a,3f10.5,' :',E18.10)") " Corr. hole for alpha, ref.:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz) 2082 if (pairfunctype==2) write(ifileid,"(a,3f10.5,' :',E18.10)") " Corr. hole for beta, ref.:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz) 2083 if (pairfunctype==4) write(ifileid,"(a,3f10.5,' :',E18.10)") " Corr. fac. for alpha, ref.:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz) 2084 if (pairfunctype==5) write(ifileid,"(a,3f10.5,' :',E18.10)") " Corr. fac. for beta, ref.:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz) 2085 if (pairfunctype==7) write(ifileid,"(a,3f10.5,' :',E18.10)") " Exc.-corr. dens. for alpha, ref:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz) 2086 if (pairfunctype==8) write(ifileid,"(a,3f10.5,' :',E18.10)") " Exc.-corr. dens. for beta, ref:",refx,refy,refz,pairfunc(refx,refy,refz,inx,iny,inz) 2087 write(ifileid,"(' Source function, ref.:',3f10.5,' :',E18.10)") refx,refy,refz,srcfunc(inx,iny,inz,srcfuncmode) 2088 if (nmo/=0) write(ifileid,"(' Wavefunction value for orbital',i10,' :',E18.10)") iorbsel,fmo(inx,iny,inz,iorbsel) 2089 if (iALIEdecomp==0) then 2090 write(ifileid,"(' Average local ionization energy:',E18.10)") avglocion(inx,iny,inz) 2091 else if (iALIEdecomp==1) then 2092 call avglociondecomp(ifileid,inx,iny,inz) 2093 end if 2094 write(ifileid,"(' User-defined real space function:',E18.10)") userfunc(inx,iny,inz) 2095 fesptmp=nucesp(inx,iny,inz) 2096 if (ifiletype==4) then 2097 write(ifileid,"(' ESP from atomic charges:',E18.10)") fesptmp 2098 else 2099 write(ifileid,"(' ESP from nuclear charges:',E18.10)") fesptmp 2100 end if 2101 if (ishowptESP==1) then 2102 fesptmpelec=eleesp(inx,iny,inz) 2103 write(ifileid,"(' ESP from electrons:',E18.10)") fesptmpelec 2104 write(ifileid,"(' Total ESP:',E18.10,' a.u. (',E14.7,' J/C,',E14.7,' kcal/mol)')") fesptmpelec+fesptmp,(fesptmpelec+fesptmp)*27.2113838D0,(fesptmpelec+fesptmp)*au2kcal 2105 end if 2106 write(ifileid,*) 2107 if (ifuncsel==1) then 2108 write(ifileid,*) "Note: Below information are for electron density" 2109 funchess=elehess 2110 funcgrad=elegrad 2111 else 2112 if (ifuncsel==3) write(ifileid,*) "Note: Below information are for Laplacian of electron density" 2113 if (ifuncsel==4) write(ifileid,*) "Note: Below information are for value of orbital wavefunction" 2114 if (ifuncsel==9) write(ifileid,*) "Note: Below information are for electron localization function" 2115 if (ifuncsel==10) write(ifileid,*) "Note: Below information are for localized orbital locator" 2116 if (ifuncsel==12) write(ifileid,*) "Note: Below information are for total ESP" 2117 if (ifuncsel==100) write(ifileid,*) "Note: Below information are for user-defined real space function" 2118 call gencalchessmat(2,ifuncsel,inx,iny,inz,funcvalue,funcgrad,funchess) 2119 end if 2120 write(ifileid,*) 2121 write(ifileid,*) "Components of gradient in x/y/z are:" 2122 write(ifileid,"(3E18.10)") funcgrad(1),funcgrad(2),funcgrad(3) 2123 write(ifileid,"(' Norm of gradient is:',E18.10)") dsqrt(sum(funcgrad**2)) 2124 write(ifileid,*) 2125 write(ifileid,*) "Components of Laplacian in x/y/z are:" 2126 write(ifileid,"(3E18.10)") funchess(1,1),funchess(2,2),funchess(3,3) 2127 write(ifileid,"(' Total:',E18.10)") funchess(1,1)+funchess(2,2)+funchess(3,3) 2128 write(ifileid,*) 2129 write(ifileid,*) "Hessian matrix:" 2130 write(ifileid,"(3E18.10)") funchess 2131 call diagmat(funchess,eigvecmat,eigval,300,1D-12) 2132 write(ifileid,"(' Eigenvalues of Hessian:',3E18.10)") eigval(1:3) 2133 write(ifileid,*) "Eigenvectors(columns) of Hessian:" 2134 write(ifileid,"(3E18.10)") ((eigvecmat(i,j),j=1,3),i=1,3) 2135 write(ifileid,"(' Determinant of Hessian:',D18.10)") detmat(funchess) 2136 if (ifuncsel==1) then !Output ellipticity for rho 2137 call sortr8(eigval) 2138 eigmax=eigval(3) 2139 eigmed=eigval(2) 2140 eigmin=eigval(1) 2141 write(ifileid,"(a,f12.6)") " Ellipticity of electron density:",eigmin/eigmed-1 2142 write(ifileid,"(a,f12.6)") " eta index:",abs(eigmin)/eigmax 2143 end if 2144 2145else !Only loaded structure, use YWT promolecule density 2146 if (ifiletype==4) then 2147 write(ifileid,"(' ESP from atomic charges:',E18.10)") nucesp(inx,iny,inz) 2148 else 2149 write(ifileid,"(' ESP from nuclear charges:',E18.10)") nucesp(inx,iny,inz) 2150 end if 2151 write(ifileid,*) 2152 call calchessmat_prodens(inx,iny,inz,elerho,elegrad,elehess) 2153 write(ifileid,"(a)") " Note: The loaded file does not contain wavefunction information, below results are evaluated from promolecule density" 2154 write(ifileid,"(' Density of electrons:',E18.10)") elerho 2155 write(ifileid,"(' Reduced density gradient:',E18.10)") RDGprodens(inx,iny,inz) 2156 write(ifileid,"(' Sign(lambda2)*rho:',E18.10)") signlambda2rho_prodens(inx,iny,inz) 2157 write(ifileid,"(' User-defined real space function:',E18.10)") userfunc(inx,iny,inz) 2158 write(ifileid,*) 2159 write(ifileid,*) "Components of gradient in x/y/z are:" 2160 write(ifileid,"(3E18.10)") elegrad(1),elegrad(2),elegrad(3) 2161 write(ifileid,"(' Norm of gradient is:',E18.10)") dsqrt(sum(elegrad**2)) 2162 write(ifileid,*) 2163 write(ifileid,*) "Components of Laplacian in x/y/z are:" 2164 write(ifileid,"(3E18.10)") elehess(1,1),elehess(2,2),elehess(3,3) 2165 write(ifileid,"(' Total:',E18.10)") elehess(1,1)+elehess(2,2)+elehess(3,3) 2166 write(ifileid,*) 2167 write(ifileid,*) "Hessian matrix:" 2168 write(ifileid,"(3E18.10)") elehess 2169 call diagmat(elehess,eigvecmat,eigval,300,1D-12) 2170 write(ifileid,"(' Eigenvalues of Hessian:',3E18.10)") eigval(1:3) 2171 write(ifileid,*) "Eigenvectors(columns) of Hessian:" 2172 write(ifileid,"(3E18.10)") ((eigvecmat(i,j),j=1,3),i=1,3) 2173end if 2174end subroutine 2175 2176 2177 2178!!!------- Decompose property at a point as contribution from various orbitals 2179subroutine decompptprop(x,y,z) 2180use defvar 2181use util 2182use function 2183implicit real*8(a-h,o-z) 2184character c2000tmp*2000 2185real*8 MOocctmp(nmo) 2186integer,allocatable :: orbidx(:) 2187write(*,*) "Select function to be studied" 2188call funclist 2189read(*,*) ifunc 2190write(*,*) "Input range of orbitals, e.g. 3,6-8,12-15" 2191read(*,"(a)") c2000tmp 2192call str2arr(c2000tmp,norb) 2193allocate(orbidx(norb)) 2194call str2arr(c2000tmp,norb,orbidx) 2195MOocctmp=MOocc 2196sumval=0 2197do itmp=1,norb 2198 iorb=orbidx(itmp) 2199 MOocc=0 2200 MOocc(iorb)=MOocctmp(iorb) 2201 tmpval=calcfuncall(ifunc,x,y,z) 2202 write(*,"(' Contribution from orbital',i6,' (occ=',f10.6,'):',1D16.8,' a.u.')") iorb,MOocc(iorb),tmpval 2203 sumval=sumval+tmpval 2204end do 2205write(*,"(' Summing up above values:',1D16.8)") sumval 2206MOocc=MOocctmp 2207end subroutine 2208 2209 2210 2211 2212!!------------------ Set up grid setting 2213!If ienableloadextpt==1, then show the option used to load external points 2214!If igridsel==100, that means user didn't set up grid here but choose to load a set of point coordinates from external plain text file 2215subroutine setgrid(ienableloadextpt,igridsel) 2216use defvar 2217implicit real*8 (a-h,o-z) 2218real*8 molxlen,molylen,molzlen,tmpx,tmpy,tmpz 2219character*200 cubefilename,pointfilename 2220character c80*80 2221integer ienableloadextpt 2222logical filealive 2223ntotlow=125000 2224ntotmed=512000 2225ntothigh=1728000 2226do while(.true.) 2227 write(*,*) "Please select a method to set up grid" 2228 write(*,"(a,f10.6,a)") " -10 Set extension distance of grid range for mode 1~4, current:",aug3D," Bohr" 2229 write(*,*) "1 Low quality grid , covering whole system, about 125000 points in total" 2230 write(*,*) "2 Medium quality grid, covering whole system, about 512000 points in total" 2231 write(*,*) "3 High quality grid , covering whole system, about 1728000 points in total" 2232 write(*,*) "4 Input the number of points or grid spacing in X,Y,Z, covering whole system" 2233 write(*,*) "5 Input original point, translation vector and the number of points" 2234 write(*,*) "6 Input center coordinate, number of points and extension distance" 2235 write(*,*) "7 The same as 6, but input two atoms, the midpoint will be defined as center" 2236 write(*,*) "8 Use grid setting of another cube file" 2237 if (ienableloadextpt==1) write(*,*) "100 Load a set of points from external file" 2238 read(*,*) igridsel 2239 if (igridsel/=-10) exit 2240 write(*,*) "Input extension distance (Bohr) e.g. 6.5" 2241 read(*,*) aug3D 2242end do 2243 2244if (igridsel==100) then !Load points rather than set up grid 2245 write(*,*) "Input the path of the file containing points, e.g. c:\ltwd.txt" 2246 write(*,*) "Note: See program manual for the format of the file" 2247 do while(.true.) 2248 read(*,"(a)") pointfilename 2249 inquire(file=pointfilename,exist=filealive) 2250 if (filealive) then 2251 open(10,file=pointfilename,status="old") 2252 read(10,*) numextpt 2253 write(*,"(a,i10,a)") ' There are',numextpt,' points' 2254 if (allocated(extpt)) deallocate(extpt) 2255 allocate(extpt(numextpt,4)) 2256 do itmp=1,numextpt 2257 read(10,*) extpt(itmp,1:3) 2258 end do 2259 close(10) 2260 exit 2261 else 2262 write(*,*) "Error: File cannot be found, input again" 2263 end if 2264 end do 2265 write(*,*) "Please wait..." 2266else 2267 molxlen=(maxval(a%x)-minval(a%x))+2*aug3D 2268 molylen=(maxval(a%y)-minval(a%y))+2*aug3D 2269 molzlen=(maxval(a%z)-minval(a%z))+2*aug3D 2270 if (molxlen==0.0D0.or.molylen==0.0D0.or.molzlen==0.0D0) then !Avoid catastrophe when aug3D=0 and system is plane 2271 write(*,"(a,/)") " WARNING: The box size in one of Caresian axis is zero, & 2272 the calculation cannot be proceeded. Therefore, the size of corresponding direction is automatically set to 3 Bohr" 2273 if (molxlen==0D0) then 2274 molxlen=3D0 2275 else if (molylen==0D0) then 2276 molylen=3D0 2277 else if (molzlen==0D0) then 2278 molzlen=3D0 2279 end if 2280 end if 2281 if (igridsel==1.or.igridsel==2.or.igridsel==3) then 2282 if (igridsel==1) dx=(molxlen*molylen*molzlen/dfloat(ntotlow))**(1.0D0/3.0D0) 2283 if (igridsel==2) dx=(molxlen*molylen*molzlen/dfloat(ntotmed))**(1.0D0/3.0D0) 2284 if (igridsel==3) dx=(molxlen*molylen*molzlen/dfloat(ntothigh))**(1.0D0/3.0D0) 2285 dy=dx 2286 dz=dx 2287 nx=nint(molxlen/dx)+1 2288 ny=nint(molylen/dy)+1 2289 nz=nint(molzlen/dz)+1 2290 orgx=minval(a%x)-aug3D 2291 orgy=minval(a%y)-aug3D 2292 orgz=minval(a%z)-aug3D 2293 else if (igridsel==4) then 2294 write(*,*) "Input the number of grid points in X,Y,Z direction, e.g. 139,59,80" 2295 write(*,"(a)") " or input the grid spacing (bohr) in X,Y,Z direction, e.g. 0.05,0.08,0.08 (if only input one value, it will be used for all directions)" 2296 read(*,"(a)") c80 2297 if (index(c80,'.')/=0) then 2298 if (index(c80,',')/=0) then 2299 read(c80,*) dx,dy,dz 2300 else 2301 read(c80,*) tmp 2302 dx=tmp 2303 dy=tmp 2304 dz=tmp 2305 end if 2306 nx=molxlen/dx+1 2307 ny=molylen/dy+1 2308 nz=molzlen/dz+1 2309 else 2310 read(c80,*) nx,ny,nz 2311 dx=molxlen/(nx-1) 2312 dy=molylen/(ny-1) 2313 dz=molzlen/(nz-1) 2314 end if 2315 orgx=minval(a%x)-aug3D 2316 orgy=minval(a%y)-aug3D 2317 orgz=minval(a%z)-aug3D 2318 else if (igridsel==5) then 2319 write(*,*) "Input X,Y,Z coordinate of original point (Bohr) e.g. 0.1,4,-1" 2320 read(*,*) orgx,orgy,orgz 2321 write(*,*) "Input X,Y,Z component of translation vector (Bohr) e.g. 0.1,0.1,0.15" 2322 read(*,*) dx,dy,dz 2323 write(*,*) "Input the number of points in X,Y,Z direction e.g. 139,59,80" 2324 read(*,*) nx,ny,nz 2325 else if (igridsel==6.or.igridsel==7) then 2326 if (igridsel==6) then 2327 write(*,*) "Input X,Y,Z coordinate of center (Angstrom)" 2328 read(*,*) cenx,ceny,cenz 2329 cenx=cenx/b2a 2330 ceny=ceny/b2a 2331 cenz=cenz/b2a 2332 else if (igridsel==7) then 2333 write(*,*) "Input index of the two atoms e.g. 2,5" 2334 write(*,*) "If the two indices are identical, box center will be placed at the nucleus" 2335 read(*,*) indatm1,indatm2 2336 cenx=(a(indatm1)%x+a(indatm2)%x)/2.0D0 2337 ceny=(a(indatm1)%y+a(indatm2)%y)/2.0D0 2338 cenz=(a(indatm1)%z+a(indatm2)%z)/2.0D0 2339 end if 2340 write(*,*) "Input the number of points in X,Y,Z direction e.g. 40,40,25" 2341 read(*,*) nx,ny,nz 2342 write(*,*) "Input the extended distance in X,Y,Z direction (Bohr) e.g. 4.0,4.0,6.5" 2343 read(*,*) aug3Dx,aug3Dy,aug3Dz 2344 orgx=cenx-aug3Dx 2345 orgy=ceny-aug3Dy 2346 orgz=cenz-aug3Dz 2347 dx=aug3Dx*2.0D0/(nx-1) 2348 dy=aug3Dy*2.0D0/(ny-1) 2349 dz=aug3Dz*2.0D0/(nz-1) 2350 else if (igridsel==8) then 2351 write(*,*) "Input filename of a cube file" 2352 do while(.true.) 2353 read(*,"(a)") cubefilename 2354 inquire(file=cubefilename,exist=filealive) 2355 if (filealive) then 2356 open(10,file=cubefilename,status="old") 2357 read(10,*) 2358 read(10,*) 2359 read(10,*) nouse,orgx,orgy,orgz 2360 read(10,*) nx,dx 2361 read(10,*) ny,rnouse,dy 2362 read(10,*) nz,rnouse,rnouse,dz 2363 close(10) 2364 exit 2365 else 2366 write(*,*) "Error: File cannot be found, input again" 2367 end if 2368 end do 2369 end if 2370 endx=orgx+dx*(nx-1) 2371 endy=orgy+dy*(ny-1) 2372 endz=orgz+dz*(nz-1) 2373 write(*,"(' Coordinate of origin in X,Y,Z is',3f12.6,' Bohr')") orgx,orgy,orgz 2374 write(*,"(' Coordinate of end point in X,Y,Z is',3f12.6,' Bohr')") endx,endy,endz 2375 write(*,"(' Grid spacing in X,Y,Z is',3f12.6,' Bohr')") dx,dy,dz 2376 write(*,"(' The number of points in X,Y,Z is',3i5,' Total:',i12)") nx,ny,nz,nx*ny*nz 2377end if 2378end subroutine 2379 2380 2381 2382!!---- Set up grid setting with fixed grid spacing, similar to setgridforbasin, but not so complicated, thus may be used for other subroutine 2383subroutine setgridfixspc 2384use defvar 2385implicit real*8 (a-h,o-z) 2386real*8 :: molxlen,molylen,molzlen 2387real*8 :: spclowqual=0.2D0,spcmedqual=0.1D0,spchighqual=0.06D0,spclunaqual=0.04D0 2388character c80tmp*80,cubefilename*200 2389do while(.true.) 2390 orgx=minval(a%x)-aug3D 2391 orgy=minval(a%y)-aug3D 2392 orgz=minval(a%z)-aug3D 2393 endx=maxval(a%x)+aug3D 2394 endy=maxval(a%y)+aug3D 2395 endz=maxval(a%z)+aug3D 2396 molxlen=endx-orgx 2397 molylen=endy-orgy 2398 molzlen=endz-orgz 2399 ntotlow=(nint(molxlen/spclowqual)+1)*(nint(molylen/spclowqual)+1)*(nint(molzlen/spclowqual)+1) 2400 ntotmed=(nint(molxlen/spcmedqual)+1)*(nint(molylen/spcmedqual)+1)*(nint(molzlen/spcmedqual)+1) 2401 ntothigh=(nint(molxlen/spchighqual)+1)*(nint(molylen/spchighqual)+1)*(nint(molzlen/spchighqual)+1) 2402 ntotluna=(nint(molxlen/spclunaqual)+1)*(nint(molylen/spclunaqual)+1)*(nint(molzlen/spclunaqual)+1) 2403 2404 write(*,*) "Please select a method for setting up grid" 2405 write(*,"(a,f10.5,a)") " -10 Set grid extension distance for mode 1~6, current:",aug3D," Bohr" 2406 write(*,"(a,f4.2,a,i14)") " 1 Low quality grid, spacing=",spclowqual," Bohr, number of grids: ",ntotlow 2407 write(*,"(a,f4.2,a,i14)") " 2 Medium quality grid, spacing=",spcmedqual," Bohr, number of grids: ",ntotmed 2408 write(*,"(a,f4.2,a,i14)") " 3 High quality grid, spacing=",spchighqual," Bohr, number of grids: ",ntothigh 2409 write(*,"(a,f4.2,a,i14)") " 4 Lunatic quality grid, spacing=",spclunaqual," Bohr, number of grids:",ntotluna 2410 write(*,*) "5 Only input grid spacing, automatically set other parameters" 2411 write(*,*) "6 Only input the number of points in X,Y,Z, automatically set other parameters" 2412 write(*,*) "7 Input original point, translation vector and the number of points" 2413 write(*,*) "8 Set center position, grid spacing and box length" 2414 write(*,*) "9 Use grid setting of another cube file" 2415 read(*,*) igridsel 2416 if (igridsel/=-10) then 2417 exit 2418 else 2419 write(*,*) "Input extension distance (Bohr) e.g. 6.5" 2420 read(*,*) aug3D 2421 end if 2422end do 2423 2424!Note: orgx,orgy,orgz,endx,endy,endz as well as molx/y/zlen for igridsel==1~6 have already been set above 2425if (igridsel==1.or.igridsel==2.or.igridsel==3.or.igridsel==4.or.igridsel==5) then 2426 if (igridsel==1) dx=spclowqual 2427 if (igridsel==2) dx=spcmedqual 2428 if (igridsel==3) dx=spchighqual 2429 if (igridsel==4) dx=spclunaqual 2430 if (igridsel==5) then 2431 write(*,*) "Input the grid spacing (bohr) e.g. 0.08" 2432 read(*,*) dx 2433 end if 2434 dy=dx 2435 dz=dx 2436 nx=nint(molxlen/dx)+1 2437 ny=nint(molylen/dy)+1 2438 nz=nint(molzlen/dz)+1 2439else if (igridsel==6) then 2440 write(*,*) "Input the number of grid points in X,Y,Z direction e.g. 139,59,80" 2441 read(*,*) nx,ny,nz 2442 dx=molxlen/(nx-1) 2443 dy=molylen/(ny-1) 2444 dz=molzlen/(nz-1) 2445else if (igridsel==7) then 2446 write(*,*) "Input X,Y,Z coordinate of original point (Bohr) e.g. 0.1,4,-1" 2447 read(*,*) orgx,orgy,orgz 2448 write(*,*) "Input X,Y,Z component of translation vector (Bohr) e.g. 0.1,0.1,0.15" 2449 read(*,*) dx,dy,dz 2450 write(*,*) "Input the number of points in X,Y,Z direction e.g. 139,59,80" 2451 read(*,*) nx,ny,nz 2452 endx=orgx+dx*(nx-1) 2453 endy=orgy+dy*(ny-1) 2454 endz=orgz+dz*(nz-1) 2455else if (igridsel==8) then 2456 write(*,*) "Input X,Y,Z coordinate of box center (in Angstrom)" 2457 write(*,*) "or input such as a8 to take the coordinate of atom 8 as box center" 2458 write(*,*) "or input such as a3,a7 to take the midpoint of atom 3 and atom 7 as box center" 2459 read(*,"(a)") c80tmp 2460 if (c80tmp(1:1)=='a') then 2461 do ich=1,len_trim(c80tmp) 2462 if (c80tmp(ich:ich)==',') exit 2463 end do 2464 if (ich==len_trim(c80tmp)+1) then 2465 read(c80tmp(2:),*) itmp 2466 cenx=a(itmp)%x 2467 ceny=a(itmp)%y 2468 cenz=a(itmp)%z 2469 else 2470 read(c80tmp(2:ich-1),*) itmp 2471 read(c80tmp(ich+2:),*) jtmp 2472 cenx=(a(itmp)%x+a(jtmp)%x)/2D0 2473 ceny=(a(itmp)%y+a(jtmp)%y)/2D0 2474 cenz=(a(itmp)%z+a(jtmp)%z)/2D0 2475 end if 2476 else 2477 read(c80tmp,*) cenx,ceny,cenz 2478 cenx=cenx/b2a 2479 ceny=ceny/b2a 2480 cenz=cenz/b2a 2481 end if 2482 write(*,*) "Input the grid spacing (bohr) e.g. 0.08" 2483 read(*,*) dx 2484 dy=dx 2485 dz=dx 2486 write(*,*) "Input the box lengths in X,Y,Z direction (Bohr) e.g. 8.0,8.0,13.5" 2487 read(*,*) molxlen,molylen,molzlen 2488 orgx=cenx-molxlen/2D0 2489 orgy=ceny-molylen/2D0 2490 orgz=cenz-molzlen/2D0 2491 endx=orgx+molxlen 2492 endy=orgy+molylen 2493 endz=orgz+molzlen 2494 nx=nint(molxlen/dx)+1 2495 ny=nint(molylen/dy)+1 2496 nz=nint(molzlen/dz)+1 2497else if (igridsel==9) then 2498 write(*,*) "Input filename of a cube file" 2499 do while(.true.) 2500 read(*,"(a)") cubefilename 2501 inquire(file=cubefilename,exist=alive) 2502 if (alive) then 2503 open(10,file=cubefilename,status="old") 2504 read(10,*) 2505 read(10,*) 2506 read(10,*) nouse,orgx,orgy,orgz 2507 read(10,*) nx,dx 2508 read(10,*) ny,rnouse,dy 2509 read(10,*) nz,rnouse,rnouse,dz 2510 close(10) 2511 exit 2512 else 2513 write(*,*) "Error: File cannot be found, input again" 2514 end if 2515 end do 2516 endx=orgx+dx*(nx-1) 2517 endy=orgy+dy*(ny-1) 2518 endz=orgz+dz*(nz-1) 2519end if 2520 2521write(*,"(' Coordinate of origin in X,Y,Z is ',3f12.6)") orgx,orgy,orgz 2522write(*,"(' Coordinate of end point in X,Y,Z is',3f12.6)") endx,endy,endz 2523write(*,"(' Spacing in X,Y,Z is',3f11.6)") dx,dy,dz 2524write(*,"(' Number of points in X,Y,Z is',3i5,' Total',i10)") nx,ny,nz,nx*ny*nz 2525end subroutine 2526 2527 2528!!!------------------------- Delete virtual orbitals higher than LUMO+10 for HF/DFT wavefunctions 2529subroutine delvirorb(infomode) 2530use defvar 2531implicit real*8 (a-h,o-z) 2532integer :: infomode,nvirsave=10 !Lowest nvirsave virtual orbitals will be reserved 2533if (ifiletype/=1.and.ifiletype/=9) return !Only works for fch and molden 2534if (idelvirorb==0) return 2535if (iuserfunc==24) return !linear response kernel require all orbital information 2536if (imodwfn==1) return !Don't make things more complicated! 2537!This routine doesn't work for post-HF instances 2538if (wfntype==0.or.wfntype==2) then !RHF, ROHF 2539 if (nmo<=naelec+nvirsave) return 2540 nmo=naelec+10 !Simply shield those virtual orbitals 2541else if (wfntype==1) then !Perserve up to LUMO+10 for alpha, and identical number of orbitals for beta 2542 if (nmo/2<=naelec+nvirsave) return !naelec is always >= nbelec 2543 nperserve=naelec+nvirsave 2544 !Cobasa and Cobasb are needn't to be modified 2545 co(naelec+11:naelec+nvirsave+nperserve,:)=co(nmo/2+1:nmo/2+nperserve,:) 2546 MOene(naelec+nvirsave+1:naelec+nvirsave+nperserve)=MOene(nmo/2+1:nmo/2+nperserve) 2547 MOocc(naelec+nvirsave+1:naelec+nvirsave+nperserve)=MOocc(nmo/2+1:nmo/2+nperserve) 2548 MOtype(naelec+nvirsave+1:naelec+nvirsave+nperserve)=MOtype(nmo/2+1:nmo/2+nperserve) 2549 nmo=2*nperserve 2550end if 2551imodwfn=1 !Will not call this routine again 2552if (infomode==1.and.(wfntype==0.or.wfntype==1.or.wfntype==2)) then 2553 write(*,"(a)") " Note: Virtual orbitals higher than LUMO+10 have been discarded for saving computational time" 2554 write(*,*) 2555end if 2556end subroutine 2557 2558 2559!!!-------- imode=1: Convert unit of grid/plane parameters from Bohr to Angstrom. =2: Convert them back 2560subroutine convgridlenunit(imode) 2561use defvar 2562implicit none 2563integer imode 2564real*8 scll 2565if (imode==1) scll=b2a 2566if (imode==2) scll=1/b2a 2567orgx=orgx*scll 2568orgy=orgy*scll 2569orgz=orgz*scll 2570orgx2D=orgx2D*scll 2571orgy2D=orgy2D*scll 2572orgz2D=orgz2D*scll 2573endx=endx*scll 2574endy=endy*scll 2575endz=endz*scll 2576dx=dx*scll 2577dy=dy*scll 2578dz=dz*scll 2579v1x=v1x*scll 2580v1y=v1y*scll 2581v1z=v1z*scll 2582v2x=v2x*scll 2583v2y=v2y*scll 2584v2z=v2z*scll 2585a1x=a1x*scll 2586a1y=a1y*scll 2587a1z=a1z*scll 2588a2x=a2x*scll 2589a2y=a2y*scll 2590a2z=a2z*scll 2591a3x=a3x*scll 2592a3y=a3y*scll 2593a3z=a3z*scll 2594d1=d1*scll 2595d2=d2*scll 2596end subroutine 2597 2598 2599!!-------- Deallocate all arrays about wavefunction except that with _org suffix, to avoid problem when load another file 2600subroutine dealloall 2601use defvar 2602implicit real*8 (a-h,o-z) 2603if (allocated(a)) deallocate(a) 2604if (allocated(b)) deallocate(b) 2605if (allocated(CO)) deallocate(CO) 2606if (allocated(MOocc)) deallocate(MOocc) 2607if (allocated(MOsym)) deallocate(MOsym) 2608if (allocated(MOene)) deallocate(MOene) 2609if (allocated(MOtype)) deallocate(MOtype) 2610if (allocated(b_EDF)) then 2611 deallocate(CO_EDF,b_EDF) 2612 nEDFprims=0 2613 nEDFelec=0 2614end if 2615!Loaded file contains basis information 2616if (allocated(shtype)) deallocate(shtype,shcen,shcon,primshexp,primshcoeff,& 2617basshell,bascen,bastype,basstart,basend,primstart,primend,primconnorm) 2618if (allocated(CObasa)) deallocate(CObasa) 2619if (allocated(CObasb)) deallocate(CObasb) 2620if (allocated(Ptot)) deallocate(Ptot) 2621if (allocated(Palpha)) deallocate(Palpha) 2622if (allocated(Pbeta)) deallocate(Pbeta) 2623if (allocated(Sbas)) deallocate(Sbas) 2624if (allocated(Dbas)) deallocate(Dbas) 2625end subroutine 2626 2627 2628 2629 2630 2631!!------- Generate atomic/fragmental Hirshfeld weight and store it to planemat, calculate free-atom/fragmental density and store it to planemattmp 2632!The atoms in the fragment is inputted as "selatm" array, nselatm is the number of its elements 2633!if itype=1, use atomic wavefunction to calculate Hirshfeld weight, and setpromol must have been invoked; if =2, use built-in atomic density to generate it 2634subroutine genHirshplanewei(selatm,nselatm,itype) 2635use defvar 2636use function 2637implicit real*8 (a-h,o-z) 2638integer selatm(nselatm),nselatm,itype 2639if (allocated(planemat)) deallocate(planemat) 2640if (allocated(planemattmp)) deallocate(planemattmp) 2641allocate(planemat(ngridnum1,ngridnum2),planemattmp(ngridnum1,ngridnum2)) 2642planemat=0D0 2643planemattmp=0D0 2644do iatm=1,ncenter_org !Calc free atomic density of each atom, get promolecular density and Hirshfeld weight of present atom 2645 iyes=0 2646 if (any(selatm==iatm)) iyes=1 2647 if (itype==1) then 2648 call dealloall 2649 call readwfn(custommapname(iatm),1) 2650 end if 2651nthreads=getNThreads() 2652!$OMP PARALLEL DO private(i,j,rnowx,rnowy,rnowz,tmpval) shared(planemat) schedule(dynamic) NUM_THREADS(nthreads) 2653 do i=1,ngridnum1 !First calculate promolecular density and store it to planemat 2654 do j=1,ngridnum2 2655 rnowx=orgx2D+(i-1)*v1x+(j-1)*v2x 2656 rnowy=orgy2D+(i-1)*v1y+(j-1)*v2y 2657 rnowz=orgz2D+(i-1)*v1z+(j-1)*v2z 2658 if (itype==1) then 2659 tmpval=fdens(rnowx,rnowy,rnowz) 2660 else 2661 tmpval=calcatmdens(iatm,rnowx,rnowy,rnowz,0) 2662 end if 2663 planemat(i,j)=planemat(i,j)+tmpval 2664 if (iyes==1) planemattmp(i,j)=planemattmp(i,j)+tmpval 2665 end do 2666 end do 2667!$OMP END PARALLEL DO 2668end do 2669if (itype==1) then 2670 call dealloall 2671 call readinfile(firstfilename,1) !Retrieve the first loaded file(whole molecule) 2672end if 2673 2674do i=1,ngridnum1 !Calculate Hirshfeld weighting function 2675 do j=1,ngridnum2 2676 if (planemat(i,j)/=0D0) then 2677 planemat(i,j)=planemattmp(i,j)/planemat(i,j) 2678 else 2679 planemat(i,j)=0D0 2680 end if 2681 end do 2682end do 2683end subroutine 2684 2685!!------- Calculate some quantities involved in shubin's project in a plane 2686!itype=1: Calculate the sum of atomic relative Shannon entropy (namely total relative Shannon entropy) 2687!itype=2: Calculate the sum of x=[rhoA-rho0A]/rhoA 2688!itype=3: Calculate the difference between total relative Shannon entropy and deformation density 2689subroutine genentroplane(itype) 2690use defvar 2691use function 2692implicit real*8 (a-h,o-z) 2693integer itype 2694real*8 planeprodens(ngridnum1,ngridnum2),planedens(ngridnum1,ngridnum2) 2695if (allocated(planemat)) deallocate(planemat) 2696allocate(planemat(ngridnum1,ngridnum2)) 2697planeprodens=0D0 2698planemat=0D0 2699!Calculate molecular density in the plane and store it to planedens 2700nthreads=getNThreads() 2701!$OMP PARALLEL DO private(i,j,rnowx,rnowy,rnowz) shared(planedens) schedule(dynamic) NUM_THREADS(nthreads) 2702do i=1,ngridnum1 2703 do j=1,ngridnum2 2704 rnowx=orgx2D+(i-1)*v1x+(j-1)*v2x 2705 rnowy=orgy2D+(i-1)*v1y+(j-1)*v2y 2706 rnowz=orgz2D+(i-1)*v1z+(j-1)*v2z 2707 planedens(i,j)=fdens(rnowx,rnowy,rnowz) 2708 end do 2709end do 2710!$OMP END PARALLEL DO 2711do jatm=1,ncenter_org !Calculate promolecular density in the plane and store it to planeprodens 2712 call dealloall 2713 call readwfn(custommapname(jatm),1) 2714nthreads=getNThreads() 2715!$OMP PARALLEL DO private(i,j,rnowx,rnowy,rnowz) shared(planeprodens) schedule(dynamic) NUM_THREADS(nthreads) 2716 do i=1,ngridnum1 2717 do j=1,ngridnum2 2718 rnowx=orgx2D+(i-1)*v1x+(j-1)*v2x 2719 rnowy=orgy2D+(i-1)*v1y+(j-1)*v2y 2720 rnowz=orgz2D+(i-1)*v1z+(j-1)*v2z 2721 planeprodens(i,j)=planeprodens(i,j)+fdens(rnowx,rnowy,rnowz) 2722 end do 2723 end do 2724!$OMP END PARALLEL DO 2725end do 2726!Calculate Hirshfeld weight, relative Shannon entropy and x=[rhoA-rho0A]/rhoA for each atom in the plane and accumulate them to planemat 2727do jatm=1,ncenter_org !Cycle each atom, calculate its contribution in the plane 2728 call dealloall 2729 call readwfn(custommapname(jatm),1) 2730nthreads=getNThreads() 2731!$OMP PARALLEL DO private(i,j,rnowx,rnowy,rnowz,rho0A,rhoA,tmpval) shared(planemat) schedule(dynamic) NUM_THREADS(nthreads) 2732 do i=1,ngridnum1 2733 do j=1,ngridnum2 2734 rnowx=orgx2D+(i-1)*v1x+(j-1)*v2x 2735 rnowy=orgy2D+(i-1)*v1y+(j-1)*v2y 2736 rnowz=orgz2D+(i-1)*v1z+(j-1)*v2z 2737 rho0A=fdens(rnowx,rnowy,rnowz) 2738 rhoA=planedens(i,j)*rho0A/planeprodens(i,j) 2739 if (itype==1.or.itype==3) tmpval=rhoA*log(rhoA/rho0A) !Relative Shannon entropy 2740 if (itype==2) tmpval=(rhoA-rho0A)/rhoA !x=[rhoA-rho0A]/rhoA 2741 planemat(i,j)=planemat(i,j)+tmpval 2742 end do 2743 end do 2744!$OMP END PARALLEL DO 2745end do 2746call dealloall 2747call readinfile(firstfilename,1) !Retrieve the first loaded file(whole molecule) 2748if (itype==3) planemat=planemat-(planedens-planeprodens) !Diff between total relative Shannon entropy and deformation density 2749end subroutine 2750 2751 2752 2753!!----- Generate atomic Hirshfeld weight and store it to cubmat 2754!The atoms in the fragment is inputted as "selatm" array, nselatm is the number of its elements 2755!if itype=1, use atomic wavefunction to calculate Hirshfeld weight, and setpromol must have been invoked; if =2, use built-in atomic density to generate it 2756subroutine genHirshcubewei(selatm,nselatm,itype) 2757use defvar 2758use function 2759implicit real*8 (a-h,o-z) 2760integer selatm(nselatm),nselatm,itype 2761if (allocated(cubmat)) deallocate(cubmat) 2762if (allocated(cubmattmp)) deallocate(cubmattmp) 2763allocate(cubmat(nx,ny,nz),cubmattmp(nx,ny,nz)) 2764cubmat=0D0 2765cubmattmp=0D0 2766do iatm=1,ncenter_org 2767 write(*,"(' Finished',i6,' /',i6)") iatm,ncenter_org 2768 if (itype==1) then 2769 call dealloall 2770 call readwfn(custommapname(iatm),1) 2771 end if 2772nthreads=getNThreads() 2773!$OMP PARALLEL DO SHARED(cubmat,cubmattmp,ifinish) PRIVATE(i,j,k,tmpx,tmpy,tmpz,tmpval) schedule(dynamic) NUM_THREADS(nthreads) 2774 do k=1,nz !First calculate promolecular density and store it to cubmat 2775 tmpz=orgz+(k-1)*dz 2776 do j=1,ny 2777 tmpy=orgy+(j-1)*dy 2778 do i=1,nx 2779 tmpx=orgx+(i-1)*dx 2780 if (itype==1) then 2781 tmpval=fdens(tmpx,tmpy,tmpz) 2782 else 2783 tmpval=calcatmdens(iatm,tmpx,tmpy,tmpz,0) 2784 end if 2785 cubmat(i,j,k)=cubmat(i,j,k)+tmpval 2786 if (any(selatm==iatm)) cubmattmp(i,j,k)=cubmattmp(i,j,k)+tmpval 2787 end do 2788 end do 2789 end do 2790!$OMP END PARALLEL DO 2791end do 2792if (itype==1) then 2793 call dealloall 2794 call readinfile(firstfilename,1) !Retrieve the first loaded file(whole molecule) 2795end if 2796 2797do k=1,nz !Calculate Hirshfeld weighting function 2798 do j=1,ny 2799 do i=1,nx 2800 if (cubmat(i,j,k)/=0D0) then 2801 cubmat(i,j,k)=cubmattmp(i,j,k)/cubmat(i,j,k) 2802 else 2803 cubmat(i,j,k)=0D0 2804 end if 2805 end do 2806 end do 2807end do 2808end subroutine 2809 2810 2811 2812 2813!!----------- Output spherically averaged atomic radial density, then used for generating promolecular density for heavy atoms 2814subroutine sphatmraddens 2815use defvar 2816use function 2817implicit real*8 (a-h,o-z) 2818real*8,allocatable :: potx(:),poty(:),potz(:),potw(:),radpos(:),sphavgval(:) 2819truncrho=1D-8 2820rlow=0D0 2821rhigh=12 2822nsphpt=2030 2823nradpt=200 !Totally 200 radial points, but the number of point is truncated at truncrho 2824allocate(potx(nsphpt),poty(nsphpt),potz(nsphpt),potw(nsphpt),radpos(nradpt),sphavgval(nradpt)) 2825call Lebedevgen(nsphpt,potx,poty,potz,potw) 2826ifinish=0 2827iprogstp=20 2828iprogcrit=iprogstp 2829write(*,*) "Calculating..." 2830nthreads=getNThreads() 2831!$OMP PARALLEL DO SHARED(sphavgval,radpos,ifinish,iprogcrit) PRIVATE(irad,radx,radr,isph,rnowx,rnowy,rnowz,tmpval) schedule(dynamic) NUM_THREADS(nthreads) 2832do irad=1,nradpt 2833 radx=cos(irad*pi/(nradpt+1)) 2834 radr=(1+radx)/(1-radx) !Becke transform 2835 radpos(irad)=radr 2836 tmpval=0 2837 do isph=1,nsphpt 2838 rnowx=potx(isph)*radr 2839 rnowy=poty(isph)*radr 2840 rnowz=potz(isph)*radr 2841 tmpval=tmpval+fdens(rnowx,rnowy,rnowz)*potw(isph) 2842 end do 2843 sphavgval(irad)=tmpval !Spherically average density 2844 ifinish=ifinish+1 2845 if (ifinish==iprogcrit) then 2846 write(*,"(' Finished:',i6,' /',i6)") ifinish,nradpt 2847 iprogcrit=iprogcrit+iprogstp 2848 end if 2849end do 2850!$OMP END PARALLEL DO 2851open(10,file="sphavgval.txt",status="replace") 2852itmp=0 2853do irad=nradpt,1,-1 2854 if (sphavgval(irad)>truncrho) itmp=itmp+1 2855end do 2856write(10,"(a,i3,a)") "else if (iele==",a(1)%index,") then !" 2857write(10,"(' npt=',i5)") itmp 2858itmp=0 2859do irad=nradpt,1,-1 2860 if (sphavgval(irad)>truncrho) then 2861 itmp=itmp+1 2862 write(10,"(' rhoarr(',i3,')=',f25.10,'D0')") itmp,sphavgval(irad) 2863 end if 2864end do 2865close(10) 2866write(*,*) "The result has been output to sphavgval.txt in current folder" 2867write(*,*) "The second column is radial distance (Bohr), the third column is value" 2868end subroutine 2869 2870 2871 2872 2873!!--- Generate single-center integration grid for Becke's integration. Not adapted according to element. Return iradcut and gridatm 2874subroutine gen1cintgrid(gridatm,iradcut) 2875use defvar 2876implicit real*8 (a-h,o-z) 2877integer iradcut 2878real*8 potx(sphpot),poty(sphpot),potz(sphpot),potw(sphpot) 2879type(content) gridatm(radpot*sphpot) 2880call Lebedevgen(sphpot,potx,poty,potz,potw) 2881iradcut=0 !Before where the radial points will be cut 2882parm=1D0 2883do i=1,radpot !Combine spherical point&weights with second kind Gauss-Chebyshev method for radial part 2884 radx=cos(i*pi/(radpot+1)) 2885 radr=(1+radx)/(1-radx)*parm !Becke transform 2886 radw=2*pi/(radpot+1)*parm**3 *(1+radx)**2.5D0/(1-radx)**3.5D0 *4*pi 2887 gridatm( (i-1)*sphpot+1:i*sphpot )%x=radr*potx 2888 gridatm( (i-1)*sphpot+1:i*sphpot )%y=radr*poty 2889 gridatm( (i-1)*sphpot+1:i*sphpot )%z=radr*potz 2890 gridatm( (i-1)*sphpot+1:i*sphpot )%value=radw*potw 2891 if (radcut/=0D0.and.iradcut==0.and.radr<radcut) iradcut=i-1 2892end do 2893end subroutine 2894!!--- Generate Becke weight for a batch of points around iatm, sharpness parameter=3 2895!!--- Input: iatm, iradcut, gridatm Return: beckeweigrid 2896subroutine gen1cbeckewei(iatm,iradcut,gridatm,beckeweigrid) 2897use defvar 2898implicit real*8 (a-h,o-z) 2899integer iatm,iradcut 2900real*8 beckeweigrid(radpot*sphpot),smat(ncenter,ncenter),Pvec(ncenter) 2901type(content) gridatm(radpot*sphpot) 2902nthreads=getNThreads() 2903!$OMP parallel do shared(beckeweigrid) private(i,rnowx,rnowy,rnowz,smat,ii,ri,jj,rj,rmiu,chi,uij,aij,tmps,Pvec) num_threads(nthreads) schedule(dynamic) 2904do i=1+iradcut*sphpot,radpot*sphpot 2905 smat=1D0 2906 rnowx=gridatm(i)%x 2907 rnowy=gridatm(i)%y 2908 rnowz=gridatm(i)%z 2909 do ii=1,ncenter 2910 ri=dsqrt( (rnowx-a(ii)%x)**2+(rnowy-a(ii)%y)**2+(rnowz-a(ii)%z)**2 ) 2911 do jj=1,ncenter 2912 if (ii==jj) cycle 2913 rj=dsqrt( (rnowx-a(jj)%x)**2+(rnowy-a(jj)%y)**2+(rnowz-a(jj)%z)**2 ) 2914 rmiu=(ri-rj)/distmat(ii,jj) 2915 !Adjust for heteronuclear 2916 chi=covr_tianlu(a(ii)%index)/covr_tianlu(a(jj)%index) 2917 uij=(chi-1)/(chi+1) 2918 aij=uij/(uij**2-1) 2919 if (aij>0.5D0) aij=0.5D0 2920 if (aij<-0.5D0) aij=-0.5D0 2921 rmiu=rmiu+aij*(1-rmiu**2) 2922 tmps=rmiu 2923 do iter=1,3 2924 tmps=1.5D0*(tmps)-0.5D0*(tmps)**3 2925 end do 2926 smat(ii,jj)=0.5D0*(1-tmps) 2927 end do 2928 end do 2929 Pvec=1D0 2930 do ii=1,ncenter 2931 Pvec=Pvec*smat(:,ii) 2932 end do 2933 beckeweigrid(i)=Pvec(iatm)/sum(Pvec) 2934end do 2935!$OMP end parallel do 2936end subroutine 2937 2938 2939 2940!!--------- Convert 1RDM in MO basis outputted by MRCC program to natural orbitals 2941!In CCDENSITIES, the density matrix is represented in MO basis 2942!When frozen core is enabled, the indices are counted from the first correlated orbital 2943subroutine MRCC_gennatorb 2944use defvar 2945use util 2946character c200tmp*200 2947real*8,allocatable :: eigvecmat(:,:),eigvalarr(:),tmparr(:) 2948do while(.true.) 2949 write(*,*) "Input the path of CCDENSITIES, e.g. C:\lovelive\CCDENSITIES" 2950! c200tmp="D:\CM\my_program\Multiwfn\x\MRCCdens\HF_m3_CCSD\CCDENSITIES" 2951 read(*,"(a)") c200tmp 2952 inquire(file=c200tmp,exist=alive) 2953 if (alive) exit 2954 write(*,*) "Cannot find the file, input again" 2955end do 2956write(*,*) 2957write(*,*) "Input the number of frozen orbitals, e.g. 3" 2958write(*,*) "If no orbitals are frozen, simply input 0" 2959write(*,"(a)") " PS: For unrestricted reference, if you input n, the n lowest alpha and n lowest beta MOs will be regarded as frozen" 2960read(*,*) nfrz 2961write(*,*) "Please wait..." 2962if (wfntype==0) then !RHF reference 2963 open(10,file=c200tmp,status="old") 2964 Ptot=0 2965 do while(.true.) 2966 read(10,*,iostat=ierror) tmp,i,j,k,l 2967 if (ierror/=0) exit 2968 if (k==0.and.l==0) then !Only load 1RDM 2969 Ptot(i+nfrz,j+nfrz)=tmp 2970 Ptot(j+nfrz,i+nfrz)=tmp 2971 end if 2972 end do 2973 close(10) 2974 do ifrz=1,nfrz 2975 Ptot(ifrz,ifrz)=2D0 2976 end do 2977 allocate(eigvecmat(nbasis,nbasis),eigvalarr(nbasis),tmparr(nbasis)) 2978 call diagsymat(Ptot,eigvecmat,eigvalarr,istat) 2979 MOocc=eigvalarr 2980 !Currently the occupation is from low to high, now invert the sequence 2981 do i=1,int(nmo/2D0) 2982 idx=i 2983 jdx=nmo+1-i 2984 tmp=MOocc(idx) 2985 MOocc(idx)=MOocc(jdx) 2986 MOocc(jdx)=tmp 2987 tmparr=eigvecmat(:,idx) 2988 eigvecmat(:,idx)=eigvecmat(:,jdx) 2989 eigvecmat(:,jdx)=tmparr 2990 end do 2991 CObasa=matmul(CObasa,eigvecmat) 2992 wfntype=3 2993 write(*,*) "Occupation number:" 2994 write(*,"(6f12.6)") MOocc 2995else if (wfntype==1) then !UHF reference 2996 !In CCDENSITIES, the sequence is: 2997 !2RDM-alpha 2998 ! 0.00000000000000000000E+00 0 0 0 0 2999 !2RDM-beta 3000 ! 0.00000000000000000000E+00 0 0 0 0 3001 !Unknown 3002 ! 0.00000000000000000000E+00 0 0 0 0 3003 !1RDM-alpha 3004 ! 0.00000000000000000000E+00 0 0 0 0 3005 !1RDM-beta 3006 ! 0.00000000000000000000E+00 0 0 0 0 3007 open(10,file=c200tmp,status="old") 3008 Palpha=0 3009 Pbeta=0 3010 itime=0 3011 do while(.true.) 3012 read(10,*) tmp,i,j,k,l 3013 if (i==0.and.j==0.and.k==0.and.l==0) then 3014 itime=itime+1 3015 if (itime==5) exit 3016 cycle 3017 end if 3018 if (itime==3) then 3019 Palpha(i+nfrz,j+nfrz)=tmp 3020 Palpha(j+nfrz,i+nfrz)=tmp 3021 else if (itime==4) then 3022 Pbeta(i+nfrz,j+nfrz)=tmp 3023 Pbeta(j+nfrz,i+nfrz)=tmp 3024 end if 3025 end do 3026 close(10) 3027 do ifrz=1,nfrz 3028 Palpha(ifrz,ifrz)=1D0 3029 Pbeta(ifrz,ifrz)=1D0 3030 end do 3031 allocate(eigvecmat(nbasis,nbasis),eigvalarr(nbasis),tmparr(nbasis)) 3032 !Alpha part 3033 call diagsymat(Palpha,eigvecmat,eigvalarr,istat) 3034 MOocc(1:nbasis)=eigvalarr 3035 do i=1,int(nbasis/2D0) 3036 idx=i 3037 jdx=nbasis+1-i 3038 tmp=MOocc(idx) 3039 MOocc(idx)=MOocc(jdx) 3040 MOocc(jdx)=tmp 3041 tmparr=eigvecmat(:,idx) 3042 eigvecmat(:,idx)=eigvecmat(:,jdx) 3043 eigvecmat(:,jdx)=tmparr 3044 end do 3045 CObasa=matmul(CObasa,eigvecmat) 3046 write(*,*) "Occupation number of Alpha part:" 3047 write(*,"(6f12.6)") MOocc(1:nbasis) 3048 !Beta part 3049 call diagsymat(Pbeta,eigvecmat,eigvalarr,istat) 3050 MOocc(nbasis+1:nmo)=eigvalarr 3051 do i=1,int(nbasis/2D0) 3052 idx=nbasis+i 3053 jdx=nmo+1-i 3054 tmp=MOocc(idx) 3055 MOocc(idx)=MOocc(jdx) 3056 MOocc(jdx)=tmp 3057 tmparr=eigvecmat(:,i) 3058 eigvecmat(:,i)=eigvecmat(:,nbasis+1-i) 3059 eigvecmat(:,nbasis+1-i)=tmparr 3060 end do 3061 CObasb=matmul(CObasb,eigvecmat) 3062 write(*,*) "Occupation number of Beta part:" 3063 write(*,"(6f12.6)") MOocc(nbasis+1:nmo) 3064 wfntype=4 3065end if 3066 3067call genP 3068MOene=0 3069write(*,*) "Done! Basis function information now correspond to natural orbital cases" 3070write(*,"(a)") " Note: If next you would like to analyze real space functions, you should export .molden file, & 3071and then reload it, so that GTF information will also correspond to natural orbitals" 3072end subroutine 3073 3074 3075 3076 3077!!----------- Generate spherical harmonic -> Cartesian basis function conversion table for d,f,g,h. 3078!iprog=1: for readfch; iprog=2: for readmolden 3079!The table comes from IJQC,54,83, which is used by Gaussian 3080!The sequence of d and f shell is also identical to .molden convention, but for g, another conversion table is used, & 3081!since in Multiwfn g cartesian shell starts from ZZZZ, but that of .molden starts from xxxx 3082subroutine gensphcartab(iprog,matd,matf,matg,math) 3083real*8 matd(6,5),matf(10,7),matg(15,9),math(21,11) 3084integer iprog 3085matd=0D0 3086matf=0D0 3087matg=0D0 3088math=0D0 3089! From 5D: D 0,D+1,D-1,D+2,D-2 3090! To 6D: 1 2 3 4 5 6 3091! XX,YY,ZZ,XY,XZ,YZ 3092! 3093! D0=-0.5*XX-0.5*YY+ZZ 3094matd(1:3,1)=(/ -0.5D0,-0.5D0,1D0 /) 3095! D+1=XZ 3096matd(5,2)=1D0 3097! D-1=YZ 3098matd(6,3)=1D0 3099! D+2=SQRT(3)/2*(XX-YY) 3100matd(1:2,4)=(/ sqrt(3D0)/2D0,-sqrt(3D0)/2D0 /) 3101! D-2=XY 3102matd(4,5)=1D0 3103 3104! From 7F: F 0,F+1,F-1,F+2,F-2,F+3,F-3 3105! To 10F: 1 2 3 4 5 6 7 8 9 10 3106! XXX,YYY,ZZZ,XYY,XXY,XXZ,XZZ,YZZ,YYZ,XYZ (Gaussian sequence, not identical to Multiwfn) 3107! 3108! F 0=-3/(2*��5)*(XXZ+YYZ)+ZZZ 3109matf(3,1)=1D0 3110matf(6,1)=-1.5D0/sqrt(5D0) 3111matf(9,1)=-1.5D0/sqrt(5D0) 3112! F+1=-��(3/8)*XXX-��(3/40)*XYY+��(6/5)*XZZ 3113matf(1,2)=-sqrt(3D0/8D0) 3114matf(4,2)=-sqrt(3D0/40D0) 3115matf(7,2)=sqrt(6D0/5D0) 3116! F-1=-��(3/40)*XXY-��(3/8)*YYY+��(6/5)*YZZ 3117matf(2,3)=-sqrt(3D0/8D0) 3118matf(5,3)=-sqrt(3D0/40D0) 3119matf(8,3)=sqrt(6D0/5D0) 3120! F+2=��3/2*(XXZ-YYZ) 3121matf(6,4)=sqrt(3D0)/2D0 3122matf(9,4)=-sqrt(3D0)/2D0 3123! F-2=XYZ 3124matf(10,5)=1D0 3125! F+3=��(5/8)*XXX-3/��8*XYY 3126matf(1,6)=sqrt(5D0/8D0) 3127matf(4,6)=-3D0/sqrt(8D0) 3128! F-3=3/��8*XXY-��(5/8)*YYY 3129matf(2,7)=-sqrt(5D0/8D0) 3130matf(5,7)=3D0/sqrt(8D0) 3131 3132if (iprog==1) then !for .fch 3133 ! From 9G: G 0,G+1,G-1,G+2,G-2,G+3,G-3,G+4,G-4 3134 ! To 15G: 1 2 3 4 5 6 7 8 3135 ! ZZZZ,YZZZ,YYZZ,YYYZ,YYYY,XZZZ,XYZZ,XYYZ 3136 ! 9 10 11 12 13 14 15 3137 ! XYYY,XXZZ,XXYZ,XXYY,XXXZ,XXXY,XXXX 3138 ! 3139 !G 0=ZZZZ+3/8*(XXXX+YYYY)-3*��(3/35)*(XXZZ+YYZZ-1/4*XXYY) 3140 ! matg(1,1)=1D0 3141 ! matg(3,1)=-3D0*sqrt(3D0/35D0) 3142 ! matg(5,1)=3D0/8D0 3143 ! matg(10,1)=-3D0*sqrt(3D0/35D0) 3144 ! matg(12,1)=3D0/4D0*sqrt(3D0/35D0) 3145 ! matg(15,1)=3D0/8D0 3146 ! !G+1=2*��(5/14)*XZZZ-3/2*��(5/14)*XXXZ-3/2/��14*XYYZ 3147 ! matg(6,2)=2D0*sqrt(5D0/14D0) 3148 ! matg(8,2)=-1.5D0/sqrt(14D0) 3149 ! matg(13,2)=-1.5D0*sqrt(5D0/14D0) 3150 ! !G-1=2*��(5/14)*YZZZ-3/2*��(5/14)*YYYZ-3/2/��14*XXYZ 3151 ! matg(2,3)=2D0*sqrt(5D0/14D0) 3152 ! matg(4,3)=-1.5D0*sqrt(5D0/14D0) 3153 ! matg(11,3)=-1.5D0/sqrt(14D0) 3154 ! !G+2=3*��(3/28)*(XXZZ-YYZZ)-��5/4*(XXXX-YYYY) 3155 ! matg(3,4)=-3D0*sqrt(3D0/28D0) 3156 ! matg(5,4)=sqrt(5D0)/4D0 3157 ! matg(10,4)=3D0*sqrt(3D0/28D0) 3158 ! matg(15,4)=-sqrt(5D0)/4D0 3159 ! !G-2=3/��7*XYZZ-��(5/28)*(XXXY+XYYY) 3160 ! matg(7,5)=3D0/sqrt(7D0) 3161 ! matg(9,5)=-sqrt(5D0/28D0) 3162 ! matg(14,5)=-sqrt(5D0/28D0) 3163 ! !G+3=��(5/8)*XXXZ-3/��8*XYYZ 3164 ! matg(8,6)=-3D0/sqrt(8D0) 3165 ! matg(13,6)=sqrt(5D0/8D0) 3166 ! !G-3=-��(5/8)*YYYZ+3/��8*XXYZ 3167 ! matg(4,7)=-sqrt(5D0/8D0) 3168 ! matg(11,7)=3D0/sqrt(8D0) 3169 ! !G+4=��35/8*(XXXX+YYYY)-3/4*��3*XXYY 3170 ! matg(5,8)=sqrt(35D0)/8D0 3171 ! matg(12,8)=-3D0/4D0*sqrt(3D0) 3172 ! matg(15,8)=sqrt(35D0)/8D0 3173 ! !G-4=��5/2*(XXXY-XYYY) 3174 ! matg(9,9)=-sqrt(5D0)/2D0 3175 ! matg(14,9)=sqrt(5D0)/2D0 3176else if (iprog==2) then !For .molden 3177 ! From 9G: G 0,G+1,G-1,G+2,G-2,G+3,G-3,G+4,G-4 3178 ! To 15G: 1 2 3 4 5 6 7 8 3179 ! xxxx,yyyy,zzzz,xxxy,xxxz,yyyx,yyyz,zzzx 3180 ! 9 10 11 12 13 14 15 3181 ! zzzy,xxyy,xxzz,yyzz,xxyz,yyxz,zzxy 3182 ! 3183 !G 0=ZZZZ+3/8*(XXXX+YYYY)-3*��(3/35)*(XXZZ+YYZZ-1/4*XXYY) 3184 matg(3,1)=1D0 3185 matg(1,1)=3D0/8D0 3186 matg(2,1)=3D0/8D0 3187 matg(11,1)=-3D0*sqrt(3D0/35D0) 3188 matg(12,1)=-3D0*sqrt(3D0/35D0) 3189 matg(10,1)=3D0/4D0*sqrt(3D0/35D0) 3190 !G+1=2*��(5/14)*XZZZ-3/2*��(5/14)*XXXZ-3/2/��14*XYYZ 3191 matg(8,2)=2D0*sqrt(5D0/14D0) 3192 matg(5,2)=-1.5D0*sqrt(5D0/14D0) 3193 matg(14,2)=-1.5D0/sqrt(14D0) 3194 !G-1=2*��(5/14)*YZZZ-3/2*��(5/14)*YYYZ-3/2/��14*XXYZ 3195 matg(9,3)=2D0*sqrt(5D0/14D0) 3196 matg(7,3)=-1.5D0*sqrt(5D0/14D0) 3197 matg(13,3)=-1.5D0/sqrt(14D0) 3198 !G+2=3*��(3/28)*(XXZZ-YYZZ)-��5/4*(XXXX-YYYY) 3199 matg(11,4)=3D0*sqrt(3D0/28D0) 3200 matg(12,4)=-3D0*sqrt(3D0/28D0) 3201 matg(1,4)=-sqrt(5D0)/4D0 3202 matg(2,4)=sqrt(5D0)/4D0 3203 !G-2=3/��7*XYZZ-��(5/28)*(XXXY+XYYY) 3204 matg(15,5)=3D0/sqrt(7D0) 3205 matg(4,5)=-sqrt(5D0/28D0) 3206 matg(6,5)=-sqrt(5D0/28D0) 3207 !G+3=��(5/8)*XXXZ-3/��8*XYYZ 3208 matg(5,6)=sqrt(5D0/8D0) 3209 matg(14,6)=-3D0/sqrt(8D0) 3210 !G-3=-��(5/8)*YYYZ+3/��8*XXYZ 3211 matg(7,7)=-sqrt(5D0/8D0) 3212 matg(13,7)=3D0/sqrt(8D0) 3213 !G+4=��35/8*(XXXX+YYYY)-3/4*��3*XXYY 3214 matg(1,8)=sqrt(35D0)/8D0 3215 matg(2,8)=sqrt(35D0)/8D0 3216 matg(10,8)=-3D0/4D0*sqrt(3D0) 3217 !G-4=��5/2*(XXXY-XYYY) 3218 matg(4,9)=sqrt(5D0)/2D0 3219 matg(6,9)=-sqrt(5D0)/2D0 3220end if 3221 3222! From 11H: H 0,H+1,H-1,H+2,H-2,H+3,H-3,H+4,H-4,H+5,H-5 3223! To 21H: 1 2 3 4 5 6 7 8 9 10 3224! ZZZZZ YZZZZ YYZZZ YYYZZ YYYYZ YYYYY XZZZZ XYZZZ XYYZZ XYYYZ 3225! 11 12 13 14 15 16 17 18 19 20 21 3226! XYYYY XXZZZ XXYZZ XXYYZ XXYYY XXXZZ XXXYZ XXXYY XXXXZ XXXXY XXXXX 3227! 3228!H 0=ZZZZZ-5/��21*(XXZZZ+YYZZZ)+5/8*(XXXXZ+YYYYZ)+��(15/7)/4*XXYYZ 3229math(1,1)=1D0 3230math(12,1)=-5D0/sqrt(21D0) 3231math(3,1)=-5D0/sqrt(21D0) 3232math(19,1)=5D0/8D0 3233math(5,1)=5D0/8D0 3234math(14,1)=sqrt(15D0/7D0)/4D0 3235!H+1=��(5/3)*XZZZZ-3*��(5/28)*XXXZZ-3/��28*XYYZZ+��15/8*XXXXX+��(5/3)/8*XYYYY+��(5/7)/4*XXXYY 3236math(7,2)=sqrt(5D0/3D0) 3237math(16,2)=-3D0*sqrt(5D0/28D0) 3238math(9,2)=-3D0/sqrt(28D0) 3239math(21,2)=sqrt(15D0)/8D0 3240math(11,2)=sqrt(5D0/3D0)/8D0 3241math(18,2)=sqrt(5D0/7D0)/4D0 3242!H-1=��(5/3)*YZZZZ-3*��(5/28)*YYYZZ-3/��28*XXYZZ+��15/8*YYYYY+��(5/3)/8*XXXXY+��(5/7)/4*XXYYY 3243math(2,3)=sqrt(5D0/3D0) 3244math(4,3)=-3D0*sqrt(5D0/28D0) 3245math(13,3)=-3D0/sqrt(28D0) 3246math(6,3)=sqrt(15D0)/8D0 3247math(20,3)=sqrt(5D0/3D0)/8D0 3248math(15,3)=sqrt(5D0/7D0)/4D0 3249!H+2=��5/2*(XXZZZ-YYZZZ)-��(35/3)/4*(XXXXZ-YYYYZ) 3250math(12,4)=sqrt(5D0)/2D0 3251math(3,4)=-sqrt(5D0)/2D0 3252math(19,4)=-sqrt(35D0/3D0)/4D0 3253math(5,4)=sqrt(35D0/3D0)/4D0 3254!H-2=��(5/3)*XYZZZ-��(5/12)*(XXXYZ+XYYYZ) 3255math(8,5)=sqrt(5D0/3D0) 3256math(17,5)=-sqrt(5D0/12D0) 3257math(10,5)=-sqrt(5D0/12D0) 3258!H+3=��(5/6)*XXXZZ-��(3/2)*XYYZZ-��(35/2)/8*(XXXXX-XYYYY)+��(5/6)/4*XXXYY 3259math(16,6)=sqrt(5D0/6D0) 3260math(9,6)=-sqrt(1.5D0) 3261math(21,6)=-sqrt(17.5D0)/8D0 3262math(11,6)=sqrt(17.5D0)/8D0 3263math(18,6)=sqrt(5D0/6D0)/4D0 3264!H-3=-��(5/6)*YYYZZ+��(3/2)*XXYZZ-��(35/2)/8*(XXXXY-YYYYY)-��(5/6)/4*XXYYY 3265math(4,7)=-sqrt(5D0/6D0) 3266math(13,7)=sqrt(1.5D0) 3267math(20,7)=-sqrt(17.5D0)/8D0 3268math(6,7)=sqrt(17.5D0)/8D0 3269math(15,7)=-sqrt(5D0/6D0)/4D0 3270!H+4=��35/8*(XXXXZ+YYYYZ)-3/4*��3*XXYYZ 3271math(19,8)=sqrt(35D0)/8D0 3272math(5,8)=sqrt(35D0)/8D0 3273math(14,8)=-0.75D0*sqrt(3D0) 3274!H-4=��5/2*(XXXYZ-XYYYZ) 3275math(17,9)=sqrt(5D0)/2D0 3276math(10,9)=-sqrt(5D0)/2D0 3277!H+5=3/8*��(7/2)*XXXXX+5/8*��(7/2)*XYYYY-5/4*��(3/2)*XXXYY 3278math(21,10)=3D0/8D0*sqrt(3.5D0) 3279math(11,10)=5D0/8D0*sqrt(3.5D0) 3280math(18,10)=-1.25D0*sqrt(1.5D0) 3281!H-5=3/8*��(7/2)*YYYYY+5/8*��(7/2)*XXXXY-5/4*��(3/2)*XXYYY 3282math(6,11)=3D0/8D0*sqrt(3.5D0) 3283math(20,11)=5D0/8D0*sqrt(3.5D0) 3284math(15,11)=-1.25D0*sqrt(1.5D0) 3285end subroutine 3286 3287 3288 3289!!-------- Randomly generate name of Sobereva's lover 3290subroutine mylover(outname) 3291integer,parameter :: nlovers=46 3292character*80 lovername(nlovers),outname 3293CALL RANDOM_SEED() 3294CALL RANDOM_NUMBER(tmp) 3295lovername(1)="K-ON\Mio_Akiyama" 3296lovername(2)="K-ON\Azusa_Nakano" 3297lovername(3)="EVA\Rei_Ayanami" 3298lovername(4)="Ore_no_Imoto\Black_Cat" 3299lovername(5)="Touhou_project\Ran_Yakumo" 3300lovername(6)="Haiyore!Nyaruko-san\Nyaruko" 3301lovername(7)="Bodacious_Space_Pirates\Kurihara_Chiaki" 3302lovername(8)="Otoboku\Mariya_Mikado" 3303lovername(9)="Amagami\Miya_Tachibana" 3304lovername(10)="Shakugan_no_Shana\Shana" 3305lovername(11)="Yuru_Yuri\Akari_Akaza" 3306lovername(12)="Natsuiro_Kiseki\Yuka_Hanaki" 3307lovername(13)="Love_Live!\Nico_Yazawa" 3308lovername(14)="Vocaloid\Miku_Hatsune" 3309lovername(15)="iDOLM@STER\Makoto_Kikuchi" 3310lovername(16)="Last_Exile\Dio_Eraclea" 3311lovername(17)="NHK_ni_Youkoso!\Misaki_Nakahara" 3312lovername(18)="Rio_Rainbow_Gate\Rio_Rollins" 3313lovername(19)="Blood-C\Saya_Kisaragi" 3314lovername(20)="Mahou_Shoujo_Madoka-Magica\Homura_Akemi" 3315lovername(21)="Saki\Hisa_Takei" 3316lovername(22)="Strawberry_Panic\Chikaru_Minamoto" 3317lovername(23)="Najica\Najica_Hiiragi" 3318lovername(24)="Blue_Drop\Hagino_Senkouji" 3319lovername(25)="Fate_Zero\Saber" 3320lovername(26)="Baka_to_Test_to_Shoukanjuu\Hideyoshi_Kinoshita" 3321lovername(27)="Watamote\Tomoko_Kuroki" 3322lovername(28)="Genshiken_Nidaime\Kenjirou_Hato" 3323lovername(29)="The_World_God_Only_Knows\Chihiro_Kosaka" 3324lovername(30)="Kan_Colle\Shimakaze" 3325lovername(31)="Wake_Up,Girls!\Miyu_Okamoto" 3326lovername(32)="Gokukoku\Kazumi_Schlierenzauer" 3327lovername(33)="Love_Live!\Nozomi_Tojo" 3328lovername(34)="Tokimeki_Memorial\Yuina_Himoo" 3329lovername(35)="MADLAX\MADLAX" 3330lovername(36)="Gun_Gale_Online\Kirito" 3331lovername(37)="Denkigai_No_Honyasan\Sennsei" 3332lovername(38)="Kan_Colle\Kongou" 3333lovername(39)="Plastic_Memories\Aira" 3334lovername(40)="Real_world\sell-moe-kun\" 3335lovername(41)="Sakurako-san_no_Ashimoto_ni_wa_Shitai_ga_Umatteiru\Sakurako" 3336lovername(42)="Hibike!_Euphonium\Reina_Kousaka" 3337lovername(43)="Planetarian\Yumemi_Hoshino" 3338lovername(44)="Lovelive_sunshine!!\Yoshiko_Tsushima" 3339lovername(45)="Lovelive_sunshine!!\You_Watanabe" 3340lovername(46)="Tiger_Mask_W\Miss_X" 3341!Dear Kanan, 3342! 3343!You are the only one I deeply love forever in the real world, 3344!although you can't be with me, and I am even unable to know your name and touch your finger. 3345!I believe I will never love anyone else in the rest of my life. 3346! 3347!I love your brilliant dance, your kawaii smile, your lovely double ponytail, and especially, your extremely pure and beautiful heart. 3348! 3349! ----- The author of Multiwfn, Tian Lu, 2015-May-19 3350outname=lovername(ceiling(tmp*nlovers)) 3351end subroutine 3352 3353 3354!!------------ Load parameters in settings.ini when boot up multiwfn 3355subroutine loadsetting 3356use defvar 3357use util 3358implicit real*8 (a-h,o-z) 3359character c80tmp*80,c200tmp*200,settingpath*80 3360!Set default color of atomic spheres 3361atm3Dclr(:,1)=0.85D0 3362atm3Dclr(:,2)=0.6D0 3363atm3Dclr(:,3)=0.5D0 3364atm3Dclr(0,:)=(/0.0D0, 0.5D0, 0.6D0 /) !Bq 3365atm3Dclr(1,:)=(/0.95D0, 0.95D0, 0.95D0/) !H 3366atm3Dclr(5,:)=(/0.95D0, 0.7D0, 0.7D0 /) !B 3367atm3Dclr(6,:)=(/0.85D0, 0.85D0, 0.55D0/) !C 3368atm3Dclr(7,:)=(/0.5D0, 0.5D0, 1.0D0 /) !N 3369atm3Dclr(8,:)=(/1.0D0, 0.2D0, 0.2D0 /) !O 3370atm3Dclr(9,:)=(/0.6D0, 0.9D0, 0.9D0 /) !F 3371atm3Dclr(15,:)=(/0.9D0, 0.4D0, 0.0D0 /) !P 3372atm3Dclr(16,:)=(/0.9D0, 0.7D0, 0.1D0 /) !S 3373atm3Dclr(17,:)=(/0.1D0, 0.9D0, 0.1D0 /) !Cl 3374 3375inquire(file="settings.ini",exist=alive) 3376if (alive.eqv..true.) then 3377 settingpath="settings.ini" 3378else if (alive.eqv..false.) then 3379 call getenv("Multiwfnpath",c80tmp) 3380 if (isys==1) then 3381 settingpath=trim(c80tmp)//"\settings.ini" 3382 else if (isys==2.or.isys==3) then 3383 settingpath=trim(c80tmp)//"/settings.ini" 3384 end if 3385 inquire(file=settingpath,exist=alive) 3386 if (alive.eqv..false.) then 3387 write(*,"(a)") " Warning: ""settings.ini"" was found neither in current folder nor in the path defined by ""Multiwfnpath"" & 3388 environment variable. Now using default settings instead" 3389 write(*,*) 3390 return 3391 end if 3392end if 3393 3394open(20,file=settingpath,status="old") 3395! Below are the parameters can affect calculation results 3396call loclabel(20,'iuserfunc=') 3397read(20,*) c80tmp,iuserfunc 3398call loclabel(20,'refxyz=') 3399read(20,*) c80tmp,refx,refy,refz 3400call loclabel(20,'iDFTxcsel=') 3401read(20,*) c80tmp,iDFTxcsel 3402call loclabel(20,'paircorrtype=') 3403read(20,*) c80tmp,paircorrtype 3404call loclabel(20,'pairfunctype=') 3405read(20,*) c80tmp,pairfunctype 3406call loclabel(20,'iautointgrid=') 3407read(20,*) c80tmp,iautointgrid 3408call loclabel(20,'radpot=') 3409read(20,*) c80tmp,radpot 3410call loclabel(20,'sphpot=') 3411read(20,*) c80tmp,sphpot 3412call loclabel(20,'radcut=') 3413read(20,*) c80tmp,radcut 3414call loclabel(20,'expcutoff=') 3415read(20,*) c80tmp,expcutoff 3416call loclabel(20,'espprecutoff=') 3417read(20,*) c80tmp,espprecutoff 3418call loclabel(20,'RDG_maxrho=') 3419read(20,*) c80tmp,RDG_maxrho 3420call loclabel(20,'RDGprodens_maxrho=') 3421read(20,*) c80tmp,RDGprodens_maxrho 3422call loclabel(20,'ELF_addminimal=') 3423read(20,*) c80tmp,ELF_addminimal 3424call loclabel(20,'ELFLOL_type=') 3425read(20,*) c80tmp,ELFLOL_type 3426call loclabel(20,'ELFLOL_cut=') 3427read(20,*) c80tmp,ELFLOL_cut 3428call loclabel(20,'iALIEdecomp=') 3429read(20,*) c80tmp,iALIEdecomp 3430call loclabel(20,'srcfuncmode=') 3431read(20,*) c80tmp,srcfuncmode 3432call loclabel(20,'atomdenscut=') 3433read(20,*) c80tmp,atomdenscut 3434call loclabel(20,'aug1D=') 3435read(20,*) c80tmp,aug1D 3436call loclabel(20,'aug2D=') 3437read(20,*) c80tmp,aug2D 3438call loclabel(20,'aug3D=') 3439read(20,*) c80tmp,aug3D 3440call loclabel(20,'num1Dpoints=') 3441read(20,*) c80tmp,num1Dpoints 3442call loclabel(20,'nprevorbgrid=') 3443read(20,*) c80tmp,nprevorbgrid 3444call loclabel(20,'bndordthres=') 3445read(20,*) c80tmp,bndordthres 3446call loclabel(20,'compthres=') 3447read(20,*) c80tmp,compthres 3448call loclabel(20,'compthresCDA=') 3449read(20,*) c80tmp,compthresCDA 3450call loclabel(20,'ispheratm=') 3451read(20,*) c80tmp,ispheratm 3452call loclabel(20,'laplfac=') 3453read(20,*) c80tmp,laplfac 3454call loclabel(20,'ipolarpara=') 3455read(20,*) c80tmp,ipolarpara 3456call loclabel(20,'ishowchgtrans=') 3457read(20,*) c80tmp,ishowchgtrans 3458call loclabel(20,'SpherIVgroup=') 3459read(20,*) c80tmp,SpherIVgroup 3460call loclabel(20,'MCvolmethod=') 3461read(20,*) c80tmp,MCvolmethod 3462call loclabel(20,'readEDF=') 3463read(20,*) c80tmp,readEDF 3464call loclabel(20,'isupplyEDF=') 3465read(20,*) c80tmp,isupplyEDF 3466call loclabel(20,'idelvirorb=') 3467read(20,*) c80tmp,idelvirorb 3468call loclabel(20,'ifchprog=') 3469read(20,*) c80tmp,ifchprog 3470call loclabel(20,'ishowptESP=') 3471read(20,*) c80tmp,ishowptESP 3472call loclabel(20,'imolsurparmode=') 3473read(20,*) c80tmp,imolsurparmode 3474call loclabel(20,'steric_addminimal=') 3475read(20,*) c80tmp,steric_addminimal 3476call loclabel(20,'steric_potcutrho=') 3477read(20,*) c80tmp,steric_potcutrho 3478call loclabel(20,'steric_potcons=') 3479read(20,*) c80tmp,steric_potcons 3480call loclabel(20,'NICSnptlim=') 3481read(20,*) c80tmp,NICSnptlim 3482call loclabel(20,'iplaneextdata=') 3483read(20,*) c80tmp,iplaneextdata 3484call loclabel(20,'igenP=') 3485read(20,*) c80tmp,igenP 3486call loclabel(20,'igenDbas=') 3487read(20,*) c80tmp,igenDbas 3488call loclabel(20,'igenMagbas=') 3489read(20,*) c80tmp,igenMagbas 3490call loclabel(20,'iloadasCart=') 3491read(20,*) c80tmp,iloadasCart 3492!Below are the parameters involved in plotting 3493call loclabel(20,'symbolsize=') 3494read(20,*) c80tmp,symbolsize 3495call loclabel(20,'pleatmlabsize=') 3496read(20,*) c80tmp,pleatmlabsize 3497call loclabel(20,'disshowlabel=') 3498read(20,*) c80tmp,disshowlabel 3499call loclabel(20,'iatom_on_contour_far=') 3500read(20,*) c80tmp,iatom_on_contour_far 3501call loclabel(20,'iatmlabtype=') 3502read(20,*) c80tmp,iatmlabtype 3503call loclabel(20,'iatmlabtype3D=') 3504read(20,*) c80tmp,iatmlabtype3D 3505call loclabel(20,'graphformat=') 3506read(20,*) c80tmp,graphformat 3507call loclabel(20,'graph1Dsize=') 3508read(20,*) c80tmp,graph1Dwidth,graph1Dheight 3509call loclabel(20,'graph2Dsize=') 3510read(20,*) c80tmp,graph2Dwidth,graph2Dheight 3511call loclabel(20,'graph3Dsize=') 3512read(20,*) c80tmp,graph3Dwidth,graph3Dheight 3513call loclabel(20,'numdigxyz=') 3514read(20,*) c80tmp,numdigx,numdigy,numdigz 3515call loclabel(20,'numdiglinexy=') 3516read(20,*) c80tmp,numdiglinex,numdigliney 3517call loclabel(20,'numdigctr=') 3518read(20,*) c80tmp,numdigctr 3519call loclabel(20,'fillcoloritpxy=') 3520read(20,*) c80tmp,fillcoloritpx,fillcoloritpy 3521call loclabel(20,'inowhiteblack=') 3522read(20,*) c80tmp,inowhiteblack 3523call loclabel(20,'isurfstyle=') 3524read(20,*) c80tmp,isurfstyle 3525call loclabel(20,'bondRGB=') 3526read(20,*) c80tmp,bondclrR,bondclrG,bondclrB 3527call loclabel(20,'atmlabRGB=') 3528read(20,*) c80tmp,atmlabclrR,atmlabclrG,atmlabclrB 3529call loclabel(20,'CP_RGB=') 3530read(20,*) c80tmp,CP3n3RGB,CP3n1RGB,CP3p1RGB,CP3p3RGB 3531call loclabel(20,'atmcolorfile=') !Set atom 3D color either according to external file or default setting 3532read(20,*) c80tmp,c200tmp 3533inquire(file=c200tmp,exist=alive) 3534if (index(c200tmp,"none")==0.and.alive) then 3535 write(*,"(' Note: Loading atom color settings from ',a)") trim(c200tmp) 3536 open(21,file=c200tmp,status="old") 3537 do iele=0,nelesupp 3538 read(21,*) inouse,atm3Dclr(iele,:) 3539 end do 3540 close(21) 3541end if 3542!Below are parameters about system 3543call loclabel(20,'nthreads=') 3544read(20,*) c80tmp,iniNThreads 3545call loclabel(20,'ompstacksize=') 3546read(20,*) c80tmp,ompstacksize 3547call loclabel(20,'gaupath=') 3548read(20,*) c80tmp,gaupath 3549call loclabel(20,'isilent=') 3550read(20,*) c80tmp,isilent 3551call loclabel(20,'isys=') 3552read(20,*) c80tmp,isys 3553call loclabel(20,'imodlayout=') 3554read(20,*) c80tmp,imodlayout 3555call loclabel(20,'outmedinfo=') 3556read(20,*) c80tmp,outmedinfo 3557call loclabel(20,'iwfntmptype=') 3558read(20,*) c80tmp,iwfntmptype 3559call loclabel(20,'ispecial=') 3560read(20,*) c80tmp,ispecial 3561!The last opened file name 3562call loclabel(20,'lastfile=') 3563read(20,"(10x,a)") lastfile 3564close(20) 3565end subroutine 3566