1!!-------- RDG analysis for MD, use promolecular approximation 2subroutine RDG_MD 3use defvar 4use util 5use function 6implicit real*8 (a-h,o-z) 7!The first index of avggrad and the first two indices of avghess correspond to components of gradient and Hessian, respectively 8real*8,allocatable :: avgdens(:,:,:),avggrad(:,:,:,:),avghess(:,:,:,:,:) 9real*8 denstmp,gradtmp(3),hesstmp(3,3),eigvecmat(3,3),eigval(3) 10real*8,allocatable :: avgRDG(:,:,:),thermflu(:,:,:),avgsl2r(:,:,:) !Final result 11real*8,allocatable :: scatterx(:),scattery(:) 12integer walltime1,walltime2 13real*8 time_end,time_endtmp 14 15write(*,*) "Input the range of the frames to be analyzed, e.g. 150,400" 16write(*,*) "Note: The index of frame starts from 1" 17read(*,*) ifpsstart,ifpsend 18call setgrid(0,igridsel) 19 20nfps=ifpsend-ifpsstart+1 21write(*,"(' Totally',i8,' frames, from frame',i8,' to',i8,' will be processed')") nfps,ifpsstart,ifpsend 22 23!Generate averaged density, averaged gradient and averaged Hessian 24allocate(avgdens(nx,ny,nz),avggrad(3,nx,ny,nz),avghess(3,3,nx,ny,nz)) 25avgdens=0D0 26avggrad=0D0 27avghess=0D0 28call walltime(walltime1) 29CALL CPU_TIME(time_begin) 30open(10,file=filename,access="sequential",status="old") 31write(*,*) 32write(*,*) "Now calculate averaged density, density gradient and density Hessian" 33do ifps=1,ifpsend 34 call readxyz(filename,1,0) 35 if (ifps<ifpsstart) cycle 36 write(*,"(' Processing frame',i8,' ...')") ifps 37 !fragatm must be defined for each frame, because each frame may have different content, and calchessmat_prodens will use it 38 nfragatmnum=ncenter 39 if (allocated(fragatm)) deallocate(fragatm) 40 allocate(fragatm(nfragatmnum)) 41 do itmp=1,nfragatmnum 42 fragatm(itmp)=itmp 43 end do 44nthreads=getNThreads() 45!$OMP PARALLEL DO SHARED(avgdens,avggrad,avghess) PRIVATE(i,j,k,tmpx,tmpy,tmpz,denstmp,gradtmp,hesstmp) schedule(dynamic) NUM_THREADS(nthreads) 46 do k=1,nz 47 tmpz=orgz+(k-1)*dz 48 do j=1,ny 49 tmpy=orgy+(j-1)*dy 50 do i=1,nx 51 tmpx=orgx+(i-1)*dx 52 call calchessmat_prodens(tmpx,tmpy,tmpz,denstmp,gradtmp,hesstmp) 53 avgdens(i,j,k)=avgdens(i,j,k)+denstmp 54 avggrad(:,i,j,k)=avggrad(:,i,j,k)+gradtmp 55 avghess(:,:,i,j,k)=avghess(:,:,i,j,k)+hesstmp 56 end do 57 end do 58 end do 59!$OMP END PARALLEL DO 60end do 61close(10) 62avgdens=avgdens/nfps 63avggrad=avggrad/nfps 64avghess=avghess/nfps 65 66!Use averaged density, averaged gradient and averaged Hessian to compute averaged RDG and averaged Sign(lambda2)*rho 67write(*,*) 68write(*,*) "Calculating averaged RDG and averaged Sign(lambda2)*rho..." 69allocate(avgRDG(nx,ny,nz),avgsl2r(nx,ny,nz)) 70do k=1,nz 71 tmpz=orgz+(k-1)*dz 72 do j=1,ny 73 tmpy=orgy+(j-1)*dy 74 do i=1,nx 75 tmpx=orgx+(i-1)*dx 76 77 avggradnormtmp=dsqrt(sum(avggrad(:,i,j,k)**2)) 78 if (RDGprodens_maxrho/=0D0.and.avgdens(i,j,k)>=RDGprodens_maxrho) then 79 avgRDG(i,j,k)=100D0 80 else if (avggradnormtmp==0D0.or.avgdens(i,j,k)==0D0) then 81 avgRDG(i,j,k)=999D0 82 else 83 avgRDG(i,j,k)=0.161620459673995D0*avggradnormtmp/avgdens(i,j,k)**(4D0/3D0) !0.161620459673995D0=1/(2*(3*pi**2)**(1/3)) 84 end if 85 86 call diagmat(avghess(:,:,i,j,k),eigvecmat,eigval,100,1D-6) 87 call sort(eigval) 88 if (eigval(2)/=0D0) then 89 avgsl2r(i,j,k)=avgdens(i,j,k)*eigval(2)/abs(eigval(2)) !At nuclei of single atom system, Hessian returned may be zero matrix 90 else 91 avgsl2r(i,j,k)=-avgdens(i,j,k) !Around nuclei, eigval(2)/abs(eigval(2)) always be negative 92 end if 93 94 end do 95 end do 96end do 97CALL CPU_TIME(time_end) 98call walltime(walltime2) 99write(*,"(' Calculation totally took up CPU time',f12.2,'s, wall clock time',i10,'s')") time_end-time_begin,walltime2-walltime1 100 101deallocate(avghess) !avghess will not be used further, so release its memory 102allocate(scatterx(nx*ny*nz),scattery(nx*ny*nz)) 103ii=1 104do k=1,nz 105 do j=1,ny 106 do i=1,nx 107 scatterx(ii)=avgsl2r(i,j,k) 108 scattery(ii)=avgRDG(i,j,k) 109 ii=ii+1 110 end do 111 end do 112end do 113!Default axis range of scatter plot 114xmin=-RDGprodens_maxrho 115xmax=RDGprodens_maxrho 116if (RDGprodens_maxrho==0.0D0) xmin=-2.0D0 117if (RDGprodens_maxrho==0.0D0) xmax=2.0D0 118ymin=0.0D0 119ymax=1.0D0 120 121write(*,*) 122do while (.true.) 123 write(*,*) "0 Return" 124 write(*,*) "1 Draw scatter graph between averaged RDG and Sign(lambda2)*rho" 125 write(*,*) "2 Save the scatter graph to file" 126 write(*,*) "3 Change range of X-axis of the scatter graph" 127 write(*,*) "4 Change range of Y-axis of the scatter graph" 128 write(*,*) "5 Export scatter points to output.txt" 129 write(*,*) "6 Export averaged RDG and Sign(lambda2)*rho to cube files in current folder" 130 write(*,*) "7 Compute and export thermal fluctuation index to cube file in current folder" 131 write(*,*) "8 Show isosurface of averaged RDG" 132 read(*,*) isel 133 134 if (isel==0) then 135 exit 136 else if (isel==1) then 137 write(*,*) "Drawing graph, please wait..." 138 else if (isel==2) then 139 isavepic=1 140 isavepic=0 141 write(*,"(a,a,a)") " Graph have been saved to ",trim(graphformat)," file with ""DISLIN"" prefix in current directory" 142 else if (isel==3) then 143 write(*,*) "Input lower limit and upper limit of X axis e.g. 0,1.5" 144 read(*,*) xmin,xmax 145 else if (isel==4) then 146 write(*,*) "Input lower limit and upper limit of Y axis e.g. 0,1.5" 147 read(*,*) ymin,ymax 148 else if (isel==5) then 149 open(10,file="output.txt",status="replace") 150 write(*,*) "Outputting output.txt in current folder..." 151 write(10,"(3f11.6,2E16.8)") ((( (orgx+(i-1)*dx),(orgy+(j-1)*dy),(orgz+(k-1)*dz),avgRDG(i,j,k),avgsl2r(i,j,k),i=1,nx),j=1,ny),k=1,nz) 152 close(10) 153 write(*,"(a)") " Finished, column 1/2/3/4/5 = X/Y/Z/averaged RDG/averaged Sign(lambda2)*rho, unit is Bohr" 154 else if (isel==6) then 155 write(*,*) "Outputting averaged reduced density gradient to avgRDG.cub in current folder" 156 open(10,file="avgRDG.cub",status="replace") 157 call outcube(avgRDG,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 158 close(10) 159 write(*,*) "Done!" 160 write(*,*) 161 write(*,*) "Outputting averaged Sign(lambda2)*rho to avgsl2r.cub in current folder" 162 open(10,file="avgsl2r.cub",status="replace") 163 call outcube(avgsl2r,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 164 close(10) 165 write(*,*) "Done!" 166 else if (isel==7) then 167 call walltime(walltime1) 168 CALL CPU_TIME(time_begin) 169 write(*,*) "Calculating thermal fluctuation index..." 170 if (allocated(thermflu)) deallocate(thermflu) 171 allocate(thermflu(nx,ny,nz)) 172 thermflu=0D0 173 open(10,file=filename,access="sequential",status="old") 174 !Calculate fluctuation square term of density and temporarily store it to thermflu array 175 do ifps=1,ifpsend 176 call readxyz(filename,1,0) 177 if (ifps<ifpsstart) cycle 178 write(*,"(' Processing frame',i8,' ...')") ifps 179nthreads=getNThreads() 180!$OMP PARALLEL DO SHARED(thermflu) PRIVATE(i,j,k,tmpx,tmpy,tmpz,denstmp) schedule(dynamic) NUM_THREADS(nthreads) 181 do k=1,nz 182 tmpz=orgz+(k-1)*dz 183 do j=1,ny 184 tmpy=orgy+(j-1)*dy 185 do i=1,nx 186 tmpx=orgx+(i-1)*dx 187 call calchessmat_prodens(tmpx,tmpy,tmpz,denstmp) 188 thermflu(i,j,k)=thermflu(i,j,k)+(denstmp-avgdens(i,j,k))**2 189 end do 190 end do 191 end do 192!$OMP END PARALLEL DO 193 end do 194 close(10) 195 thermflu=dsqrt(thermflu/nfps)/avgdens 196 CALL CPU_TIME(time_end) 197 call walltime(walltime2) 198 write(*,"(' Totally took up CPU time',f12.2,'s, wall clock time',i10,'s',/)") time_end-time_begin,walltime2-walltime1 199 write(*,*) "Outputting thermal fluctuation index to thermflu.cub in current folder" 200 open(10,file="thermflu.cub",status="replace") 201 call outcube(thermflu,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 202 close(10) 203 write(*,*) "Done!" 204 else if (isel==8) then 205 write(*,*) "Please wait..." 206 sur_value=0.3D0 207 if (allocated(cubmat)) deallocate(cubmat) 208 allocate(cubmat(nx,ny,nz)) 209 cubmat=avgRDG 210 deallocate(cubmat) 211 end if 212 write(*,*) 213end do 214 215end subroutine 216 217 218!!----- Calculate atomic and bond dipole moments in Hilbert space 219!For derivation, see Ideas of Quantum Chemistry, p634 220subroutine atmbonddip 221use defvar 222use util 223implicit real*8 (a-h,o-z) 224real*8 xdipmat(nbasis,nbasis),ydipmat(nbasis,nbasis),zdipmat(nbasis,nbasis),Ptottmp(nbasis,nbasis) 225character c80tmp*80 226!!!!Beware that the dipole moment integral has taken the negative sign of electron charge into account! 227if (igenDbas==0) then !Haven't calculated dipole moment integral matrix, so reload the input file and calculate it 228 Ptottmp=Ptot !Backup Ptot, which may have already been modified by users (via modifying occ), otherwise it will be flushed when loading 229 igenDbas=1 230 write(*,*) "Reloading input file and meantime generating dipole moment integral matrix..." 231 call dealloall 232 call readinfile(firstfilename,1) 233 Ptot=Ptottmp 234end if 235xdipmat=Dbas(1,:,:) 236ydipmat=Dbas(2,:,:) 237zdipmat=Dbas(3,:,:) 238 239! call showmatgau(xdipmat,"dip x",0,"f14.8",6) 240! call showmatgau(ydipmat,"dip y",0,"f14.8",6) 241! call showmatgau(zdipmat,"dip z",0,"f14.8",6) 242 243!Calculate total dipole moment 244xnucdip=sum(a(:)%charge*a(:)%x) 245ynucdip=sum(a(:)%charge*a(:)%y) 246znucdip=sum(a(:)%charge*a(:)%z) 247write(*,"(' Molecular nuclear dipole moment (a.u.):')") 248write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xnucdip,ynucdip,znucdip,dsqrt(xnucdip**2+ynucdip**2+znucdip**2) 249xeledip=sum(Ptot*xdipmat) 250yeledip=sum(Ptot*ydipmat) 251zeledip=sum(Ptot*zdipmat) 252write(*,"(' Molecular electron dipole moment (a.u.):')") 253write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xeledip,yeledip,zeledip,dsqrt(xeledip**2+yeledip**2+zeledip**2) 254xmoldip=xnucdip+xeledip 255ymoldip=ynucdip+yeledip 256zmoldip=znucdip+zeledip 257write(*,"(' Molecular dipole moment (a.u.):')") 258write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xmoldip,ymoldip,zmoldip,dsqrt(xmoldip**2+ymoldip**2+zmoldip**2) 259 260do while(.true.) 261 write(*,*) 262 write(*,*) " ----- Atomic and bond dipole moments in Hilbert space -----" 263 write(*,*) "0 Return" 264 write(*,*) "1 Output atomic dipole moment of specific atom" 265 write(*,*) "2 Output bond dipole moment of specific atom pair" 266 write(*,*) "3 Output atomic overall dipole moment of specific atom (Mulliken partition)" 267 write(*,*) "10 Export entire dipole moment matrix" 268 read(*,*) isel 269 if (isel==0) then 270 exit 271 else if (isel==1) then 272 do while(.true.) 273 write(*,*) "Input the atom index, e.g. 5" 274 write(*,*) "Note: Input 0 can return, input -1 can output result for all atoms" 275 read(*,*) isel2 276 if (isel2==0) then 277 exit 278 else if (isel2==-1) then 279 iatmstart=1 280 iatmend=ncenter 281 else 282 if (isel2>ncenter) then 283 write(*,*) "Atom index exceeded valid range! Input again" 284 cycle 285 end if 286 iatmstart=isel2 287 iatmend=isel2 288 end if 289 do iatm=iatmstart,iatmend 290 write(*,"(' Result of atom',i8,' (',a2,')')") iatm,a(iatm)%name 291 istart=basstart(iatm) 292 iend=basend(iatm) 293 atmelepop=sum(Ptot(istart:iend,istart:iend)*Sbas(istart:iend,istart:iend)) 294 xatmelediptot=sum(Ptot(istart:iend,istart:iend)*xdipmat(istart:iend,istart:iend)) 295 yatmelediptot=sum(Ptot(istart:iend,istart:iend)*ydipmat(istart:iend,istart:iend)) 296 zatmelediptot=sum(Ptot(istart:iend,istart:iend)*zdipmat(istart:iend,istart:iend)) 297 xatmeledip=xatmelediptot+atmelepop*a(iatm)%x 298 yatmeledip=yatmelediptot+atmelepop*a(iatm)%y 299 zatmeledip=zatmelediptot+atmelepop*a(iatm)%z 300 xatmnucdip=a(iatm)%charge*a(iatm)%x 301 yatmnucdip=a(iatm)%charge*a(iatm)%y 302 zatmnucdip=a(iatm)%charge*a(iatm)%z 303 xatmdiptot=xatmnucdip+xatmelediptot 304 yatmdiptot=yatmnucdip+yatmelediptot 305 zatmdiptot=zatmnucdip+zatmelediptot 306 write(*,"(' Atomic local population number:',f12.6)") atmelepop 307 write(*,"(' Atomic dipole moment (a.u.):')") 308 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xatmeledip,yatmeledip,zatmeledip,dsqrt(xatmeledip**2+yatmeledip**2+zatmeledip**2) 309 write(*,"(' Contribution to system dipole moment due to nuclear charge (a.u.):')") 310 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xatmnucdip,yatmnucdip,zatmnucdip,dsqrt(xatmnucdip**2+yatmnucdip**2+zatmnucdip**2) 311 write(*,"(' Contribution to system dipole moment due to electron (a.u.):')") 312 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xatmelediptot,yatmelediptot,zatmelediptot,dsqrt(xatmelediptot**2+yatmelediptot**2+zatmelediptot**2) 313 write(*,"(' Contribution to system dipole moment (a.u.):')") 314 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xatmdiptot,yatmdiptot,zatmdiptot,dsqrt(xatmdiptot**2+yatmdiptot**2+zatmdiptot**2) 315 write(*,*) 316 end do 317 end do 318 else if (isel==2) then 319 do while(.true.) 320 write(*,*) "Input the index of two atoms, e.g. 5,8" 321 write(*,*) "Note: Input q can return. Input b can output result for all bonds" 322 read(*,"(a)") c80tmp 323 if (index(c80tmp,'q')/=0) then 324 exit 325 else if (index(c80tmp,'b')/=0) then 326 write(*,*) "Notice that the bonds are determined according to distance criterion" 327 write(*,*) 328 bondcritval=1.15D0 329 else 330 read(c80tmp,*) iatomsel1,iatomsel2 331 if (iatomsel1>ncenter.or.iatomsel2>ncenter) then 332 write(*,*) "Atom index exceeded valid range! Input again" 333 cycle 334 end if 335 if (iatomsel1>iatomsel2) then 336 itmp=iatomsel2 337 iatomsel2=iatomsel1 338 iatomsel1=itmp 339 end if 340 end if 341 do iatm=1,ncenter 342 do jatm=iatm+1,ncenter 343 bonddist=dsqrt((a(iatm)%x-a(jatm)%x)**2+(a(iatm)%y-a(jatm)%y)**2+(a(iatm)%z-a(jatm)%z)**2) 344 if (index(c80tmp,'b')/=0) then 345 if (bonddist>( covr(a(iatm)%index)+covr(a(jatm)%index) )*bondcritval) cycle 346 else 347 if (iatm/=iatomsel1.or.jatm/=iatomsel2) cycle 348 end if 349 xcen=(a(iatm)%x+a(jatm)%x)/2D0 350 ycen=(a(iatm)%y+a(jatm)%y)/2D0 351 zcen=(a(iatm)%z+a(jatm)%z)/2D0 352 write(*,"(' Result between atom',i7,' (',a2,') and atom',i7,' (',a2,'), distance:',f10.5,' Ang')") iatm,a(iatm)%name,jatm,a(jatm)%name,bonddist*b2a 353 istart=basstart(iatm) 354 iend=basend(iatm) 355 jstart=basstart(jatm) 356 jend=basend(jatm) 357 bondpop=2*sum(Ptot(istart:iend,jstart:jend)*Sbas(istart:iend,jstart:jend)) !The matrix is symmetical, so multiplied by 2 358 xbonddiptot=2*sum(Ptot(istart:iend,jstart:jend)*xdipmat(istart:iend,jstart:jend)) 359 ybonddiptot=2*sum(Ptot(istart:iend,jstart:jend)*ydipmat(istart:iend,jstart:jend)) 360 zbonddiptot=2*sum(Ptot(istart:iend,jstart:jend)*zdipmat(istart:iend,jstart:jend)) 361 xbonddip=xbonddiptot+bondpop*xcen 362 ybonddip=ybonddiptot+bondpop*ycen 363 zbonddip=zbonddiptot+bondpop*zcen 364 write(*,"(' Bond population number (Overlap population):',f12.6)") bondpop 365 write(*,"(' Bond dipole moment (a.u.):')") 366 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xbonddip,ybonddip,zbonddip,dsqrt(xbonddip**2+ybonddip**2+zbonddip**2) 367 write(*,"(' Contribution to system dipole moment (a.u.):')") 368 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xbonddiptot,ybonddiptot,zbonddiptot,dsqrt(xbonddiptot**2+ybonddiptot**2+zbonddiptot**2) 369 write(*,*) 370 end do 371 end do 372 end do 373 else if (isel==3) then 374 do while(.true.) 375 write(*,*) "Input the atom index, e.g. 5" 376 write(*,*) "Note: Input 0 can return, input -1 can output result for all atoms" 377 read(*,*) isel2 378 if (isel2==0) then 379 exit 380 else if (isel2==-1) then 381 iatmstart=1 382 iatmend=ncenter 383 else 384 if (isel2>ncenter) then 385 write(*,*) "Atom index exceeded valid range! Input again" 386 cycle 387 end if 388 iatmstart=isel2 389 iatmend=isel2 390 end if 391 do iatm=iatmstart,iatmend 392 write(*,"(' Result of atom',i8,' (',a2,')')") iatm,a(iatm)%name 393 atmelepop=0D0 394 xatmelediptot=0D0 395 yatmelediptot=0D0 396 zatmelediptot=0D0 397 istart=basstart(iatm) 398 iend=basend(iatm) 399 do jatm=1,ncenter 400 jstart=basstart(jatm) 401 jend=basend(jatm) 402 atmelepop=atmelepop+sum(Ptot(istart:iend,jstart:jend)*Sbas(istart:iend,jstart:jend)) 403 xatmelediptot=xatmelediptot+sum(Ptot(istart:iend,jstart:jend)*xdipmat(istart:iend,jstart:jend)) 404 yatmelediptot=yatmelediptot+sum(Ptot(istart:iend,jstart:jend)*ydipmat(istart:iend,jstart:jend)) 405 zatmelediptot=zatmelediptot+sum(Ptot(istart:iend,jstart:jend)*zdipmat(istart:iend,jstart:jend)) 406 end do 407 xatmeledip=xatmelediptot+atmelepop*a(iatm)%x 408 yatmeledip=yatmelediptot+atmelepop*a(iatm)%y 409 zatmeledip=zatmelediptot+atmelepop*a(iatm)%z 410 xatmnucdip=a(iatm)%charge*a(iatm)%x 411 yatmnucdip=a(iatm)%charge*a(iatm)%y 412 zatmnucdip=a(iatm)%charge*a(iatm)%z 413 xatmdiptot=xatmnucdip+xatmelediptot 414 yatmdiptot=yatmnucdip+yatmelediptot 415 zatmdiptot=zatmnucdip+zatmelediptot 416 write(*,"(' Atomic Mulliken population number:',f12.6)") atmelepop 417 write(*,"(' Atomic overall dipole moment (a.u.):')") 418 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xatmeledip,yatmeledip,zatmeledip,dsqrt(xatmeledip**2+yatmeledip**2+zatmeledip**2) 419 write(*,"(' Contribution to system dipole moment due to nuclear charge (a.u.):')") 420 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xatmnucdip,yatmnucdip,zatmnucdip,dsqrt(xatmnucdip**2+yatmnucdip**2+zatmnucdip**2) 421 write(*,"(' Contribution to system dipole moment due to electron (a.u.):')") 422 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xatmelediptot,yatmelediptot,zatmelediptot,dsqrt(xatmelediptot**2+yatmelediptot**2+zatmelediptot**2) 423 write(*,"(' Contribution to system dipole moment (a.u.):')") 424 write(*,"(' X=',f12.6,' Y=',f12.6,' Z=',f12.6,' Norm=',f12.6)") xatmdiptot,yatmdiptot,zatmdiptot,dsqrt(xatmdiptot**2+yatmdiptot**2+zatmdiptot**2) 425 write(*,*) 426 end do 427 end do 428 else if (isel==10) then 429 open(10,file="dipmatx.txt",status="replace") 430 call showmatgau(Ptot*xdipmat,"",1,"f14.8",10) 431 close(10) 432 open(10,file="dipmaty.txt",status="replace") 433 call showmatgau(Ptot*ydipmat,"",1,"f14.8",10) 434 close(10) 435 open(10,file="dipmatz.txt",status="replace") 436 call showmatgau(Ptot*zdipmat,"",1,"f14.8",10) 437 close(10) 438 write(*,"(a)") " X, Y and Z components of electron dipole moment matrix have been outputted to dipmatx, dipmaty and dipmatz.txt in current folder, respectively" 439 end if 440end do 441end subroutine 442 443 444!!------------ Generate cube file for multiple orbitals 445subroutine genmultiorbcube 446use defvar 447use util 448implicit real*8 (a-h,o-z) 449integer orbsellist(nmo) 450integer tmparr(nmo+1) 451character c1000tmp*1000,cubname*20 452real*8,allocatable :: orbcubmat(:,:,:,:) 453write(*,"(a)") " Input orbital index. e.g. 1,3-6,8,10-11 denotes 1,3,4,5,6,8,10,11" 454write(*,*) "Input q can return" 455read(*,"(a)") c1000tmp 456if (index(c1000tmp,'q')/=0) return 457call str2arr(c1000tmp,norbsel,orbsellist) 458if ( any(orbsellist(1:norbsel)<1) .or. any(orbsellist(1:norbsel)>nmo) ) then 459 write(*,*) "Error: The orbitals you selected exceeded valid range!" 460 return 461end if 462call setgrid(0,itmp) 463if (allocated(cubmat)) deallocate(cubmat) 464allocate(cubmat(nx,ny,nz)) 465write(*,*) 466write(*,*) "1 Output the grid data of these orbitals as separate cube files" 467write(*,*) "2 Output the grid data of these orbitals as a single cube file" 468read(*,*) ioutmode 469 470if (ioutmode==1) then 471 do iorbidx=1,norbsel 472 iorb=orbsellist(iorbidx) 473 write(cubname,"('orb',i6.6,'.cub')") iorb 474 write(*,"(' Calculating and exporting orbital',i6)") iorb 475 call savecubmat(4,1,iorb) 476 open(10,file=cubname,status="replace") 477 call outcube(cubmat,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 478 close(10) 479 write(*,"(' Orbital',i7,' has been exported to ',a,' in current folder',/)") iorb,trim(cubname) 480 end do 481 482else if (ioutmode==2) then 483 allocate(orbcubmat(nx,ny,nz,norbsel)) 484 do iorbidx=1,norbsel 485 iorb=orbsellist(iorbidx) 486 write(*,"(a,i6,a)") " Calculating grid data for orbital",iorb,"..." 487 call savecubmat(4,1,iorb) 488 orbcubmat(:,:,:,iorbidx)=cubmat(:,:,:) 489 end do 490 where (abs(orbcubmat)<=1D-99) orbcubmat=0D0 !Diminish too small value, otherwise the symbol "E" cannot be shown by 1PE13.5 format e.g. 9.39376-116, 491 write(*,*) 492 write(*,*) "Exporting cube file, please wait..." 493 open(10,file="orbital.cub",status="replace") 494 write(10,"(' Generated by Multiwfn')") 495 write(10,"(' Totally ',i12,' grid points')") nx*ny*nz 496 write(10,"(i5,3f12.6)") -ncenter,orgx,orgy,orgz 497 write(10,"(i5,3f12.6)") nx,dx,0.0,0.0 498 write(10,"(i5,3f12.6)") ny,0.0,dy,0.0 499 write(10,"(i5,3f12.6)") nz,0.0,0.0,dz 500 do i=1,ncenter 501 write(10,"(i5,4f12.6)") a(i)%index,a(i)%charge,a(i)%x,a(i)%y,a(i)%z 502 end do 503 tmparr(1)=norbsel 504 tmparr(2:norbsel+1)=orbsellist(1:norbsel) 505 write(10,"(10i5)") tmparr(1:norbsel+1) 506 do ix=1,nx 507 do iy=1,ny 508 write(10,"(6(1PE13.5))",advance="no") ((orbcubmat(ix,iy,iz,iorbidx),iorbidx=1,norbsel),iz=1,nz) 509 write(10,*) 510 end do 511 end do 512 close(10) 513 write(*,*) "The grid data of the orbitals have been stored to orbital.cub in current folder" 514 deallocate(orbcubmat) 515end if 516 517deallocate(cubmat) 518end subroutine 519 520 521 522 523!!---------- Generate grid data of iso-chemical shielding surfaces (ICSS) or related quantities 524subroutine ICSS 525use defvar 526use util 527implicit real*8 (a-h,o-z) 528character c200tmp*200,gauinpfile*200,gauoutfile*200,selectyn,suffix*4 529character,allocatable :: gauinpcontent(:)*79 530!Set grid for calculating NICS 531call setgrid(0,itmp) 532numbqper=NICSnptlim-ncenter 533write(*,"(' The number of Bq per batch:',i10)") numbqper 534npttot=nx*ny*nz 535nfile=ceiling(dfloat(npttot)/numbqper) 536!Generate Gaussian input file 537write(*,*) 538write(*,*) "If skip generating Gaussian input file of NMR task? (y/n)" 539read(*,*) selectyn 540if (selectyn=='n'.or.selectyn=='N') then 541 write(*,*) "Input the path of template Gaussian input file, e.g. c:\ltwd.gjf" 542 do while(.true.) 543 read(*,"(a)") gauinpfile 544 inquire(file=gauinpfile,exist=alive) 545 if (alive) exit 546 write(*,*) "Cannot find corresponding files, input again" 547 end do 548 open(10,file=gauinpfile,status="old") 549 numgauline=totlinenum(10,2) 550 allocate(gauinpcontent(numgauline)) 551 numblank=0 552 iendcoord=numgauline !Which line is the last line recording coordinates 553 do i=1,numgauline 554 read(10,"(a)") gauinpcontent(i) 555 if (gauinpcontent(i)==" ") then 556 numblank=numblank+1 557 if (numblank==3) iendcoord=i-1 558 end if 559 if (index(gauinpcontent(i),'#')/=0) then 560 gauinpcontent(i)=trim(gauinpcontent(i))//" NMR" 561 if (index(gauinpcontent(i),'conn')==0) gauinpcontent(i)=trim(gauinpcontent(i))//" geom=connectivity" 562 end if 563! write(*,"(a)") gauinpcontent(i) 564 end do 565 close(10) 566 567 gauinpfile="NICS" 568 do ifile=1,nfile 569 write(c200tmp,"(a,i4.4,a)") trim(gauinpfile),ifile,".gjf" 570 open(10,file=c200tmp,status="replace") 571 write(*,"(a,a,a)") " Outputting ",trim(c200tmp)," to current folder..." 572 !Write head part 573 do i=1,iendcoord 574 if (ifile>1.and.index(gauinpcontent(i),'#')/=0) then 575 write(10,"(a)") trim(gauinpcontent(i))//" guess=read" 576 else 577 write(10,"(a)") gauinpcontent(i) 578 end if 579 end do 580 !Write Bq information 581 itmp=0 582 do i=1,nx 583 do j=1,ny 584 do k=1,nz 585 itmp=itmp+1 586 rnowx=orgx+(i-1)*dx 587 rnowy=orgy+(j-1)*dy 588 rnowz=orgz+(k-1)*dz 589 if (itmp<=(ifile-1)*numbqper) cycle 590 if (itmp>ifile*numbqper) exit 591 write(10,"('Bq ',3f12.6)") rnowx*b2a,rnowy*b2a,rnowz*b2a 592 end do 593 end do 594 end do 595 write(10,*) 596 !Write connectivity explicitly 597 if (ifile/=nfile) then 598 do i=1,NICSnptlim 599 write(10,"(i7)") i 600 end do 601 else if (ifile==nfile) then 602 do i=1,npttot-(ifile-1)*numbqper+ncenter 603 write(10,"(i7)") i 604 end do 605 end if 606 !Write remaining part 607 do i=iendcoord+1,numgauline 608 write(10,"(a)") gauinpcontent(i) 609 end do 610 close(10) 611 end do 612end if 613write(*,*) "Now please run these input files by Gaussian" 614write(*,*) 615!Load NICS from Gaussian output file 616if (allocated(cubmat)) deallocate(cubmat) 617allocate(cubmat(nx,ny,nz)) 618write(*,*) "Input the path of Gaussian output file of NMR task" 619write(*,"(a)") " Assume that you input ""c:\ltwd\NICS"", then c:\ltwd\NICS0001.out, c:\ltwd\NICS0002.out, c:\ltwd\NICS0003.out... will be loaded" 620suffix=".out" 621do while(.true.) 622 read(*,"(a)") gauoutfile 623 inquire(file=trim(gauoutfile)//"0001"//suffix,exist=alive) 624 if (alive) exit 625 if (.not.alive) then 626 suffix=".log" 627 inquire(file=trim(gauoutfile)//"0001"//suffix,exist=alive) 628 end if 629 if (alive) exit 630 write(*,*) "Cannot find corresponding files, input again" 631end do 632100 write(*,*) "Load which term?" 633write(*,*) "1: Isotropic 2: Anisotropy 3: XX component 4: YY component 5: ZZ component" 634read(*,*) iload 635do ifile=1,nfile 636 write(c200tmp,"(a,i4.4,a)") trim(gauoutfile),ifile,suffix 637 open(10,file=c200tmp,status="old") 638 write(*,"(' Loading ',a,'...')") trim(c200tmp) 639 call loclabel(10,"GIAO Magnetic shielding tensor") 640 read(10,*) 641 !Detect format. The NMR output format changes since G09 D.01 to leave more space for atomic index 642 read(10,"(a80)") c200tmp 643 backspace(10) 644 iformat=1 645 if (c200tmp(25:25)=='=') iformat=2 !Since G09 D.01 646 647 do i=1,ncenter !Skip atom's result 648 read(10,*) 649 read(10,*) 650 read(10,*) 651 read(10,*) 652 read(10,*) 653 end do 654 itmp=0 655 do i=1,nx 656 do j=1,ny 657 do k=1,nz 658 itmp=itmp+1 659 if (itmp<=(ifile-1)*numbqper) cycle 660 if (itmp>ifile*numbqper) exit 661 if (iload==1.or.iload==2) then 662! read(10,"(a80)") c200tmp 663! write(*,"(a80)") c200tmp 664! backspace(10) 665 if (iformat==1) then 666 if (iload==1) read(10,"(22x,f10.4)") cubmat(i,j,k) 667 if (iload==2) read(10,"(48x,f10.4)") cubmat(i,j,k) 668 else if (iformat==2) then 669 if (iload==1) read(10,"(26x,f10.4)") cubmat(i,j,k) 670 if (iload==2) read(10,"(52x,f10.4)") cubmat(i,j,k) 671 end if 672 read(10,*) 673 read(10,*) 674 read(10,*) 675 read(10,*) 676 else 677 read(10,*) 678 if (iload==3) then 679 read(10,"(8x,f10.4)") cubmat(i,j,k) 680 read(10,*) 681 read(10,*) 682 else if (iload==4) then 683 read(10,*) 684 read(10,"(24x,f10.4)") cubmat(i,j,k) 685 read(10,*) 686 else if (iload==5) then 687 read(10,*) 688 read(10,*) 689 read(10,"(42x,f10.4)") cubmat(i,j,k) 690 end if 691 read(10,*) 692 end if 693 end do 694 end do 695 end do 696 close(10) 697end do 698write(*,*) "Loading finished!" 699do while(.true.) 700 write(*,*) 701 write(*,*) "-1 Load another property" 702 write(*,*) "0 Return to main menu" 703 if (iload==1) write(*,*) "1 Visualize iso-chemical shielding surface" 704 if (iload==2) write(*,*) "1 Visualize aniso-chemical shielding surface" 705 if (iload==3) write(*,*) "1 Visualize XX-chemical shielding surface" 706 if (iload==4) write(*,*) "1 Visualize YY-chemical shielding surface" 707 if (iload==5) write(*,*) "1 Visualize ZZ-chemical shielding surface" 708 if (iload==1) write(*,*) "2 Export the grid data to ICSS.cub current folder" 709 if (iload==2) write(*,*) "2 Export the grid data to AICSS.cub current folder" 710 if (iload==3) write(*,*) "2 Export the grid data to ICSSXX.cub current folder" 711 if (iload==4) write(*,*) "2 Export the grid data to ICSSYY.cub current folder" 712 if (iload==5) write(*,*) "2 Export the grid data to ICSSZZ.cub current folder" 713 read(*,*) isel 714 if (isel==-1) then 715 goto 100 716 else if (isel==0) then 717 exit 718 else if (isel==1) then 719 else if (isel==2) then 720 if (iload==1) open(10,file="ICSS.cub",status="replace") 721 if (iload==2) open(10,file="AICSS.cub",status="replace") 722 if (iload==3) open(10,file="ICSSXX.cub",status="replace") 723 if (iload==4) open(10,file="ICSSYY.cub",status="replace") 724 if (iload==5) open(10,file="ICSSZZ.cub",status="replace") 725 call outcube(cubmat,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 726 close(10) 727 write(*,"(a)") " The cube file has been exported to current folder" 728 end if 729end do 730end subroutine 731 732 733 734!!----- Plot radial distribution function for a real space function 735subroutine plotraddis 736use defvar 737use function 738implicit real*8 (a-h,o-z) 739real*8,allocatable :: potx(:),poty(:),potz(:),potw(:),radval(:),radpos(:),intradval(:),sphavgval(:) 740ifunc=1 741cenx=0D0 742ceny=0D0 743cenz=0D0 744rlow=0D0 745rhigh=5D0/b2a !5 Angstrom 746nsphpt=2030 747nradpt=500 748do while(.true.) 749 write(*,*) 750 write(*,*) " ====== Plot radial distribution function for a real space function ======" 751 write(*,*) "-1 Exit" 752 write(*,*) "0 Calculate radial distribution function and its integration curve" 753 write(*,"(a,i4)") " 1 Select real space function, current:",ifunc 754 write(*,"(a,3f10.4,' Ang')") " 2 Set sphere center, current",cenx*b2a,ceny*b2a,cenz*b2a 755 write(*,"(a,2f9.4,' Ang')") " 3 Set lower and upper limit of radial distance, current:",rlow*b2a,rhigh*b2a 756 write(*,"(a,i6)") " 4 Set the number of integration point in each shell, current:",nsphpt 757 write(*,"(a,i6)") " 5 Set the number of radial points, current:",nradpt 758 read(*,*) isel 759 if (isel==-1) then 760 return 761 else if (isel==1) then 762 call selfunc_interface(ifunc) 763 else if (isel==2) then 764 write(*,*) "Input sphere center (Angstrom), e.g. 0.0,1.2,-0.4" 765 read(*,*) cenx,ceny,cenz 766 cenx=cenx/b2a !Convert to Bohr 767 ceny=ceny/b2a 768 cenz=cenz/b2a 769 else if (isel==3) then 770 write(*,*) "Input lower and upper limit (Angstrom), e.g. 0.0,8.0" 771 read(*,*) rlow,rhigh 772 rlow=rlow/b2a !Convert to Bohr 773 rhigh=rhigh/b2a 774 else if (isel==4) then 775 write(*,"(a)") " Input the number of integration point in each shell, the value must be one of & 776 110/170/230/266/302/434/590/770/974/1454/1730/2030/2354/2702/3074/3470/3890/4334/4802/5294/5810" 777 read(*,*) nsphpt 778 else if (isel==5) then 779 write(*,*) "Input the number of radial points, e.g. 800" 780 read(*,*) nradpt 781 else if (isel==0) then 782 allocate(potx(nsphpt),poty(nsphpt),potz(nsphpt),potw(nsphpt)) 783 allocate(radval(nradpt),radpos(nradpt),intradval(nradpt),sphavgval(nradpt)) !radval records RDF, radpot records r position, intradval records integration of RDF 784 call Lebedevgen(nsphpt,potx,poty,potz,potw) 785 radval=0D0 786 radstp=(rhigh-rlow)/(nradpt-1) 787 ifinish=0 788 iprogstp=20 789 iprogcrit=iprogstp 790 write(*,*) "Calculating..." 791nthreads=getNThreads() 792!$OMP PARALLEL DO SHARED(radval,radpos,ifinish,iprogcrit) PRIVATE(irad,radnow,isph,rnowx,rnowy,rnowz,tmpval) schedule(dynamic) NUM_THREADS(nthreads) 793 do irad=1,nradpt 794 radnow=rlow+(irad-1)*radstp 795 radpos(irad)=radnow 796 tmpval=0 797 do isph=1,nsphpt 798 rnowx=potx(isph)*radnow+cenx 799 rnowy=poty(isph)*radnow+ceny 800 rnowz=potz(isph)*radnow+cenz 801 tmpval=tmpval+calcfuncall(ifunc,rnowx,rnowy,rnowz)*potw(isph) 802 end do 803 radval(irad)=4*pi*tmpval*radnow**2 !Multiplied by 4*pi is because the Lebedev integration routine produces unity rather than 4*pi 804 sphavgval(irad)=tmpval !Spherically average function 805 ifinish=ifinish+1 806 if (ifinish==iprogcrit) then 807 write(*,"(' Finished:',i6,' /',i6)") ifinish,nradpt 808 iprogcrit=iprogcrit+iprogstp 809 end if 810 end do 811!$OMP END PARALLEL DO 812 813 !Calculate integration of RDF 814 intradval(1)=0D0 815 do irad=2,nradpt 816 intradval(irad)=intradval(irad-1)+radval(irad-1)*radstp 817 end do 818 write(*,"(a,f22.10)") " Integrating the RDF in the specified range is",intradval(nradpt) 819 valrange=maxval(radval)-minval(radval) 820 valrangeint=maxval(intradval)-minval(intradval) 821 ilenunit1D=2 822 do while(.true.) 823 write(*,*) 824 if (ilenunit1D==1) write(*,*) "-1 Switch the length unit for plotting, current: Bohr" 825 if (ilenunit1D==2) write(*,*) "-1 Switch the length unit for plotting, current: Angstrom" 826 write(*,*) "0 Return" 827 write(*,*) "1 Plot the radial distribution function" 828 write(*,*) "2 Plot integration curve of the RDF" 829 write(*,*) "3 Save the radial distribution function map in current folder" 830 write(*,*) "4 Save integration curve of the RDF map in current folder" 831 write(*,*) "5 Export the radial distribution function to RDF.txt in current folder" 832 write(*,*) "6 Export integration curve of the RDF to intRDF.txt in current folder" 833 write(*,*) "7 Export the spherically averaged function" 834 read(*,*) isel 835 if (isel==-1) then 836 if (ilenunit1D==1) then 837 ilenunit1D=2 838 else if (ilenunit1D==2) then 839 ilenunit1D=1 840 end if 841 else if (isel==0) then 842 deallocate(potx,poty,potz,potw,radval,radpos,intradval,sphavgval) 843 exit 844 else if (isel==1.or.isel==3) then 845 ylow=minval(radval)-0.1D0*valrange 846 yhigh=maxval(radval)+0.1D0*valrange 847 if (isel==1) then 848 else 849 write(*,"(a,a,a)") " Graph have been saved to ",trim(graphformat)," file with ""DISLIN"" prefix in current directory" 850 end if 851 else if (isel==2.or.isel==4) then 852 ylow=minval(intradval)-0.1D0*valrangeint 853 yhigh=maxval(intradval)+0.1D0*valrangeint 854 if (isel==2) then 855 else if (isel==4) then 856 write(*,"(a,a,a)") " Graph have been saved to ",trim(graphformat)," file with ""DISLIN"" prefix in current directory" 857 end if 858 else if (isel==5) then 859 open(10,file="RDF.txt",status="replace") 860 do irad=1,nradpt 861 write(10,"(i7,f12.4,f22.10)") irad,radpos(irad)*b2a,radval(irad) 862 end do 863 close(10) 864 write(*,*) "The result has been output to RDF.txt in current folder" 865 write(*,*) "The second column is radial distance (Angstrom), the third column is value" 866 else if (isel==6) then 867 open(10,file="intRDF.txt",status="replace") 868 do irad=1,nradpt 869 write(10,"(i7,f12.4,f22.10)") irad,radpos(irad)*b2a,intradval(irad) 870 end do 871 close(10) 872 write(*,*) "The result has been output to intRDF.txt in current folder" 873 write(*,*) "The second column is radial distance (Angstrom), the third column is value" 874 else if (isel==7) then 875 open(10,file="sphavgval.txt",status="replace") 876 do irad=1,nradpt 877 write(10,"(i7,f12.6,f18.7)") irad,radpos(irad),sphavgval(irad) 878 end do 879 close(10) 880 write(*,*) "The result has been output to sphavgval.txt in current folder" 881 write(*,*) "The second column is radial distance (Bohr), the third column is value" 882 end if 883 end do 884 885 end if 886end do 887end subroutine 888 889 890 891!!--------- Analyze correspondence between orbitals in two wavefunctions 892subroutine orbcorres 893use defvar 894use function 895use util 896implicit real*8 (a-h,o-z) 897real*8,allocatable :: convmat(:,:) !(i,j) is the coefficient of j MO of the second wavefunction in i MO of current wavefunction 898real*8,allocatable :: MOvalgrd(:,:),MOvalgrd2(:,:) !MOvalgrd(j,n),MOvalgrd2(j,n) means the the value of the nth MO of the first/second wavefunction at the jth grid 899real*8,allocatable :: comparr(:),beckeweigrid(:) 900integer,allocatable :: comparridx(:) 901character filename2*200,c80tmp*80 902type(content),allocatable :: gridatm(:),gridatmorg(:) 903if (iautointgrid==1) then 904 sphpotold=sphpot 905 radpotold=radpot 906 sphpot=230 907 radpot=50 908end if 909allocate(gridatm(radpot*sphpot),gridatmorg(radpot*sphpot),beckeweigrid(radpot*sphpot)) 910 911do isep=nmo,1,-1 912 if (MOtype(isep)==1) exit 913end do 914if (wfntype==1.or.wfntype==4) write(*,"(' Note: The orbitals from',i6,' to',i6,' are alpha; from',i6,' to',i6,' are beta')") 1,isep,isep+1,nmo 915write(*,"(a)") " Input the range of the orbitals of present wavefunction to be considered, e.g. 2,9. & 916If press ENTER directly, all orbitals will be taken into account" 917read(*,"(a)") c80tmp 918if (c80tmp==" ") then 919 istart1=1 920 iend1=nmo 921else 922 read(c80tmp,*) istart1,iend1 923end if 924 925write(*,*) 926write(*,*) "Input filename of the second wavefunction, e.g. c:\ltwd.fch" 927do while(.true.) 928 read(*,"(a)") filename2 929 inquire(file=filename2,exist=alive) 930 if (alive) exit 931 write(*,*) "Cannot find the file, input again" 932end do 933call dealloall 934call readinfile(filename2,1) !Get some knowledge about the second wavefunction 935nmo2=nmo !The number of MOs in the wfn2 936iwfntype2=wfntype 937do isep=nmo2,1,-1 938 if (MOtype(isep)==1) exit 939end do 940write(*,*) 941if (iwfntype2==1.or.iwfntype2==4) write(*,"(' Note: The orbitals from',i6,' to',i6,' are alpha; from',i6,' to',i6,' are beta')") 1,isep,isep+1,nmo 942write(*,"(a)") " Input the range of the orbitals of the second wavefunction to be considered, e.g. 2,9. & 943If press ENTER directly, all orbitals will be taken into account" 944read(*,"(a)") c80tmp 945if (c80tmp==" ") then 946 istart2=1 947 iend2=nmo2 948else 949 read(c80tmp,*) istart2,iend2 950end if 951 952call dealloall 953call readinfile(firstfilename,1) 954allocate(MOvalgrd(radpot*sphpot,nmo),MOvalgrd2(radpot*sphpot,nmo2),convmat(nmo,nmo2)) 955convmat=0D0 956 957write(*,"(' Radial points:',i5,' Angular points:',i5,' Total:',i10,' per center')") radpot,sphpot,radpot*sphpot 958write(*,*) "Calculating, please wait..." 959call gen1cintgrid(gridatmorg,iradcut) 960 961call walltime(iwalltime1) 962CALL CPU_TIME(time_begin) 963 964do iatm=1,ncenter 965 write(*,"(' Progress: ',i5,' /',i5)") iatm,ncenter 966 gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule 967 gridatm%y=gridatmorg%y+a(iatm)%y 968 gridatm%z=gridatmorg%z+a(iatm)%z 969 970 !Calculate value of all MOs of the first and second wavefunction at all grids 971 call dealloall 972 call readinfile(filename2,1) !Load wfn2 973nthreads=getNThreads() 974!$OMP parallel do shared(MOvalgrd2) private(ipt) num_threads(nthreads) schedule(dynamic) 975 do ipt=1+iradcut*sphpot,radpot*sphpot 976 call orbderv(1,istart2,iend2,gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z,MOvalgrd2(ipt,:)) 977 end do 978!$OMP end parallel do 979 call dealloall 980 call readinfile(firstfilename,1) !Retrieve to wfn1 981nthreads=getNThreads() 982!$OMP parallel do shared(MOvalgrd) private(ipt) num_threads(nthreads) schedule(dynamic) 983 do ipt=1+iradcut*sphpot,radpot*sphpot 984 call orbderv(1,istart1,iend1,gridatm(ipt)%x,gridatm(ipt)%y,gridatm(ipt)%z,MOvalgrd(ipt,:)) 985 end do 986!$OMP end parallel do 987 988 !Calculate Becke weight at all grids 989 call gen1cbeckewei(iatm,iradcut,gridatm,beckeweigrid) 990 991nthreads=getNThreads() 992!$OMP parallel do shared(convmat) private(imo,jmo,tmpval,ipt) num_threads(nthreads) schedule(dynamic) 993 do imo=istart1,iend1 994 do jmo=istart2,iend2 995 tmpval=0D0 996 do ipt=1+iradcut*sphpot,radpot*sphpot 997 tmpval=tmpval+beckeweigrid(ipt)*gridatmorg(ipt)%value*MOvalgrd(ipt,imo)*MOvalgrd2(ipt,jmo) 998 end do 999 convmat(imo,jmo)=convmat(imo,jmo)+tmpval 1000 end do 1001 end do 1002!$OMP end parallel do 1003end do 1004CALL CPU_TIME(time_end) 1005call walltime(iwalltime2) 1006write(*,"(' Calculation took up CPU time',f12.2,'s, wall clock time',i10,'s',/)") time_end-time_begin,iwalltime2-iwalltime1 1007 1008! call showmatgau(convmat,"convmat",0,"f12.3") 1009 1010devmax=0D0 1011idevmax=1 1012! write(*,*) "The sum of composition of each orbital of current wavefunction" 1013if ((wfntype==0.or.wfntype==2).and.(iwfntype2==0.or.iwfntype2==2)) then !Both of the two wavefunctions are R or RO types 1014 do imo=istart1,iend1 !Check normalization 1015 totcomp=sum(convmat(imo,:)**2)*100D0 1016 ! write(*,"(i6,':',f12.5,' %')") imo,totcomp 1017 if (abs(totcomp-100D0)>devmax) then 1018 devmax=abs(totcomp-100D0) 1019 idevmax=imo 1020 end if 1021 end do 1022end if 1023write(*,"(' The maximum deviation to normalization condition is',f8.3,' % (Orbital',i6,')')") devmax,idevmax 1024write(*,"(a)") " Note: The first column below is the index of the orbitals in present wavefunction, the largest five contributions from & 1025the orbitals in the second wavefunction are shown at right side. If the dominative index is inconsistent to the first column, the row will be marked by asterisk" 1026write(*,*) 1027allocate(comparr(nmo2),comparridx(nmo2)) 1028do imo=istart1,iend1 1029 !Sort the composition array from small to large 1030 comparr(:)=convmat(imo,:) 1031 do itmp=1,nmo2 1032 comparridx(itmp)=itmp 1033 end do 1034 do i=istart2,iend2 1035 do j=i+1,iend2 1036 if (abs(comparr(i))>abs(comparr(j))) then 1037 temp=comparr(i) 1038 comparr(i)=comparr(j) 1039 comparr(j)=temp 1040 itemp=comparridx(i) 1041 comparridx(i)=comparridx(j) 1042 comparridx(j)=itemp 1043 end if 1044 end do 1045 end do 1046 if (comparridx(iend2)/=imo) then 1047 write(*,"('*',i5,': ')",advance="no") imo 1048 else 1049 write(*,"(' ',i5,': ')",advance="no") imo 1050 end if 1051 do i=iend2,iend2-4,-1 1052 write(*,"(i5,'(',f6.2,'%)')",advance="no") comparridx(i),comparr(i)**2*100D0 1053 end do 1054 write(*,*) 1055end do 1056 1057write(*,*) 1058do while(.true.) 1059 write(*,*) "Input the orbital index to print detail compositions and coefficients, e.g. 5" 1060 write(*,*) "input -1 can output all overlap integrals between the chosen orbitals" 1061 write(*,*) "Input 0 can exit" 1062 read(*,*) imo 1063 if (imo==0) then 1064 exit 1065 else if (imo==-1) then 1066 open(10,file="convmat.txt",status="replace") 1067 do i=istart1,iend1 1068 do j=istart2,iend2 1069 write(10,"(2i7,f12.6)") i,j,convmat(i,j) 1070 end do 1071 end do 1072 close(10) 1073 write(*,*) "The overlap integrals have been outputted to convmat.txt in current folder" 1074 write(*,"(a,/)") " The first and second columns correspond to the orbital indices in present and in the second wavefunctions" 1075 else if (imo<istart1.or.imo>iend1) then 1076 write(*,"(a,i6,a,i6)") "Error: Exceed valid range! The value should within",istart1," and",iend1 1077 else 1078 tmpval=0D0 1079 do jmo=istart2,iend2 1080 tmpval=tmpval+convmat(imo,jmo)**2*100D0 1081 write(*,"(i6,' Contribution:',f10.3,' % Coefficient:',f12.6)") jmo,convmat(imo,jmo)**2*100D0,convmat(imo,jmo) 1082 end do 1083 write(*,"(' Total:',f10.3,' %')") tmpval 1084 write(*,*) 1085 end if 1086end do 1087 1088if (iautointgrid==1) then 1089 radpot=radpotold 1090 sphpot=sphpotold 1091end if 1092end subroutine 1093 1094 1095 1096 1097 1098!!-------- Parse the output of (hyper)polarizability task of Gaussian to make it more understandable 1099subroutine parseGauPolar 1100use defvar 1101use util 1102implicit real*8 (a-h,o-z) 1103real*8 dipole(3),poltens(3,3),hypoltens(3,3,3) !hypoltens2(3,3,3,3) 1104real*8 eigvecmat(3,3),eigval(3),freqval(100000) 1105character c200tmp*200,sepchar,c210tmp*210 1106character*20 :: form,formau="(a,f16.6)",formother="(a,1PE16.6)" 1107integer :: irdfreq=0,ides=6,iunit=1 1108poltens=0D0 1109hypoltens=0D0 1110write(*,*) "Note: This function only works for ""polar"" tasks of Gaussian 09 with #P" 1111do while(.true.) 1112 write(*,"(a,/)") " Select the type of your Gaussian task by option 1~6. 2,4,6 only give polarizability (alpha), & 1113 the others also give hyperpolarizability (beta). -1 can be chosen only if CPHF=RdFreq or polar=DCSHG is used" 1114 if (iunit==1) write(*,*) "-3 Set the unit in the output, current: a.u." 1115 if (iunit==2) write(*,*) "-3 Set the unit in the output, current: SI" 1116 if (iunit==3) write(*,*) "-3 Set the unit in the output, current: esu" 1117 if (ides==6) write(*,*) "-2 Set the output destination, current: Output to screen" 1118 if (ides==11) write(*,*) "-2 Set the output destination, current: polar.txt in current folder" 1119 if (irdfreq==1) write(*,*) "-1 Toggle if load frequency-dependent result for option 1, current: Yes" 1120 if (irdfreq==0) write(*,*) "-1 Toggle if load frequency-dependent result for option 1, current: No" 1121 write(*,*) "0 Return" 1122 write(*,*) "1 ""Polar"" + analytic 3-order deriv. (HF/DFT/Semi-empirical)" 1123 write(*,*) "2 ""Polar"" + analytic 2-order deriv. (MP2...)" 1124 write(*,*) "3 ""Polar=Cubic"" + analytic 2-order deriv." 1125 write(*,*) "4 ""Polar"" + analytic 1-order deriv. (CISD,QCISD,CCSD,MP3,MP4(SDQ)...)" 1126 write(*,*) "5 ""Polar=DoubleNumer"" or ""Polar=EnOnly"" + analytic 1-order deriv." 1127 write(*,*) "6 ""Polar"" + energy only (CCSD(T),QCISD(T),MP4(SDTQ),MP5...)" 1128 read(*,*) isel 1129 if (isel==-1) then 1130 if (irdfreq==1) then 1131 irdfreq=0 1132 else 1133 irdfreq=1 1134 end if 1135 else if (isel==-2) then 1136 write(*,*) "6: Output to screen" 1137 write(*,*) "11: Output to polar.txt in current folder" 1138 read(*,*) ides 1139 else if (isel==-3) then 1140 write(*,*) "1: Atomic unit" 1141 write(*,*) "2: SI unit (C^2*m^2/J for alpha, C^3*m^3/J for beta)" 1142 write(*,*) "3: esu" 1143 read(*,*) iunit 1144 else if (isel==0) then 1145 return 1146 else 1147 exit 1148 end if 1149end do 1150if (irdfreq==1.and.isel>=2) then 1151 write(*,*) "ERROR: Frequency-dependent values are only available for HF/DFT/semi-empirical!" 1152 return 1153end if 1154 1155open(10,file=filename,status="old") 1156if (ides==11) open(ides,file="polar.txt",status="replace") 1157 1158if (iunit==1) then 1159 write(ides,*) "Note: All units shown below are in a.u." 1160 form=formau 1161else if (iunit==2) then 1162 write(ides,*) "Note: All units shown below are in SI unit (C^2*m^2/J for alpha, C^3*m^3/J for beta)" 1163 form=formother 1164else if (iunit==3) then 1165 write(ides,*) "Note: All units shown below are in esu unit" 1166 form=formother 1167end if 1168write(ides,*) 1169 1170 1171! Dipole moment part (miu) 1172call loclabel(10,"Dipole moment (field-independent basis, Debye)",ifound) 1173if (ifound==1) then 1174 read(10,*) 1175 read(10,*) c200tmp,xtmp,c200tmp,ytmp,c200tmp,ztmp 1176 dipole(1)=xtmp/au2debye !Convert to a.u. 1177 dipole(2)=ytmp/au2debye 1178 dipole(3)=ztmp/au2debye 1179 if (iunit==2) dipole=dipole*8.47835D-30 1180 if (iunit==3) dipole=dipole*2.54175D-18 1181 dipnorm=dsqrt(sum(dipole**2)) 1182 write(ides,*) "Dipole moment:" 1183 if (iunit==1) then 1184 write(ides,"(' X,Y,Z=',3f12.6,' Norm=',f12.6)") dipole(:),dipnorm 1185 else 1186 write(ides,"(' X,Y,Z=',3(1PE15.6),' Norm=',1PE15.6)") dipole(:),dipnorm 1187 end if 1188 write(ides,*) 1189end if 1190 1191!Selecting the result at which frequency will be loaded 1192if (irdfreq==1) then 1193 nfreqval=0 1194 rewind(10) 1195 do while(.true.) 1196 call loclabel(10,"-- Alpha(-w,w) frequency",ifound,0) !Not rewind 1197 if (ifound==1) then 1198 read(10,"(46x,f12.6)") tmpval 1199 if ( any(freqval(1:nfreqval)==tmpval) ) exit !The output in stage 2 and stage 3 are identical, so determine if has entered stage 3 1200 nfreqval=nfreqval+1 1201 freqval(nfreqval)=tmpval 1202 else 1203 exit 1204 end if 1205 end do 1206 write(*,*) "Load which one? Input the index" 1207 do i=1,nfreqval 1208 if (freqval(i)>0) write(*,"(i8,' w=',f12.6,' (',f12.2,'nm )')") i,freqval(i),1240.7011D0/(freqval(i)*au2eV) 1209 if (freqval(i)==0) write(*,"(i8,' w=',f12.6,' ( Static )')") i,freqval(i) 1210 end do 1211 read(*,*) ifreq 1212 if (freqval(ifreq)>0) write(*,"(' Note: All Alpha and Beta below correspond to w=',f12.6,' (',f12.2,'nm )',/)") freqval(ifreq),1240.7011D0/(freqval(ifreq)*au2eV) 1213 if (freqval(ifreq)==0) write(*,"(' Note: All Alpha and Beta below correspond to w=',f12.6,' ( Static )',/)") freqval(ifreq) 1214 rewind(10) !Move to the beginning 1215end if 1216 1217 1218!!!! Polarizability part (Alpha) 1219if (isel==1.or.isel==4.or.isel==5) then 1220 if (isel==1) then 1221 if (irdfreq==0) then 1222 call loclabel(10,"SCF Polarizability",ifound) 1223 else if (irdfreq==1) then 1224 do i=1,ifreq 1225 call loclabel(10,"SCF Polarizability",ifound,0) !Not rewind 1226 read(10,*) 1227 end do 1228 backspace(10) 1229 end if 1230 else if (isel==4.or.isel==5) then 1231 call loclabel(10,"Isotropic polarizability",ifound) 1232 end if 1233 call skiplines(10,2) 1234 read(10,*) rnouse,poltens(1,1) 1235 read(10,*) rnouse,poltens(2,1),poltens(2,2) 1236 read(10,*) rnouse,poltens(3,1),poltens(3,2),poltens(3,3) 1237else if (isel==2.or.isel==3) then 1238 call loclabel(10,"Exact polarizability") 1239 read(10,"(23x,6f8.3)") poltens(1,1),poltens(2,1),poltens(2,2),poltens(3,1),poltens(3,2),poltens(3,3) 1240else if (isel==6) then !Find result from archive part 1241 sepchar="\" 1242 if (isys==1) sepchar="|" 1243 do while(.true.) 1244 read(10,"(1x,a70)") c210tmp(1:70) !Combine two lines 1245 read(10,"(1x,a70)") c210tmp(71:140) 1246 read(10,"(1x,a70)") c210tmp(141:210) 1247 if (index(c210tmp(1:140),sepchar//"Polar=")/=0) exit 1248 backspace(10) 1249 backspace(10) 1250 end do 1251 do istart=1,204 1252 if (c210tmp(istart:istart+6)==sepchar//"Polar=") exit 1253 end do 1254 do iend=istart+1,210 1255 if (c210tmp(iend:iend)=="|".or.c200tmp(iend:iend)=="\") exit 1256 end do 1257 if (istart <= 204) then 1258 c210tmp(1:istart+6)=" " !Clean other information so that the data can be read in free format 1259 else 1260 c210tmp(1:210)=" " !Clean other information so that the data can be read in free format 1261 !0 is usually stderr 1262 write(0,*) "Wrong input! Please double check your Gaussian polar output!" 1263 endif 1264 c210tmp(iend:)=" " 1265 read(c210tmp,*) poltens(1,1),poltens(2,1),poltens(2,2),poltens(3,1),poltens(3,2),poltens(3,3) 1266end if 1267poltens(1,2)=poltens(2,1) 1268poltens(1,3)=poltens(3,1) 1269poltens(2,3)=poltens(3,2) 1270!Convert to other unit 1271if (iunit==2) poltens=poltens*1.6488D-41 1272if (iunit==3) poltens=poltens*1.4819D-25 1273if (irdfreq==0) write(ides,*) "Static polarizability:" 1274if (irdfreq==1) write(ides,*) "Frequency-dependent polarizability:" 1275write(ides,form) " XX=",poltens(1,1) 1276write(ides,form) " XY=",poltens(1,2) 1277write(ides,form) " YY=",poltens(2,2) 1278write(ides,form) " XZ=",poltens(1,3) 1279write(ides,form) " YZ=",poltens(2,3) 1280write(ides,form) " ZZ=",poltens(3,3) 1281alphaiso=(poltens(1,1)+poltens(2,2)+poltens(3,3))/3D0 1282write(ides,form) ' Isotropic average polarizability:',alphaiso 1283term1=(poltens(1,1)-poltens(2,2))**2 + (poltens(1,1)-poltens(3,3))**2 + (poltens(2,2)-poltens(3,3))**2 1284term2=6*(poltens(1,2)**2+poltens(1,3)**2+poltens(2,3)**2) 1285alphaani1=dsqrt((term1+term2)/2D0) 1286write(ides,form) ' Polarizability anisotropy (definition 1):',alphaani1 1287alphaani2=dsqrt(term1/2D0) 1288write(ides,form) ' Polarizability anisotropy (definition 2):',alphaani2 1289call diagmat(poltens,eigvecmat,eigval,300,1D-10) 1290call sortr8(eigval) 1291if (iunit==1) then 1292 write(ides,"(a,3f13.5)") ' Eigenvalues of polarizability tensor:',eigval 1293else 1294 write(ides,"(a,3(1PE13.5))") ' Eigenvalues of polarizability tensor:',eigval 1295end if 1296alphaani3=eigval(3)-(eigval(1)+eigval(2))/2D0 1297write(ides,form) ' Polarizability anisotropy (definition 3):',alphaani3 1298write(ides,*) 1299 1300 1301!!!! First hyperpolarizability part (Beta) 1302if (isel==1.or.isel==3.or.isel==5) then 1303 if (irdfreq==0) then 1304 if (isel==1) then 1305 call loclabel(10,"SCF Static Hyperpolarizability",ifound) 1306 else if (isel==3) then 1307 call loclabel(10,"Final packed hyperpolarizability",ifound) 1308 else if (isel==5) then 1309 call loclabel(10,"Static Hyperpolarizability",ifound) 1310 end if 1311 call skiplines(10,3) 1312 read(10,*) rnouse,hypoltens(1,1,1) !XXX 1313 call skiplines(10,2) 1314 read(10,*) rnouse,hypoltens(1,1,2) !XXY 1315 read(10,*) rnouse,hypoltens(1,2,2),hypoltens(2,2,2) !XYY,YYY 1316 call skiplines(10,2) 1317 read(10,*) rnouse,hypoltens(1,1,3) !XXZ 1318 read(10,*) rnouse,hypoltens(1,2,3),hypoltens(2,2,3) !XYZ,YYZ 1319 read(10,*) rnouse,hypoltens(1,3,3),hypoltens(2,3,3),hypoltens(3,3,3) !XZZ,YZZ,ZZZ 1320 !XYX=YXX =XXY 1321 hypoltens(1,2,1)=hypoltens(1,1,2) 1322 hypoltens(2,1,1)=hypoltens(1,1,2) 1323 !YXY=YYX =XYY 1324 hypoltens(2,1,2)=hypoltens(1,2,2) 1325 hypoltens(2,2,1)=hypoltens(1,2,2) 1326 !ZXX=XZX =XXZ 1327 hypoltens(3,1,1)=hypoltens(1,1,3) 1328 hypoltens(1,3,1)=hypoltens(1,1,3) 1329 !XZY=YXZ=ZYX=YXZ=YZX =XYZ 1330 hypoltens(1,3,2)=hypoltens(1,2,3) 1331 hypoltens(2,1,3)=hypoltens(1,2,3) 1332 hypoltens(3,2,1)=hypoltens(1,2,3) 1333 hypoltens(2,1,3)=hypoltens(1,2,3) 1334 hypoltens(2,3,1)=hypoltens(1,2,3) 1335 !ZYY=YZY =YYZ 1336 hypoltens(3,2,2)=hypoltens(2,2,3) 1337 hypoltens(2,3,2)=hypoltens(2,2,3) 1338 !ZXZ,ZZX =XZZ 1339 hypoltens(3,1,3)=hypoltens(1,3,3) 1340 hypoltens(3,3,1)=hypoltens(1,3,3) 1341 !ZYZ=ZZY =YZZ 1342 hypoltens(3,2,3)=hypoltens(2,3,3) 1343 hypoltens(3,3,2)=hypoltens(2,3,3) 1344 write(ides,"(a,/)") " Note: It is well known that the sign of hyperpolarizability of Gaussian 09 should be multiplied by -1, the outputs below have already been corrected" 1345 hypoltens=-hypoltens 1346 !Convert to other unit 1347 if (iunit==2) hypoltens=hypoltens*3.20636D-53 1348 if (iunit==3) hypoltens=hypoltens*8.63922D-33 1349 write(ides,*) "Static first hyperpolarizability:" 1350 write(ides,form) " XXX=",hypoltens(1,1,1) 1351 write(ides,form) " XXY=",hypoltens(1,1,2) 1352 write(ides,form) " XYY=",hypoltens(1,2,2) 1353 write(ides,form) " YYY=",hypoltens(2,2,2) 1354 write(ides,form) " XXZ=",hypoltens(1,1,3) 1355 write(ides,form) " XYZ=",hypoltens(1,2,3) 1356 write(ides,form) " YYZ=",hypoltens(2,2,3) 1357 write(ides,form) " XZZ=",hypoltens(1,3,3) 1358 write(ides,form) " YZZ=",hypoltens(2,3,3) 1359 write(ides,form) " ZZZ=",hypoltens(3,3,3) 1360 write(ides,*) 1361 betaX=hypoltens(1,1,1)+hypoltens(1,2,2)+hypoltens(1,3,3) 1362 betaY=hypoltens(2,2,2)+hypoltens(2,1,1)+hypoltens(2,3,3) 1363 betaZ=hypoltens(3,3,3)+hypoltens(3,2,2)+hypoltens(3,1,1) 1364 if (iunit==1) then 1365 write(ides,"(' Beta_X=',f15.5,' Beta_Y=',f15.5,' Beta_Z=',f15.5)") betaX,betaY,betaZ 1366 else 1367 write(ides,"(' Beta_X=',1PE15.5,' Beta_Y=',1PE15.5,' Beta_Z=',1PE15.5)") betaX,betaY,betaZ 1368 end if 1369 write(ides,form) " Magnitude of first hyperpolarizability:",dsqrt(betaX**2+betaY**2+betaZ**2) 1370 betaprj=(betaX*dipole(1)+betaY*dipole(2)+betaZ*dipole(3))/dipnorm 1371 write(ides,form) " Projection of beta on dipole moment:",betaprj 1372 write(ides,form) " Beta || :",betaprj/5D0*3D0 1373 write(ides,form) " Beta ||(z) :",betaZ/5D0*3D0 1374 write(ides,form) " Beta _|_(z) :",betaZ/5D0 1375 1376 !---For easier inspection in special use 1377! write(*,*) 1378! write(*,*) "=============" 1379! write(ides,form) ' Isotropic average polarizability:',alphaiso 1380! write(ides,form) ' Polarizability anisotropy (definition 2):',alphaani2 1381! write(ides,form) " Magnitude of first hyperpolarizability:",dsqrt(betaX**2+betaY**2+betaZ**2) 1382! write(ides,form) " Beta || :",betaprj/5D0*3D0 1383! read(*,*) 1384 1385 else if (irdfreq==1) then !Frequency-dependent hyperpolarizability, only available for HF/DFT/semi-empirical 1386 write(*,*) "Loading which type of hyperpolarizability?" 1387 write(*,*) "1: Beta(-w;w,0) 2: Beta(-2w;w,w)" 1388 write(*,*) "Note: Option 2 is meaningless if ""DCSHG"" keyword was not used" 1389 read(*,*) ibeta 1390 rewind(10) 1391 1392 if (ibeta==1) then !Beta(-w;w,0) case 1393 call loclabel(10,"-- Beta(-w,w,0) frequency",ifound) 1394 do i=1,ifreq 1395 call loclabel(10,"-- Beta(-w,w,0) frequency",ifound,0) 1396 read(10,*) 1397 end do 1398 read(10,*) 1399 read(10,*) c200tmp,hypoltens(1,1,1) 1400 read(10,*) c200tmp,hypoltens(2,1,1) 1401 read(10,*) c200tmp,hypoltens(2,2,1) 1402 read(10,*) c200tmp,hypoltens(3,1,1) 1403 read(10,*) c200tmp,hypoltens(3,2,1) 1404 read(10,*) c200tmp,hypoltens(3,3,1) 1405 hypoltens(1,2,1)=hypoltens(2,1,1) 1406 hypoltens(1,3,1)=hypoltens(3,1,1) 1407 hypoltens(2,3,1)=hypoltens(3,2,1) 1408 read(10,*) c200tmp,hypoltens(1,1,2) 1409 read(10,*) c200tmp,hypoltens(2,1,2) 1410 read(10,*) c200tmp,hypoltens(2,2,2) 1411 read(10,*) c200tmp,hypoltens(3,1,2) 1412 read(10,*) c200tmp,hypoltens(3,2,2) 1413 read(10,*) c200tmp,hypoltens(3,3,2) 1414 hypoltens(1,2,2)=hypoltens(2,1,2) 1415 hypoltens(1,3,2)=hypoltens(3,1,2) 1416 hypoltens(2,3,2)=hypoltens(3,2,2) 1417 read(10,*) c200tmp,hypoltens(1,1,3) 1418 read(10,*) c200tmp,hypoltens(2,1,3) 1419 read(10,*) c200tmp,hypoltens(2,2,3) 1420 read(10,*) c200tmp,hypoltens(3,1,3) 1421 read(10,*) c200tmp,hypoltens(3,2,3) 1422 read(10,*) c200tmp,hypoltens(3,3,3) 1423 hypoltens(1,2,3)=hypoltens(2,1,3) 1424 hypoltens(1,3,3)=hypoltens(3,1,3) 1425 hypoltens(2,3,3)=hypoltens(3,2,3) 1426 write(ides,"(a,/)") " Note: It is well known that the sign of hyperpolarizability of Gaussian 09 should be multiplied by -1, the outputs below have already been corrected" 1427 hypoltens=-hypoltens 1428 !Convert to other unit 1429 if (iunit==2) hypoltens=hypoltens*3.20636D-53 1430 if (iunit==3) hypoltens=hypoltens*8.63922D-33 1431 write(ides,*) "Frequency-dependent first hyperpolarizability Beta(-w;w,0)" 1432 write(ides,form) " XXX= ",hypoltens(1,1,1) 1433 write(ides,form) " XYX= YXX=",hypoltens(1,2,1) 1434 write(ides,form) " YYX= ",hypoltens(2,2,1) 1435 write(ides,form) " XZX= ZXX=",hypoltens(1,3,1) 1436 write(ides,form) " YZX= ZYX=",hypoltens(2,3,1) 1437 write(ides,form) " ZZX= ",hypoltens(3,3,1) 1438 write(ides,form) " XXY= ",hypoltens(1,1,2) 1439 write(ides,form) " XYY= YXY=",hypoltens(1,2,2) 1440 write(ides,form) " YYY= ",hypoltens(2,2,2) 1441 write(ides,form) " XZY= ZXY=",hypoltens(1,3,2) 1442 write(ides,form) " YZY= ZYY=",hypoltens(2,3,2) 1443 write(ides,form) " ZZY= ",hypoltens(3,3,2) 1444 write(ides,form) " XXZ= ",hypoltens(1,1,3) 1445 write(ides,form) " XYZ= YXZ=",hypoltens(1,2,3) 1446 write(ides,form) " YYZ= ",hypoltens(2,2,3) 1447 write(ides,form) " XZZ= ZXZ=",hypoltens(1,3,3) 1448 write(ides,form) " YZZ= ZYZ=",hypoltens(2,3,3) 1449 write(ides,form) " ZZZ= ",hypoltens(3,3,3) 1450 write(ides,*) 1451 1452 else if (ibeta==2) then !used DCSHG, parsing Beta(-2w;w,w) 1453 call loclabel(10,"-- Beta(w,w,-2w) frequency",ifound) 1454 do i=1,ifreq 1455 call loclabel(10,"-- Beta(w,w,-2w) frequency",ifound,0) 1456 read(10,*) 1457 end do 1458 read(10,*) 1459 read(10,*) c200tmp,hypoltens(1,1,1) 1460 read(10,*) c200tmp,hypoltens(1,1,2) 1461 hypoltens(1,2,1)=hypoltens(1,1,2) 1462 read(10,*) c200tmp,hypoltens(1,2,2) 1463 read(10,*) c200tmp,hypoltens(1,1,3) 1464 hypoltens(1,3,1)=hypoltens(1,1,3) 1465 read(10,*) c200tmp,hypoltens(1,2,3) 1466 hypoltens(1,3,2)=hypoltens(1,2,3) 1467 read(10,*) c200tmp,hypoltens(1,3,3) 1468 read(10,*) c200tmp,hypoltens(2,1,1) 1469 read(10,*) c200tmp,hypoltens(2,2,1) 1470 hypoltens(2,1,2)=hypoltens(2,2,1) 1471 read(10,*) c200tmp,hypoltens(2,2,2) 1472 read(10,*) c200tmp,hypoltens(2,1,3) 1473 hypoltens(2,3,1)=hypoltens(2,1,3) 1474 read(10,*) c200tmp,hypoltens(2,2,3) 1475 hypoltens(2,3,2)=hypoltens(2,2,3) 1476 read(10,*) c200tmp,hypoltens(2,3,3) 1477 read(10,*) c200tmp,hypoltens(3,1,1) 1478 read(10,*) c200tmp,hypoltens(3,1,2) 1479 hypoltens(3,2,1)=hypoltens(3,1,2) 1480 read(10,*) c200tmp,hypoltens(3,2,2) 1481 read(10,*) c200tmp,hypoltens(3,1,3) 1482 hypoltens(3,3,1)=hypoltens(3,1,3) 1483 read(10,*) c200tmp,hypoltens(3,2,3) 1484 hypoltens(3,3,2)=hypoltens(3,2,3) 1485 read(10,*) c200tmp,hypoltens(3,3,3) 1486 write(ides,"(a,/)") " Note: It is well known that the sign of hyperpolarizability of Gaussian 09 should be multiplied by -1, the outputs below have already been corrected" 1487 hypoltens=-hypoltens 1488 !Convert to other unit 1489 if (iunit==2) hypoltens=hypoltens*3.20636D-53 1490 if (iunit==3) hypoltens=hypoltens*8.63922D-33 1491 if (ibeta==2) write(ides,*) "Frequency-dependent first hyperpolarizability Beta(-2w;w,w)" 1492 write(ides,form) " XXX= ",hypoltens(1,1,1) 1493 write(ides,form) " YXX= ",hypoltens(2,1,1) 1494 write(ides,form) " ZXX= ",hypoltens(3,1,1) 1495 write(ides,form) " XYX= XXY=",hypoltens(1,2,1) 1496 write(ides,form) " YYX= YXY=",hypoltens(2,2,1) 1497 write(ides,form) " ZYX= ZXY=",hypoltens(3,2,1) 1498 write(ides,form) " XYY= ",hypoltens(1,2,2) 1499 write(ides,form) " YYY= ",hypoltens(2,2,2) 1500 write(ides,form) " ZYY= ",hypoltens(3,2,2) 1501 write(ides,form) " XZX= XXZ=",hypoltens(1,3,1) 1502 write(ides,form) " YZX= YXZ=",hypoltens(2,3,1) 1503 write(ides,form) " ZZX= ZXZ=",hypoltens(3,3,1) 1504 write(ides,form) " XZY= XYZ=",hypoltens(1,3,2) 1505 write(ides,form) " YZY= YYZ=",hypoltens(2,3,2) 1506 write(ides,form) " ZZY= ZYZ=",hypoltens(3,3,2) 1507 write(ides,form) " XZZ= ",hypoltens(1,3,3) 1508 write(ides,form) " YZZ= ",hypoltens(2,3,3) 1509 write(ides,form) " ZZZ= ",hypoltens(3,3,3) 1510 write(ides,*) 1511 end if 1512 1513 betaX=0 1514 betaY=0 1515 betaZ=0 1516 do j=1,3 1517 betaX=betaX+(hypoltens(1,j,j)+hypoltens(j,j,1)+hypoltens(j,1,j))/3 1518 betaY=betaY+(hypoltens(2,j,j)+hypoltens(j,j,2)+hypoltens(j,2,j))/3 1519 betaZ=betaZ+(hypoltens(3,j,j)+hypoltens(j,j,3)+hypoltens(j,3,j))/3 1520 end do 1521 if (iunit==1) then 1522 write(ides,"(' Beta_X=',f15.5,' Beta_Y=',f15.5,' Beta_Z=',f15.5)") betaX,betaY,betaZ 1523 else 1524 write(ides,"(' Beta_X=',1PE15.5,' Beta_Y=',1PE15.5,' Beta_Z=',1PE15.5)") betaX,betaY,betaZ 1525 end if 1526 write(ides,form) " Magnitude of first hyperpolarizability:",dsqrt(betaX**2+betaY**2+betaZ**2) 1527 betaprj=(betaX*dipole(1)+betaY*dipole(2)+betaZ*dipole(3))/dipnorm 1528 write(ides,form) " Projection of beta on dipole moment:",betaprj 1529 write(ides,form) " Beta || :",betaprj*3D0/5D0 1530 write(ides,form) " Beta ||(z) :",betaZ*3D0/5D0 1531 beta_per=0 1532 do j=1,3 1533 beta_per=beta_per+(2*hypoltens(3,j,j)+2*hypoltens(j,j,3)-3*hypoltens(j,3,j))/5 1534 end do 1535 write(ides,form) " Beta _|_(z) :",beta_per 1536 write(*,*) 1537 end if 1538end if 1539 1540close(10) 1541if (ides==11) close(ides) 1542end subroutine 1543 1544 1545 1546!!--------- Sum-over-states (SOS) calculation for (hyper)polarizability 1547!Programmed based on the formulae in Sasagane et al. J. Chem. Phys., 99, 3738 (1993) 1548subroutine SOS 1549use defvar 1550use util 1551implicit real*8 (a-h,o-z) 1552character transmodestr*80,c80tmp*80,c200tmp*80 1553character :: dirlab(3)=(/ "X","Y","Z" /) 1554real*8,allocatable :: trandip(:,:,:) !Transition dipole moment between i and j in X,Y,Z. 0 corresponds to ground state 1555real*8,allocatable :: excene(:) !Excitation energy 1556real*8 :: alpha(3,3),beta(3,3,3),gamma(3,3,3,3),delta(3,3,3,3,3) 1557real*8 eigval(3),eigvecmat(3,3),tmpw(5) 1558real*8,allocatable :: freqlist(:,:) !Store the frequency to be calculated for beta and gamma 1559integer tmpdir(5),arrb(6,3),arrg(24,4),arrd(120,5),dir1,dir2,dir3,dir4,dir5 1560 1561write(*,*) "Loading data..." 1562open(10,file=filename,status="old") 1563call loclabel(10,"Gaussian",igauout) 1564rewind(10) 1565if (igauout==1) then !Load excitation energies and <0|r|n> 1566 write(*,*) "This is a Gaussian output file" 1567 call loclabel(10,"Excitation energies and oscillator strengths:") 1568 read(10,*) 1569 nstates=0 !The number of transition modes 1570 do while(.true.) 1571 call loclabel(10,"Excited State",ifound,0) 1572 if (ifound==1) then 1573 nstates=nstates+1 1574 read(10,*) 1575 else 1576 exit 1577 end if 1578 end do 1579 allocate(trandip(0:nstates,0:nstates,3),excene(0:nstates)) 1580 trandip=0 1581 excene=0 1582 call loclabel(10,"Ground to excited state transition electric dipole moments") 1583 read(10,*) 1584 read(10,*) 1585 do istat=1,nstates 1586 read(10,*) inouse,trandip(0,istat,:) !Transition dipole moment from ground state (0) to excited "istate" 1587 end do 1588 call loclabel(10,"Excitation Energies [eV] at current iteration:",ifound) 1589 if (ifound==1) then !Read the excitation energies shown in the last step of iteration process will be more accurate, but #P must be used 1590 ncyc=1 1591 do while(.true.) !Find the position of the last iteration 1592 call loclabel(10,"Excitation Energies [eV] at current iteration:",ifound,0) 1593 if (ifound==0) exit 1594 read(10,*) 1595 ncyc=ncyc+1 1596 end do 1597 rewind(10) 1598 do icyc=1,ncyc 1599 call loclabel(10,"Excitation Energies [eV] at current iteration:",ifound,0) 1600 end do 1601 read(10,*) 1602 do istat=1,nstates 1603 read(10,"(14x)",advance="no") 1604 read(10,*) excene(istat) 1605 end do 1606 else !Read excitation energies from normal output 1607 call loclabel(10,"Excitation energies and oscillator strengths:") 1608 do istat=1,nstates 1609 call loclabel(10,"Excited State",ifound,0) 1610 read(10,"(a)") transmodestr 1611 do i=10,70 1612 if (transmodestr(i:i+1)=="eV") exit 1613 end do 1614 read(transmodestr(i-10:i-1),*) excene(istat) !Read as eV 1615 end do 1616 end if 1617 call loclabel(10,"Dipole moment (field-independent basis, Debye)") 1618 read(10,*) 1619 read(10,*) c80tmp,xtmp,c80tmp,ytmp,c80tmp,ztmp 1620 trandip(0,0,1)=xtmp/au2debye 1621 trandip(0,0,2)=ytmp/au2debye 1622 trandip(0,0,3)=ztmp/au2debye 1623 ionlyalpha=1 1624else !Load excitation energies and all <m|r|n> from plain text file 1625 ionlyalpha=0 1626 read(10,*) nstates 1627 if (nstates<0) ionlyalpha=1 !Only calculate polarizability, not hyperpolarizability, so will not read transition dipole moments among excited states 1628 nstates=abs(nstates) 1629 allocate(trandip(0:nstates,0:nstates,3),excene(0:nstates)) 1630 do i=1,nstates !Read as eV 1631 read(10,*) inouse,excene(i) 1632 end do 1633 do i=0,nstates !i-i corresponds to dipole moment of corresponding state 1634 do j=i,nstates 1635 read(10,*) inouse,inouse,trandip(i,j,:) 1636 end do 1637 if (ionlyalpha==1) exit 1638 end do 1639end if 1640close(10) 1641excene(0)=0 1642excene=excene/au2eV 1643!Transition dipole moments are loaded as upper trigonal matrix, now we convert it as symmetry matrix 1644do i=0,nstates 1645 do j=i+1,nstates 1646 trandip(j,i,:)=trandip(i,j,:) 1647 end do 1648end do 1649 1650!! Output some information 1651! write(*,*) " State# Exc.ene.(a.u.) Transition dipole moment in X,Y,Z (a.u.)" 1652! do istat=1,nstates 1653! write(*,"(i8,f20.6,3f15.6)") istat,excene(istat),trandip(0,istat,:) 1654! end do 1655do istat=1,nstates 1656 write(*,"(' State',i6,' Excitation energy:',f12.5,' a.u.',f14.6,' eV')") istat,excene(istat),excene(istat)*au2eV 1657end do 1658write(*,"(' There are',i6,' excited states')") nstates 1659write(*,"(' Dipole moment of ground state:'3f12.5,' a.u.')") trandip(0,0,:) 1660write(*,*) "NOTE: All units used in this function is a.u." 1661!Gaussian output file is impossible to provide <m|r|n>, even if alltransitiondensities is used for CIS, it doesn't output <m|r|m> 1662do while(.true.) 1663write(*,*) 1664write(*,*) "0 Return" 1665write(*,*) "1 Calculate polarizability (alpha)" 1666if (ionlyalpha==0) write(*,*) "2 Calculate first hyperpolarizability (beta)" 1667if (ionlyalpha==0) write(*,*) "3 Calculate second hyperpolarizability (gamma)" 1668if (ionlyalpha==0) write(*,*) "4 Calculate third hyperpolarizability (delta)" 1669write(*,*) "5 Show the variation of alpha w.r.t. the number of states in consideration" 1670if (ionlyalpha==0) write(*,*) "6 Show the variation of beta w.r.t. the number of states in consideration" 1671if (ionlyalpha==0) write(*,*) "7 Show the variation of gamma w.r.t. the number of states in consideration" 1672write(*,*) "15 Calculate alpha in a range of frequencies" 1673if (ionlyalpha==0) write(*,*) "16 Calculate beta in a range of frequencies" 1674if (ionlyalpha==0) write(*,*) "17 Calculate gamma in a range of frequencies" 1675read(*,*) isel 1676 1677if (isel==0) then 1678 return 1679else if (isel==1.or.isel==5.or.isel==15) then ! Calculate polarizability 1680 if (isel==1.or.isel==5) then 1681 write(*,*) "Input frequency of external field w for alpha(-w;w)" 1682 write(*,*) "e.g. 0.25" 1683 write(*,*) "Note: Negative value means using nm as unit, e.g. -693.5" 1684 read(*,*) freqbeg 1685 if (freqbeg<0) freqbeg=1240.7011D0/au2eV/abs(freqbeg) 1686 freqend=freqbeg 1687 freqstep=1 1688 if (freqbeg/=0) then 1689 wavlen=1240.7011D0/(freqbeg*au2eV) 1690 write(*,"(' Wavelength of w:',f12.3,' nm',/)") wavlen 1691 end if 1692 else if (isel==15) then 1693 write(*,*) "Input lower and upper limits as well as stepsize of w for alpha(-w;w)" 1694 write(*,*) "e.g. 0.2,0.5,0.01" 1695 read(*,*) freqbeg,freqend,freqstep 1696 end if 1697 1698 if (isel==1) then 1699 istart=nstates 1700 iend=nstates 1701 else if (isel==5) then 1702 istart=1 1703 iend=nstates 1704 open(10,file="alpha_n.txt",status="replace") 1705 else if (isel==15) then 1706 istart=nstates 1707 iend=nstates 1708 open(10,file="alpha_w.txt",status="replace") 1709 end if 1710 1711 write(*,*) "Please wait..." 1712 do numstat=istart,iend 1713 do freq=freqbeg,freqend,freqstep 1714 do idir=1,3 1715 do jdir=1,3 1716 tmpval=0 1717 do istat=1,numstat 1718 tmpval1=trandip(0,istat,idir)*trandip(istat,0,jdir)/(excene(istat)-freq) 1719 tmpval2=trandip(0,istat,jdir)*trandip(istat,0,idir)/(excene(istat)+freq) 1720 tmpval=tmpval+tmpval1+tmpval2 1721 end do 1722 alpha(idir,jdir)=tmpval 1723 end do 1724 end do 1725 alphaiso=(alpha(1,1)+alpha(2,2)+alpha(3,3))/3D0 1726 term1=(alpha(1,1)-alpha(2,2))**2 + (alpha(1,1)-alpha(3,3))**2 + (alpha(2,2)-alpha(3,3))**2 1727 term2=6*(alpha(1,2)**2+alpha(1,3)**2+alpha(2,3)**2) 1728 alphaani1=dsqrt((term1+term2)/2D0) 1729 call diagmat(alpha,eigvecmat,eigval,300,1D-10) 1730 call sortr8(eigval) 1731 alphaani2=eigval(3)-(eigval(1)+eigval(2))/2D0 1732 1733 if (isel==1) then 1734 write(*,*) "Polarizability tensor:" 1735 write(*,*) " 1 2 3" 1736 do idir=1,3 1737 write(*,"(i3,3f15.6)") idir,alpha(idir,:) 1738 end do 1739 write(*,"(' Isotropic average polarizability:',f15.6)") alphaiso 1740 write(*,"(' Polarizability anisotropy (definition 1):',f15.6)") alphaani1 1741 write(*,"(' Eigenvalues:',3f15.6)") eigval(:) 1742 write(*,"(' Polarizability anisotropy (definition 2):',f15.6)") alphaani2 1743 else if (isel==5) then 1744 write(10,"(i6,9f15.6)") numstat,alphaiso,alphaani1,alphaani2,alpha(1,1),alpha(2,2),alpha(3,3),alpha(1,2),alpha(1,3),alpha(2,3) 1745 else if (isel==15) then 1746 write(10,"(f12.6,9(1PE14.5))") freq,alphaiso,alphaani1,alphaani2,alpha(1,1),alpha(2,2),alpha(3,3),alpha(1,2),alpha(1,3),alpha(2,3) 1747 end if 1748 end do 1749 end do 1750 if (isel==5.or.isel==15) then 1751 close(10) 1752 if (isel==5) write(*,*) "Done! The result has been outputted to alpha_n.txt in current folder" 1753 if (isel==15) write(*,*) "Done! The result has been outputted to alpha_w.txt in current folder" 1754 write(*,*) "The correspondence between columns and information in this file is as follows" 1755 if (isel==5) write(*,*) "Column 1: The number of states in consideration" 1756 if (isel==15) write(*,*) "Column 1: Frequency of external field" 1757 write(*,*) "Column 2: Isotropic average polarizability" 1758 write(*,*) "Column 3: Polarizability anisotropy (definition 1)" 1759 write(*,*) "Column 4: Polarizability anisotropy (definition 2)" 1760 write(*,*) "Column 5: XX" 1761 write(*,*) "Column 6: YY" 1762 write(*,*) "Column 7: ZZ" 1763 write(*,*) "Column 8: XY" 1764 write(*,*) "Column 9: XZ" 1765 write(*,*) "Column 10: YZ" 1766 end if 1767 1768! Calculate first hyperpolarizability 1769else if (isel==2.or.isel==6.or.isel==16) then 1770 if (allocated(freqlist)) deallocate(freqlist) 1771 if (isel==2.or.isel==6) then 1772 nfreq=1 1773 allocate(freqlist(1,2)) 1774 write(*,*) "Input frequency of external field w1 and w2 for beta(-(w1+w2);w1,w2)" 1775 write(*,*) "e.g. 0.25,0.13" 1776 write(*,*) "Note: Negative values mean using nm as unit, e.g. -693,0" 1777 read(*,*) freqlist(1,:) 1778 where(freqlist<0) freqlist=1240.7011D0/au2eV/abs(freqlist) 1779 if (freqlist(1,1)/=0) then 1780 wavlen1=1240.7011D0/(freqlist(1,1)*au2eV) 1781 write(*,"(' Wavelength of w1:',f12.3,' nm')") wavlen1 1782 end if 1783 if (freqlist(1,2)/=0) then 1784 wavlen2=1240.7011D0/(freqlist(1,2)*au2eV) 1785 write(*,"(' Wavelength of w2:',f12.3,' nm')") wavlen2 1786 end if 1787 write(*,*) 1788 else if (isel==16) then 1789 write(*,*) "Input the file recording frequency list, e.g. C:\freqlist.txt" 1790 write(*,"(a)") " The file should contain two columns, corresponding to frequency of w1 and w2 in a.u., respectively" 1791 do while(.true.) 1792 read(*,"(a)") c200tmp 1793 inquire(file=c200tmp,exist=alive) 1794 if (alive) exit 1795 write(*,*) "Cannot find the file, input again" 1796 end do 1797 open(10,file=c200tmp,status="old") 1798 nfreq=totlinenum(10,1) 1799 allocate(freqlist(nfreq,2)) 1800 do ifreq=1,nfreq 1801 read(10,*) freqlist(ifreq,:) 1802 end do 1803 close(10) 1804 write(*,*) "The frequencies loaded:" 1805 do ifreq=1,nfreq 1806 write(*,"(' #',i5,' w1=',f10.5,' a.u.',f10.3,' nm w2=',f10.5' a.u.',f10.3,' nm')") & 1807 ifreq,freqlist(ifreq,1),1240.7011D0/(freqlist(ifreq,1)*au2eV),freqlist(ifreq,2),1240.7011D0/(freqlist(ifreq,2)*au2eV) 1808 end do 1809 write(*,*) 1810 end if 1811 1812 if (isel==2) then 1813 istart=nstates 1814 iend=nstates 1815 else if (isel==6) then 1816 istart=1 1817 iend=nstates 1818 open(10,file="beta_n.txt",status="replace") 1819 else if (isel==16) then 1820 istart=nstates 1821 iend=nstates 1822 open(10,file="beta_w.txt",status="replace") 1823 end if 1824 1825 write(*,*) "Please wait..." 1826 call fullarrange(arrb,6,3) !Generate full arrangement matrix (3!=6 :3) for beta, each row corresponds to one permutation, e.g. 231 1827 do numstat=istart,iend 1828 do ifreq=1,nfreq 1829 freq1=freqlist(ifreq,1) 1830 freq2=freqlist(ifreq,2) 1831 freqtot=freq1+freq2 1832 tmpw(1:3)=(/ -freqtot,freq1,freq2 /) 1833 do idir=1,3 1834 do jdir=1,3 1835 do kdir=1,3 1836 1837 tmpval=0 1838 tmpdir(1:3)=(/ idir,jdir,kdir /) 1839 do iper=1,6 !Do permutation, arrb(1,:)=1,2,3 1840 dir1=tmpdir(arrb(iper,1)) 1841 dir2=tmpdir(arrb(iper,2)) 1842 dir3=tmpdir(arrb(iper,3)) 1843 w0=tmpw(arrb(iper,1)) 1844 w2=tmpw(arrb(iper,3)) 1845 do istat=1,numstat 1846 do jstat=1,numstat 1847 cen=trandip(istat,jstat,dir2) 1848 if (istat==jstat) cen=cen-trandip(0,0,dir2) 1849 tmpval=tmpval+trandip(0,istat,dir1)*cen*trandip(jstat,0,dir3)/(excene(istat)+w0)/(excene(jstat)-w2) 1850 end do 1851 end do 1852 end do 1853 beta(idir,jdir,kdir)=tmpval 1854 1855 !Below is the code manually considering each permutation, for teaching purpose, but foolish 1856! tmpval=0 1857! do istat=1,numstat 1858! do jstat=1,numstat 1859! ! Consider six permutations, i=-freqtot, j=freq1, k=freq2, the denominator is (+1th),(-3th) 1860! ! 1_3 1861! ! ijk; -freqtot,-freq2 1862! ! ikj; -freqtot,-freq1 1863! ! jik; +freq1 ,-freq2 1864! ! jki; +freq1 ,+freqtot 1865! ! kij; +freq2 ,-freq1 1866! ! kji; +freq2 ,+freqtot 1867! t1c=trandip(istat,jstat,jdir) 1868! t2c=trandip(istat,jstat,kdir) 1869! t3c=trandip(istat,jstat,idir) 1870! t4c=trandip(istat,jstat,kdir) 1871! t5c=trandip(istat,jstat,idir) 1872! t6c=trandip(istat,jstat,jdir) 1873! if (istat==jstat) then 1874! t1c=t1c-trandip(0,0,jdir) 1875! t2c=t2c-trandip(0,0,kdir) 1876! t3c=t3c-trandip(0,0,idir) 1877! t4c=t4c-trandip(0,0,kdir) 1878! t5c=t5c-trandip(0,0,idir) 1879! t6c=t6c-trandip(0,0,jdir) 1880! end if 1881! t1=trandip(0,istat,idir)*t1c*trandip(jstat,0,kdir) /(excene(istat)-freqtot)/(excene(jstat)-freq2) 1882! t2=trandip(0,istat,idir)*t2c*trandip(jstat,0,jdir) /(excene(istat)-freqtot)/(excene(jstat)-freq1) 1883! t3=trandip(0,istat,jdir)*t3c*trandip(jstat,0,kdir) /(excene(istat)+freq1) /(excene(jstat)-freq2) 1884! t4=trandip(0,istat,jdir)*t4c*trandip(jstat,0,idir) /(excene(istat)+freq1) /(excene(jstat)+freqtot) 1885! t5=trandip(0,istat,kdir)*t5c*trandip(jstat,0,jdir) /(excene(istat)+freq2) /(excene(jstat)-freq1) 1886! t6=trandip(0,istat,kdir)*t6c*trandip(jstat,0,idir) /(excene(istat)+freq2) /(excene(jstat)+freqtot) 1887! tmpval=tmpval+t1+t2+t3+t4+t5+t6 1888! end do 1889! end do 1890! beta(idir,jdir,kdir)=tmpval 1891 end do 1892 end do 1893 end do 1894 1895 betaX=0 1896 betaY=0 1897 betaZ=0 1898 do j=1,3 1899 betaX=betaX+(beta(1,j,j)+beta(j,j,1)+beta(j,1,j))/3 1900 betaY=betaY+(beta(2,j,j)+beta(j,j,2)+beta(j,2,j))/3 1901 betaZ=betaZ+(beta(3,j,j)+beta(j,j,3)+beta(j,3,j))/3 1902 end do 1903 betatot=dsqrt(betaX**2+betaY**2+betaZ**2) 1904 dipx=trandip(0,0,1) 1905 dipy=trandip(0,0,2) 1906 dipz=trandip(0,0,3) 1907 dipnorm=dsqrt(dipx**2+dipy**2+dipz**2) 1908 betaprj=(betaX*dipx+betaY*dipy+betaZ*dipz)/dipnorm 1909 beta_per=0 1910 do j=1,3 1911 beta_per=beta_per+(2*beta(3,j,j)+2*beta(j,j,3)-3*beta(j,3,j))/5 1912 end do 1913 1914 if (isel==2) then 1915 write(*,*) "First hyperpolarizability tensor:" 1916 do jdir=1,3 1917 do kdir=1,3 1918 do idir=1,3 1919 write(*,"(2x,3a,'=',f17.5,2x)",advance="no") dirlab(idir),dirlab(jdir),dirlab(kdir),beta(idir,jdir,kdir) 1920 if (idir==3) write(*,*) 1921 end do 1922 end do 1923 end do 1924 write(*,"(/,' Beta_X:',f17.5,' Beta_Y:',f17.5,' Beta_Z:',f17.5)") betaX,betaY,betaZ 1925 write(*,"(a,f17.5)") " Magnitude of beta:",betatot 1926 write(*,"(a,f17.5)") " Projection of beta on dipole moment:",betaprj 1927 write(*,"(a,f17.5)") " Beta || :",betaprj*3D0/5D0 1928 write(*,"(a,f17.5)") " Beta ||(z) :",betaZ*3D0/5D0 1929 write(*,"(a,f17.5)") " Beta _|_(z) :",beta_per 1930 else if (isel==6) then 1931 write(10,"(i6,8(1PE14.5))") numstat,betaX,betaY,betaZ,betatot,betaprj,betaprj*3D0/5D0,betaZ*3D0/5D0,beta_per 1932 else if (isel==16) then 1933 write(10,"(2f12.6,8(1PE14.5))") freq1,freq2,betaX,betaY,betaZ,betatot,betaprj,betaprj*3D0/5D0,betaZ*3D0/5D0,beta_per 1934 end if 1935 end do 1936 end do 1937 1938 if (isel==6.or.isel==16) then 1939 close(10) 1940 if (isel==6) then 1941 write(*,*) "Done! The result has been outputted to beta_n.txt in current folder" 1942 write(*,*) "The correspondence between columns and information in this file is as follows" 1943 write(*,*) "Column 1: The number of states in consideration" 1944 write(*,*) "Column 2: Beta_X" 1945 write(*,*) "Column 3: Beta_Y" 1946 write(*,*) "Column 4: Beta_Z" 1947 write(*,*) "Column 5: Magnitude of hyperpolarizability" 1948 write(*,*) "Column 6: Hyperpolarizability component along dipole moment direction" 1949 write(*,*) "Column 7: Beta ||" 1950 write(*,*) "Column 8: Beta ||(z)" 1951 write(*,*) "Column 9: Beta _|_(z)" 1952 else if (isel==16) then 1953 write(*,*) "Done! The result has been outputted to beta_w.txt in current folder" 1954 write(*,*) "The correspondence between columns and information in this file is as follows" 1955 write(*,*) "Column 1: Frequency of the first external field (w1)" 1956 write(*,*) "Column 2: Frequency of the second external field (w2)" 1957 write(*,*) "Column 3: Beta_X" 1958 write(*,*) "Column 4: Beta_Y" 1959 write(*,*) "Column 5: Beta_Z" 1960 write(*,*) "Column 6: Magnitude of hyperpolarizability" 1961 write(*,*) "Column 7: Hyperpolarizability component along dipole moment direction" 1962 write(*,*) "Column 8: Beta ||" 1963 write(*,*) "Column 9: Beta ||(z)" 1964 write(*,*) "Column 10: Beta _|_(z)" 1965 end if 1966 end if 1967 1968! Calculate second hyperpolarizability (gamma) 1969else if (isel==3.or.isel==7.or.isel==17) then 1970 if (allocated(freqlist)) deallocate(freqlist) 1971 if (isel==3.or.isel==7) then 1972 nfreq=1 1973 allocate(freqlist(1,3)) 1974 write(*,*) "Input frequency of external fields w1, w2, w3 for gamma(-(w1+w2+w3);w1,w2,w3)" 1975 write(*,*) "e.g. 0.13,0.13,0" 1976 write(*,*) "Note: Negative values mean using nm as unit, e.g. -693,0,0" 1977 read(*,*) freqlist(1,:) 1978 where(freqlist<0) freqlist=1240.7011D0/au2eV/abs(freqlist) 1979 do i=1,3 1980 if (freqlist(1,i)/=0) then 1981 wavlen=1240.7011D0/(freqlist(1,i)*au2eV) 1982 write(*,"(' Wavelength of w',i1,':',f12.3,' nm')") i,wavlen 1983 end if 1984 end do 1985 write(*,*) 1986 else if (isel==17) then 1987 write(*,*) "Input the file recording frequency list, e.g. C:\freqlist.txt" 1988 write(*,"(a)") " The file should contain three columns, corresponding to frequency of w1, w2 and w3 in a.u., respectively" 1989 do while(.true.) 1990 read(*,"(a)") c200tmp 1991 inquire(file=c200tmp,exist=alive) 1992 if (alive) exit 1993 write(*,*) "Cannot find the file, input again" 1994 end do 1995 open(10,file=c200tmp,status="old") 1996 nfreq=totlinenum(10,1) 1997 allocate(freqlist(nfreq,3)) 1998 do ifreq=1,nfreq 1999 read(10,*) freqlist(ifreq,:) 2000 end do 2001 close(10) 2002 write(*,*) "The frequencies loaded:" 2003 do ifreq=1,nfreq 2004 write(*,"(' #',i5,' w1=',f10.5,' a.u. w2=',f10.5,' a.u. w3=',f10.5,' a.u.')") ifreq,freqlist(ifreq,:) 2005 end do 2006 write(*,*) 2007 end if 2008 2009 if (isel==3.or.isel==17) then 2010 write(*,"(' Consider how many states? Should be <=',i6)") nstates 2011 read(*,*) istart 2012 iend=istart 2013 nstatstep=1 2014 else if (isel==7) then 2015 write(*,*) "Input upper limit and stepsize of the number of states" 2016 write(*,*) "e.g. 150,2" 2017 read(*,*) iend,nstatstep 2018 istart=1 2019 if (iend>nstates) iend=nstates 2020 end if 2021 if (isel==7) open(10,file="gamma_n.txt",status="replace") 2022 if (isel==17) open(10,file="gamma_w.txt",status="replace") 2023 2024 write(*,*) "Please wait..." 2025 call walltime(iwalltime1) 2026 call fullarrange(arrg,24,4) !Generate full arrangement matrix (4!=24 :4) for gamma, each row corresponds to one permutation, e.g. 2341 2027! do i=1,24 2028! write(*,"(i4,4i1)") arrg(i,:) 2029! end do 2030 iprog=0 2031 do numstat=istart,iend,nstatstep 2032 iprog=iprog+nstatstep 2033 if (isel==7) write(*,"(' Finished:',i5,' /',i5)") iprog,iend-istart+1 2034 do ifreq=1,nfreq 2035 freq1=freqlist(ifreq,1) 2036 freq2=freqlist(ifreq,2) 2037 freq3=freqlist(ifreq,3) 2038 freqtot=freq1+freq2+freq3 2039 tmpw(1:4)=(/ -freqtot,freq1,freq2,freq3 /) 2040 2041 do idir=1,3 !Cycle direction component 2042 do jdir=1,3 2043 do kdir=1,3 2044 do ldir=1,3 2045 2046 gamma1=0 2047 gamma2=0 2048 tmpdir(1:4)=(/ idir,jdir,kdir,ldir /) 2049 do iper=1,24 !Do permutation, arrg(1,:)=1,2,3,4 2050 dir1=tmpdir(arrg(iper,1)) 2051 dir2=tmpdir(arrg(iper,2)) 2052 dir3=tmpdir(arrg(iper,3)) 2053 dir4=tmpdir(arrg(iper,4)) 2054 w0=tmpw(arrg(iper,1)) 2055 w1=tmpw(arrg(iper,2)) 2056 w2=tmpw(arrg(iper,3)) 2057 w3=tmpw(arrg(iper,4)) 2058nthreads=getNThreads() 2059!$OMP PARALLEL SHARED(gamma1,gamma2) PRIVATE(istat,jstat,kstat,t1c,t2c,p1,p2,g1t,g2t) NUM_THREADS(nthreads) 2060 g1t=0 2061 g2t=0 2062!$OMP DO schedule(dynamic) 2063 do istat=1,numstat 2064 do jstat=1,numstat 2065 !Gamma 1 2066 do kstat=1,numstat 2067 t1c=trandip(istat,jstat,dir2) 2068 if (istat==jstat) t1c=t1c-trandip(0,0,dir2) 2069 t2c=trandip(jstat,kstat,dir3) 2070 if (jstat==kstat) t2c=t2c-trandip(0,0,dir3) 2071 p1=trandip(0,istat,dir1)*t1c*t2c*trandip(kstat,0,dir4) 2072 p2=(excene(istat)+w0)*(excene(jstat)-w2-w3)*(excene(kstat)-w3) 2073 g1t=g1t+p1/p2 2074 end do 2075 !Gamma 2 2076 p1=trandip(0,istat,dir1)*trandip(istat,0,dir2)*trandip(0,jstat,dir3)*trandip(jstat,0,dir4) 2077 p2=(excene(istat)+w0)*(excene(istat)-w1)*(excene(jstat)-w3) !(excene(jstat)-w3) can also be (excene(jstat)+w2), they are equivalent 2078 g2t=g2t+p1/p2 2079 end do 2080 end do 2081!$OMP END DO 2082!$OMP CRITICAL 2083 gamma1=gamma1+g1t 2084 gamma2=gamma2+g2t 2085!$OMP END CRITICAL 2086!$OMP END PARALLEL 2087 end do 2088 gamma(idir,jdir,kdir,ldir)=gamma1-gamma2 2089 2090 end do 2091 end do 2092 end do 2093 end do 2094 2095 gammaX=0 2096 gammaY=0 2097 gammaZ=0 2098 do i=1,3 2099 gammaX=gammaX+gamma(1,i,i,1)+gamma(1,i,1,i)+gamma(1,1,i,i) 2100 gammaY=gammaY+gamma(2,i,i,2)+gamma(2,i,2,i)+gamma(2,2,i,i) 2101 gammaZ=gammaZ+gamma(3,i,i,3)+gamma(3,i,3,i)+gamma(3,3,i,i) 2102 end do 2103 gammaX=gammaX/15 2104 gammaY=gammaY/15 2105 gammaZ=gammaZ/15 2106 gammatot=dsqrt(gammaX**2+gammaY**2+gammaZ**2) 2107 gammaavg1=gammaX+gammaY+gammaZ 2108 gammaavg2=( gamma(1,1,1,1)+gamma(2,2,2,2)+gamma(3,3,3,3) + gamma(1,1,2,2)+gamma(1,1,3,3)+gamma(2,2,3,3) + gamma(2,2,1,1)+gamma(3,3,1,1)+gamma(3,3,2,2) )/5 2109 2110 if (isel==3) then 2111 write(*,*) "Second hyperpolarizability tensor:" 2112 do jdir=1,3 2113 do kdir=1,3 2114 do ldir=1,3 2115 do idir=1,3 2116 write(*,"(2x,4a,'=',1PE14.5,2x)",advance="no") dirlab(idir),dirlab(jdir),dirlab(kdir),dirlab(ldir),gamma(idir,jdir,kdir,ldir) 2117 if (idir==3) write(*,*) 2118 end do 2119 end do 2120 end do 2121 end do 2122 write(*,"(/,' Gamma_X:',1PE14.5,' Gamma_Y:',1PE14.5,' Gamma_Z:',1PE14.5)") gammaX,gammaY,gammaZ 2123 write(*,"(a,1PE14.5)") " Magnitude of gamma:",gammatot 2124 write(*,"(a,1PE14.5)") " Average of gamma (definition 1):",gammaavg1 2125 write(*,"(a,1PE14.5)") " Average of gamma (definition 2):",gammaavg2 2126 write(*,*) 2127 else if (isel==7) then 2128 write(10,"(i6,6(1PE14.5))") numstat,gammaX,gammaY,gammaZ,gammatot,gammaavg1,gammaavg2 2129 else if (isel==17) then 2130 write(10,"(3f12.6,6(1PE14.5))") freq1,freq2,freq3,gammaX,gammaY,gammaZ,gammatot,gammaavg1,gammaavg2 2131 end if 2132 end do !end cycle freqlist 2133 end do !end cycle the number of states 2134 2135 call walltime(iwalltime2) 2136 write(*,"(' Calculation took up wall clock time',i10,'s')") iwalltime2-iwalltime1 2137 if (isel==7.or.isel==17) then 2138 close(10) 2139 if (isel==7) then 2140 write(*,*) "Done! The result has been outputted to gamma_n.txt in current folder" 2141 write(*,*) "The correspondence between columns and information in this file is as follows" 2142 write(*,*) "Column 1: The number of states in consideration" 2143 write(*,*) "Column 2: Gamma_X" 2144 write(*,*) "Column 3: Gamma_Y" 2145 write(*,*) "Column 4: Gamma_Z" 2146 write(*,*) "Column 5: Magnitude of gamma" 2147 write(*,*) "Column 6: Average of gamma (definition 1)" 2148 write(*,*) "Column 7: Average of gamma (definition 2)" 2149 else if (isel==17) then 2150 write(*,*) "Done! The result has been outputted to gamma_w.txt in current folder" 2151 write(*,*) "The correspondence between columns and information in this file is as follows" 2152 write(*,*) "Column 1: Frequency of the first external field (w1)" 2153 write(*,*) "Column 2: Frequency of the second external field (w2)" 2154 write(*,*) "Column 3: Frequency of the third external field (w3)" 2155 write(*,*) "Column 4: Gamma_X" 2156 write(*,*) "Column 5: Gamma_Y" 2157 write(*,*) "Column 6: Gamma_Z" 2158 write(*,*) "Column 7: Magnitude of gamma" 2159 write(*,*) "Column 8: Average of gamma (definition 1)" 2160 write(*,*) "Column 9: Average of gamma (definition 2)" 2161 end if 2162 end if 2163 2164! Calculate third hyperpolarizability 2165else if (isel==4) then 2166 write(*,*) "Input w1,w2,w3,w4 for delta(-(w1+w2+w3+w4);w1,w2,w3,w4)" 2167 write(*,*) "e.g. 0.13,0.13,0,-0.13" 2168 write(*,*) "Note: Negative values mean using nm as unit, e.g. -693,0,0" 2169 read(*,*) freq1,freq2,freq3,freq4 2170 if (freq1<0) freq1=1240.7011D0/au2eV/abs(freq1) 2171 if (freq2<0) freq2=1240.7011D0/au2eV/abs(freq2) 2172 if (freq3<0) freq3=1240.7011D0/au2eV/abs(freq3) 2173 if (freq4<0) freq4=1240.7011D0/au2eV/abs(freq4) 2174 freqtot=freq1+freq2+freq3+freq4 2175 tmpw(1:5)=(/ -freqtot,freq1,freq2,freq3,freq4 /) 2176 if (freq1/=0) then 2177 wavlen1=1240.7011D0/(freq1*au2eV) 2178 write(*,"(' Wavelength of w1:',f12.3,' nm')") wavlen1 2179 end if 2180 if (freq2/=0) then 2181 wavlen2=1240.7011D0/(freq2*au2eV) 2182 write(*,"(' Wavelength of w2:',f12.3,' nm')") wavlen2 2183 end if 2184 if (freq3/=0) then 2185 wavlen3=1240.7011D0/(freq3*au2eV) 2186 write(*,"(' Wavelength of w3:',f12.3,' nm')") wavlen3 2187 end if 2188 if (freq4/=0) then 2189 wavlen4=1240.7011D0/(freq4*au2eV) 2190 write(*,"(' Wavelength of w4:',f12.3,' nm')") wavlen4 2191 end if 2192 write(*,*) 2193 write(*,"(' Consider how many states? Should <=',i6)") nstates 2194 read(*,*) numstat 2195 2196 write(*,*) "Please wait patiently..." 2197 call walltime(iwalltime1) 2198 call fullarrange(arrd,120,5) !Generate full arrangement matrix (4!=24 :4) for gamma, each row corresponds to one permutation, e.g. 2341 2199 iprog=0 2200 do idir=1,3 !Cycle direction component 2201 do jdir=1,3 2202 do kdir=1,3 2203 do ldir=1,3 2204 iprog=iprog+1 2205 write(*,"(' Finished:',i4,' /',i4)") iprog,81 2206 do mdir=1,3 2207 2208 delta1=0 2209 delta2=0 2210 delta3=0 2211 tmpdir(1:5)=(/ idir,jdir,kdir,ldir,mdir /) 2212 do iper=1,120 !Do permutation, arrd(1,:)=1,2,3,4,5 2213 dir1=tmpdir(arrd(iper,1)) 2214 dir2=tmpdir(arrd(iper,2)) 2215 dir3=tmpdir(arrd(iper,3)) 2216 dir4=tmpdir(arrd(iper,4)) 2217 dir5=tmpdir(arrd(iper,5)) 2218 w0=tmpw(arrd(iper,1)) 2219 w1=tmpw(arrd(iper,2)) 2220 w2=tmpw(arrd(iper,3)) 2221 w3=tmpw(arrd(iper,4)) 2222 w4=tmpw(arrd(iper,5)) 2223nthreads=getNThreads() 2224!$OMP PARALLEL SHARED(delta1,delta2,delta3) PRIVATE(istat,jstat,kstat,lstat,t1c,t2c,t3c,p1,p2,p3,p4,d1t,d2t,d3t) NUM_THREADS(nthreads) 2225 d1t=0 2226 d2t=0 2227 d3t=0 2228!$OMP DO schedule(dynamic) 2229 do istat=1,numstat 2230 do jstat=1,numstat 2231 do kstat=1,numstat 2232 !Delta 1 2233 do lstat=1,numstat 2234 t1c=trandip(istat,jstat,dir2) 2235 if (istat==jstat) t1c=t1c-trandip(0,0,dir2) 2236 t2c=trandip(jstat,kstat,dir3) 2237 if (jstat==kstat) t2c=t2c-trandip(0,0,dir3) 2238 t3c=trandip(kstat,lstat,dir4) 2239 if (kstat==lstat) t3c=t3c-trandip(0,0,dir4) 2240 p1=trandip(0,istat,dir1)*t1c*t2c*t3c*trandip(lstat,0,dir5) 2241 p2=(excene(istat)+w0)*(excene(jstat)+w0+w1)*(excene(kstat)-w3-w4)*(excene(lstat)-w4) 2242 d1t=d1t+p1/p2 2243 end do 2244 !Delta 2 2245 t3c=trandip(jstat,kstat,dir4) 2246 if (jstat==kstat) t3c=t3c-trandip(0,0,dir4) 2247 p1=trandip(0,istat,dir1)*trandip(istat,0,dir2)*trandip(0,jstat,dir3)*t3c*trandip(kstat,0,dir5) 2248 p2=(excene(jstat)+w2)*(excene(kstat)-w4) 2249 p3=1/(excene(istat)+w0)+1/(excene(istat)-w1) 2250 p4=1/(excene(jstat)-w3-w4)+1/(excene(kstat)+w2+w3) 2251 d2t=d2t+p1/p2*p3*p4 2252 !Delta 3, p1 is identical to delta 2 counterpart 2253 p2=(excene(istat)+w0)*(excene(istat)-w1) 2254 p3=1/(excene(jstat)-w3-w4)/(excene(kstat)-w4) 2255 p4=1/(excene(jstat)+w2)/(excene(kstat)+w2+w3) 2256 d3t=d3t+p1/p2*(p3+p4) 2257 end do 2258 end do 2259 end do 2260!$OMP END DO 2261!$OMP CRITICAL 2262 delta1=delta1+d1t 2263 delta2=delta2+d2t 2264 delta3=delta3+d3t 2265!$OMP END CRITICAL 2266!$OMP END PARALLEL 2267 end do 2268 delta(idir,jdir,kdir,ldir,mdir)=delta1-delta2/2-delta3/2 2269 2270 end do 2271 end do 2272 end do 2273 end do 2274 end do 2275 2276 write(*,*) 2277 write(*,*) "Third hyperpolarizability tensor:" 2278 do jdir=1,3 2279 do kdir=1,3 2280 do ldir=1,3 2281 do mdir=1,3 2282 do idir=1,3 2283 write(*,"(2x,5a,'=',1PE14.5,2x)",advance="no") & 2284 dirlab(idir),dirlab(jdir),dirlab(kdir),dirlab(ldir),dirlab(mdir),delta(idir,jdir,kdir,ldir,mdir) 2285 if (idir==3) write(*,*) 2286 end do 2287 end do 2288 end do 2289 end do 2290 end do 2291 call walltime(iwalltime2) 2292 write(*,"(' Calculation took up wall clock time',i10,'s')") iwalltime2-iwalltime1 2293 2294end if 2295end do 2296end subroutine 2297 2298 2299!!---------- Calculate average bond length between two elements and average coordinate number 2300subroutine atmavgdist 2301use defvar 2302use util 2303implicit real*8 (a-h,o-z) 2304character elesel1*2,elesel2*2,selectyn 2305if (all(a%name==a(1)%name)) then !Only contain one kind of atom 2306 elesel1=a(1)%name 2307 elesel2=a(1)%name 2308else 2309 write(*,*) "Input two elements for which their average bond length will be calculated" 2310 write(*,*) "For example, B,Al" 2311 read(*,*) elesel1,elesel2 2312 call lc2uc(elesel1(1:1)) 2313 call uc2lc(elesel1(2:2)) 2314 call lc2uc(elesel2(1:1)) 2315 call uc2lc(elesel2(2:2)) 2316end if 2317write(*,*) "Input distance cutoff in Angstrom, e.g. 3.2" 2318read(*,*) discrit 2319discrit=discrit/b2a 2320 2321call gendistmat 2322avgdist=0 2323iwithin=0 2324distmax=0 2325distmin=1000000 2326do iatm=1,ncenter 2327 do jatm=iatm+1,ncenter 2328 if ((a(iatm)%name==elesel1.and.a(jatm)%name==elesel2).or.(a(jatm)%name==elesel1.and.a(iatm)%name==elesel2)) then 2329 dist=distmat(iatm,jatm) 2330 if (dist<=discrit) then 2331 iwithin=iwithin+1 2332 write(*,"(i6,'# ',i6,'(',a') --',i6,'(',a,') Length:',f12.6,' Angstrom')") iwithin,iatm,a(iatm)%name,jatm,a(jatm)%name,dist*b2a 2333 avgdist=avgdist+dist 2334 if ((iatm==1.and.jatm==iatm+1).or.dist>distmax) then 2335 distmax=dist 2336 idistmax1=iatm 2337 idistmax2=jatm 2338 end if 2339 if ((iatm==1.and.jatm==iatm+1).or.dist<distmin) then 2340 distmin=dist 2341 idistmin1=iatm 2342 idistmin2=jatm 2343 end if 2344 end if 2345 end if 2346 end do 2347end do 2348if (iwithin>0) then 2349 avgdist=avgdist/iwithin 2350 write(*,"(' Average bond length between ',a,1x,a,' is',f12.6,' Angstrom')") elesel1,elesel2,avgdist*b2a 2351 write(*,"(' Minimum length is',f12.6,' Angstrom, between',i6,'(',a,') and',i6,'(',a,')')") distmin*b2a,idistmin1,elesel1,idistmin2,elesel2 2352 write(*,"(' Maximum length is',f12.6,' Angstrom, between',i6,'(',a,') and',i6,'(',a,')')") distmax*b2a,idistmax1,elesel1,idistmax2,elesel2 2353else 2354 write(*,*) "No bond that satisfied your criterion is found" 2355 return 2356end if 2357write(*,*) 2358write(*,*) "If also calculate average coordinate number? (y/n)" 2359read(*,*) selectyn 2360if (selectyn=='y'.or.selectyn=='Y') then 2361 ncoordtot=0 2362 ntmp=0 2363 do iatm=1,ncenter 2364 if (a(iatm)%name/=elesel1) cycle 2365 ncoord=0 2366 do jatm=1,ncenter 2367 if (iatm==jatm.or.a(jatm)%name/=elesel2) cycle 2368 if (distmat(iatm,jatm)<=discrit) ncoord=ncoord+1 2369 end do 2370 write(*,"(' The coordinate number of',i6,'(',a,') due to ',a,' - ',a,' bond:',i5)") iatm,elesel1,elesel1,elesel2,ncoord 2371 ncoordtot=ncoordtot+ncoord 2372 ntmp=ntmp+1 2373 end do 2374 write(*,"(/,' The average coordinate number of ',a,' due to ',a,' - ',a,' bond:',f10.5)") elesel1,elesel1,elesel2,dfloat(ncoordtot)/ntmp 2375end if 2376end subroutine 2377 2378 2379 2380!!!------- Calculate electric/magnetic/velocity... integral between orbitals 2381subroutine outorbint 2382use defvar 2383use util 2384implicit real*8 (a-h,o-z) 2385real*8,allocatable :: GTFint(:),GTFvecint(:,:) 2386real*8 vecint(3),vecinttmp(3),intval 2387write(*,*) "Output which kind of integral between orbitals?" 2388write(*,*) "1: Electric dipole moment 2: Magnetic dipole moment 3: Velocity" 2389write(*,*) "4: Kinetic energy 5: Overlap" 2390read(*,*) itype 2391if (wfntype==0.or.wfntype==1.or.wfntype==2) then 2392 write(*,*) "Output the integrals between which orbitals?" 2393 if (wfntype==0.or.wfntype==1.or.wfntype==2) then 2394 write(*,*) "1 Between all occupied orbitals" 2395 write(*,*) "2 Between all occupied and all unoccupied orbitals" 2396 end if 2397 write(*,*) "3 Between all orbitals" 2398 write(*,*) "4 Between the same orbitals" 2399 write(*,*) "5 Between specific range of orbitals" 2400 write(*,*) "6 Between two orbitals" 2401 read(*,*) irange 2402end if 2403if (irange==5) then 2404 write(*,*) "Input orbital range of the first index, e.g. 25,30" 2405 read(*,*) ibeg,iend 2406 write(*,*) "Input orbital range of the second index, e.g. 25,30" 2407 read(*,*) jbeg,jend 2408else if (irange==6) then 2409100 write(*,*) "Input index for the two orbitals, e.g. 144,340" 2410 write(*,*) "Input 0,0 can exit" 2411 read(*,*) iMOsel,jMOsel 2412 if (iMOsel==0.and.jMOsel==0) return 2413end if 2414 2415call walltime(iwalltime1) 2416 2417nsize=nprims*(nprims+1)/2 2418if (itype==1.or.itype==2.or.itype==3) then 2419 allocate(GTFvecint(3,nsize)) 2420else 2421 allocate(GTFint(nsize)) 2422end if 2423if (itype==1) then 2424 call genGTFDmat(GTFvecint,nsize) 2425else if (itype==2) then 2426 call genGTFMmat(GTFvecint,nsize) !Notice that this only generate lower triangular matrix, (i,j)=-(j,i) 2427else if (itype==3) then 2428 call genGTFVelmat(GTFvecint,nsize) 2429else if (itype==4) then 2430 call genGTFTmat(GTFint,nsize) 2431else if (itype==5) then 2432 call genGTFSmat(GTFint,nsize) 2433end if 2434 2435if (irange/=6) open(10,file="orbint.txt",status="replace") 2436do imo=1,nmo 2437 do jmo=1,nmo 2438 if (irange==1) then 2439 if (MOocc(imo)==0D0.or.MOocc(jmo)==0D0) cycle 2440 else if (irange==2) then 2441 if (MOocc(imo)==0D0.or.MOocc(jmo)/=0D0) cycle 2442 else if (irange==4) then 2443 if (imo/=jmo) cycle 2444 else if (irange==5) then 2445 if (imo<ibeg.or.imo>iend.or.jmo<jbeg.or.jmo>jend) cycle 2446 else if (irange==6) then 2447 if (imo/=iMOsel.or.jmo/=jMOsel) cycle 2448 end if 2449 2450 if (itype==1.or.itype==2.or.itype==3) then !Vector integral 2451 vecint=0D0 2452nthreads=getNThreads() 2453!$OMP PARALLEL SHARED(vecint) PRIVATE(iGTF,jGTF,ides,vecinttmp) NUM_THREADS(nthreads) 2454 vecinttmp=0D0 2455!$OMP DO schedule(dynamic) 2456 do iGTF=1,nprims 2457 do jGTF=1,nprims 2458 if (iGTF>=jGTF) then 2459 ides=iGTF*(iGTF-1)/2+jGTF 2460 else 2461 ides=jGTF*(jGTF-1)/2+iGTF 2462 end if 2463 if ((itype==2.or.itype==3).and.iGTF>jGTF) then !Magnetic and velocity operators are Hermitean, so inverse sign 2464 vecinttmp=vecinttmp-co(imo,iGTF)*co(jmo,jGTF)*GTFvecint(:,ides) 2465 else 2466 vecinttmp=vecinttmp+co(imo,iGTF)*co(jmo,jGTF)*GTFvecint(:,ides) 2467 end if 2468 end do 2469 end do 2470!$OMP END DO 2471!$OMP CRITICAL 2472 vecint=vecint+vecinttmp 2473!$OMP END CRITICAL 2474!$OMP END PARALLEL 2475 if (irange==6) then 2476 write(*,"(' Result:',3f18.10)") vecint(:) 2477 else 2478 write(10,"(2i8,4x,3f18.10)") imo,jmo,vecint(:) 2479 end if 2480 else !Scalar integral 2481 intval=0D0 2482 do iGTF=1,nprims 2483 do jGTF=1,nprims 2484 if (iGTF>=jGTF) then 2485 ides=iGTF*(iGTF-1)/2+jGTF 2486 else 2487 ides=jGTF*(jGTF-1)/2+iGTF 2488 end if 2489 intval=intval+co(imo,iGTF)*co(jmo,jGTF)*GTFint(ides) 2490 end do 2491 end do 2492 if (irange==6) then 2493 write(*,"(' Result:',f18.10)") intval 2494 else 2495 write(10,"(2i8,4x,f18.10)") imo,jmo,intval 2496 end if 2497 end if 2498 end do 2499 if (irange/=6) write(*,"(' Finished',i7,' /',i7)") imo,nmo 2500end do 2501if (irange==6) then 2502 if (allocated(GTFvecint)) deallocate(GTFvecint) 2503 if (allocated(GTFint)) deallocate(GTFint) 2504 write(*,*) 2505 goto 100 2506else 2507 close(10) 2508 call walltime(iwalltime2) 2509 write(*,"(' Done! Calculation took up wall clock time',i10,'s',/)") iwalltime2-iwalltime1 2510 write(*,*) "The integrals has been outputted to orbint.txt in current folder" 2511 if (itype==1.or.itype==2.or.itype==3) then 2512 write(*,"(a)") " The first and the second columns correspond to orbital indices. The last three columns correspond to the integral in X/Y/Z (a.u.)" 2513 else 2514 write(*,"(a)") " The first and the second columns correspond to orbital indices. The last column corresponds to the integral value (a.u.)" 2515 end if 2516end if 2517end subroutine 2518 2519 2520 2521!!----------- Calculate center, first and second moments of a real space function 2522subroutine funcmoment 2523use defvar 2524use util 2525use function 2526implicit real*8 (a-h,o-z) 2527real*8 intval,moment1(3),moment2(3,3),moment2nuc(3,3),funcval(radpot*sphpot),beckeweigrid(radpot*sphpot),eigvecmat(3,3),eigval(3) 2528type(content) gridatm(radpot*sphpot),gridatmorg(radpot*sphpot) 2529character selectyn 2530ifunc=1 2531cenx=0 2532ceny=0 2533cenz=0 2534do while(.true.) 2535 write(*,*) "0 Return" 2536 write(*,*) "1 Calculate the first and second moments of the function" 2537 write(*,*) "2 Calculate the center and integral of the function" 2538 write(*,"(a,i5)") " 3 Select the function to be studied, current:",ifunc 2539 write(*,"(a,3f11.5,' Ang')") " 4 Set the center for option 1, current:",cenx*b2a,ceny*b2a,cenz*b2a 2540 read(*,*) isel 2541 2542 if (isel==0) then 2543 return 2544 else if (isel==3) then 2545 call selfunc_interface(ifunc) 2546 cycle 2547 else if (isel==4) then 2548 write(*,*) "Input X,Y,Z of the center in Angstrom, e.g. 2.0,0,1.5" 2549 read(*,*) cenx,ceny,cenz 2550 cenx=cenx/b2a !To Bohr 2551 ceny=ceny/b2a 2552 cenz=cenz/b2a 2553 cycle 2554 end if 2555 2556 write(*,"(' Radial points:',i5,' Angular points:',i5,' Total:',i10,' per center')") radpot,sphpot,radpot*sphpot 2557 call gen1cintgrid(gridatmorg,iradcut) 2558 2559 call walltime(iwalltime1) 2560 CALL CPU_TIME(time_begin) 2561 2562 intval=0 2563 moment1=0 2564 moment2=0 2565 realcenx=0 2566 realceny=0 2567 realcenz=0 2568 do iatm=1,ncenter 2569 write(*,"(' Processing center',i6,'(',a2,') /',i6)") iatm,a(iatm)%name,ncenter 2570 gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule 2571 gridatm%y=gridatmorg%y+a(iatm)%y 2572 gridatm%z=gridatmorg%z+a(iatm)%z 2573nthreads=getNThreads() 2574!$OMP parallel do shared(funcval) private(i) num_threads(nthreads) 2575 do i=1+iradcut*sphpot,radpot*sphpot 2576 funcval(i)=calcfuncall(ifunc,gridatm(i)%x,gridatm(i)%y,gridatm(i)%z) 2577 end do 2578!$OMP end parallel do 2579 call gen1cbeckewei(iatm,iradcut,gridatm,beckeweigrid) 2580 do i=1+iradcut*sphpot,radpot*sphpot 2581 tmpval=funcval(i)*gridatmorg(i)%value*beckeweigrid(i) 2582 xtmp=gridatm(i)%x-cenx 2583 ytmp=gridatm(i)%y-ceny 2584 ztmp=gridatm(i)%z-cenz 2585 intval=intval+tmpval 2586 moment1(1)=moment1(1)+xtmp*tmpval 2587 moment1(2)=moment1(2)+ytmp*tmpval 2588 moment1(3)=moment1(3)+ztmp*tmpval 2589 moment2(1,1)=moment2(1,1)+xtmp*xtmp*tmpval 2590 moment2(2,2)=moment2(2,2)+ytmp*ytmp*tmpval 2591 moment2(3,3)=moment2(3,3)+ztmp*ztmp*tmpval 2592 moment2(1,2)=moment2(1,2)+xtmp*ytmp*tmpval 2593 moment2(2,3)=moment2(2,3)+ytmp*ztmp*tmpval 2594 moment2(1,3)=moment2(1,3)+xtmp*ztmp*tmpval 2595 realcenx=realcenx+gridatm(i)%x*tmpval 2596 realceny=realceny+gridatm(i)%y*tmpval 2597 realcenz=realcenz+gridatm(i)%z*tmpval 2598 end do 2599 end do 2600 CALL CPU_TIME(time_end) 2601 call walltime(iwalltime2) 2602 write(*,"(' Calculation took up CPU time',f12.2,'s, wall clock time',i10,'s',/)") time_end-time_begin,iwalltime2-iwalltime1 2603 2604 if (isel==1) then 2605 moment2(3,1)=moment2(1,3) 2606 moment2(2,1)=moment2(1,2) 2607 moment2(3,2)=moment2(2,3) 2608 write(*,*) "Note: All units below are in a.u." 2609 write(*,"(/,' Integral:',1PE16.8,/)") intval 2610 write(*,"(' The first moment:')") 2611 write(*,"(' X= ',1PE16.8,' Y= ',1PE16.8,' Z= ',1PE16.8)") moment1 2612 write(*,"(' Norm= ',1PE16.8,/)") sum(moment1**2) 2613 write(*,"(' The second moment:')") 2614 write(*,"(' XX=',1PE16.8,' XY=',1PE16.8,' XZ=',1PE16.8)") moment2(1,:) 2615 write(*,"(' YX=',1PE16.8,' YY=',1PE16.8,' YZ=',1PE16.8)") moment2(2,:) 2616 write(*,"(' ZX=',1PE16.8,' ZY=',1PE16.8,' ZZ=',1PE16.8)") moment2(3,:) 2617 2618 call diagmat(moment2,eigvecmat,eigval,300,1D-10) 2619 call sortr8(eigval) 2620 write(*,"(a,3(1PE16.8))") ' Eigenvalues:',eigval 2621 write(*,"(' Anisotropy:',1PE16.8,/)") eigval(3)-(eigval(1)+eigval(2))/2D0 2622 write(*,"(' Radius of gyration:',1PE16.8,/)") dsqrt((moment2(1,1)+moment2(2,2)+moment2(3,3))/intval) 2623 2624 if (ifunc==1) then 2625 moment2nuc=0 2626 do iatm=1,ncenter 2627 xtmp=a(iatm)%x-cenx 2628 ytmp=a(iatm)%y-ceny 2629 ztmp=a(iatm)%z-cenz 2630 tmpval=a(iatm)%charge 2631 moment2nuc(1,1)=moment2nuc(1,1)+xtmp*xtmp*tmpval 2632 moment2nuc(2,2)=moment2nuc(2,2)+ytmp*ytmp*tmpval 2633 moment2nuc(3,3)=moment2nuc(3,3)+ztmp*ztmp*tmpval 2634 moment2nuc(1,2)=moment2nuc(1,2)+xtmp*ytmp*tmpval 2635 moment2nuc(2,3)=moment2nuc(2,3)+ytmp*ztmp*tmpval 2636 moment2nuc(1,3)=moment2nuc(1,3)+xtmp*ztmp*tmpval 2637 end do 2638 moment2nuc(3,1)=moment2nuc(1,3) 2639 moment2nuc(2,1)=moment2nuc(1,2) 2640 moment2nuc(3,2)=moment2nuc(2,3) 2641 write(*,*) 2642 write(*,"(' The quadrupole moment of nuclear charges:')") 2643 write(*,"(' XX=',f16.8,' XY=',f16.8,' XZ=',f16.8)") moment2nuc(1,:) 2644 write(*,"(' YX=',f16.8,' YY=',f16.8,' YZ=',f16.8)") moment2nuc(2,:) 2645 write(*,"(' ZX=',f16.8,' ZY=',f16.8,' ZZ=',f16.8)") moment2nuc(3,:) 2646 write(*,*) 2647 write(*,"(' The quadrupole moment of the system:')") 2648 write(*,"(' XX=',f16.8,' XY=',f16.8,' XZ=',f16.8)") moment2nuc(1,:)-moment2(1,:) 2649 write(*,"(' YX=',f16.8,' YY=',f16.8,' YZ=',f16.8)") moment2nuc(2,:)-moment2(2,:) 2650 write(*,"(' ZX=',f16.8,' ZY=',f16.8,' ZZ=',f16.8)") moment2nuc(3,:)-moment2(3,:) 2651 end if 2652 2653 else if (isel==2) then 2654 write(*,"(' Integral:',1PE16.8,/)") intval 2655 realcenx=realcenx/intval 2656 realceny=realceny/intval 2657 realcenz=realcenz/intval 2658 write(*,"(' The center of the function:')") 2659 write(*,"(' X=',f16.8,' Y=',f16.8,' Z=',f16.8,' Angstrom',/)") realcenx*b2a,realceny*b2a,realcenz*b2a 2660 write(*,*) "Use this center for succeeding calculations? (y/n)" 2661 read(*,*) selectyn 2662 if (selectyn=='y') then 2663 cenx=realcenx 2664 ceny=realceny 2665 cenz=realcenz 2666 end if 2667 end if 2668end do 2669end subroutine 2670 2671 2672 2673!!----------- Calculate energy index (EI) or bond polarity index (BPI) 2674!!J. Phys. Chem., 94, 5602-5607 and J. Phys. Chem.,96, 157-164 2675subroutine calcEIBPI 2676use defvar 2677implicit real*8 (a-h,o-z) 2678if (wfntype==3.or.wfntype==4) then 2679 write(*,*) "Error: Post-HF wavefunction is not supported yet!" 2680 return 2681end if 2682ninnerele=0 2683do i=1,ncenter 2684 if (int(a(i)%charge)/=a(i)%index) cycle !For the atom used ECP, skip it 2685 if (a(i)%index>2.and.a(i)%index<=10) ninnerele=ninnerele+2 2686 if (a(i)%index>10.and.a(i)%index<=18) ninnerele=ninnerele+10 2687 if (a(i)%index>18.and.a(i)%index<=36) ninnerele=ninnerele+18 2688 if (a(i)%index>36.and.a(i)%index<=54) ninnerele=ninnerele+36 2689 if (a(i)%index>54.and.a(i)%index<=86) ninnerele=ninnerele+54 2690 if (a(i)%index>86) ninnerele=ninnerele+86 2691end do 2692write(*,"(' The number of inner electrons is assumed to be',i5,/)") ninnerele 2693do while(.true.) 2694 write(*,*) "Calculate EI index for which atom? e.g. 5" 2695 write(*,*) "Input 0 can exit" 2696 read(*,*) iatm 2697 if (iatm==0) exit 2698 val1=0D0 2699 val2=0D0 2700 do imo=ninnerele/2+1,nbasis 2701 if (MOocc(imo)==0D0) exit 2702 compos=0 2703 do ibas=basstart(iatm),basend(iatm) 2704 compos=compos+CObasa(ibas,imo)**2 2705 do jbas=1,nbasis 2706 if (jbas==ibas) cycle 2707 compos=compos+CObasa(ibas,imo)*CObasa(jbas,imo)*Sbas(ibas,jbas) 2708 end do 2709 end do 2710! write(*,"(i5,3f12.6)") imo,MOocc(imo),MOene(imo),compos 2711 val1=val1+MOocc(imo)*MOene(imo)*compos 2712 val2=val2+MOocc(imo)*compos 2713 end do 2714 if (wfntype==1) then !beta part 2715 do imo=nbasis+ninnerele/2+1,nmo 2716 if (MOocc(imo)==0D0) exit 2717 compos=0 2718 do ibas=basstart(iatm),basend(iatm) 2719 compos=compos+CObasb(ibas,imo-nbasis)**2 2720 do jbas=1,nbasis 2721 if (jbas==ibas) cycle 2722 compos=compos+CObasb(ibas,imo-nbasis)*CObasb(jbas,imo-nbasis)*Sbas(ibas,jbas) 2723 end do 2724 end do 2725! write(*,"(i5,3f12.6)") imo,MOocc(imo),MOene(imo),compos 2726 val1=val1+MOocc(imo)*MOene(imo)*compos 2727 val2=val2+MOocc(imo)*compos 2728 end do 2729 end if 2730 write(*,"(' The numerator: ',f12.6,' a.u.')") val1 2731 write(*,"(' The denominator:',f12.6,' a.u.')") val2 !Corresponding to Mulliken occupation number of this atom in valence MOs 2732 write(*,"(' The EI index: ',f12.6,' a.u.',/)") val1/val2 2733end do 2734end subroutine 2735 2736 2737 2738!!--------- Calculate electronic coupling matrix element by FCD method 2739!FCDtest.fch :1-16,17-31,32-47 2740subroutine FCD 2741use defvar 2742use util 2743implicit real*8 (a-h,o-z) 2744integer :: itranstype=1 !1=hole transfer, 2=electron transfer 2745integer :: iFCDorbtype=1 !For UHF, 1=alpha 2=beta 2746integer :: imethod=1 !1=FCD 2=GMH 2747integer,allocatable :: donor(:),bridge(:),acceptor(:),donorMO(:),bridgeMO(:),acceptorMO(:) 2748real*8 :: compthresFCD=0.89D0,miu1(3),miu2(3),miu12(3) 2749real*8 :: MOcomp(nbasis,3) !Composition of donor,bridge,acceptor in each MO 2750character c2000tmp*2000,c80tmp*80 2751if (wfntype/=0) then 2752 write(*,*) "Error: This module is only applicable to RHF/RKS wavefunction" 2753 return 2754end if 2755ndonoratm=0 2756nbridgeatm=0 2757nacceptoratm=0 2758ndonorMO=1 2759nbridgeMO=0 2760nacceptorMO=1 2761 2762!For debug, GTG base stack example 2763ndonoratm=16 2764nbridgeatm=15 2765nacceptoratm=16 2766allocate(donor(ndonoratm)) 2767allocate(bridge(nbridgeatm)) 2768allocate(acceptor(nacceptoratm)) 2769do i=1,16 2770 donor(i)=i 2771 if (i/=16) bridge(i)=i+16 2772 acceptor(i)=i+31 2773end do 2774ndonorMO=1 2775nbridgeMO=1 2776nacceptorMO=1 2777!For debug, GG base stack example 2778! ndonoratm=16 2779! nbridgeatm=0 2780! nacceptoratm=16 2781! allocate(donor(ndonoratm)) 2782! allocate(bridge(nbridgeatm)) 2783! allocate(acceptor(nacceptoratm)) 2784! do i=1,16 2785! donor(i)=i 2786! acceptor(i)=i+16 2787! end do 2788! ndonorMO=1 2789! nbridgeMO=0 2790! nacceptorMO=1 2791 2792write(*,*) "Citation of FCD method: J. Chem. Phys., 117, 5607 (2002)" 2793130 do while(.true.) 2794 write(*,*) 2795 write(*,*) " ========== Fragment charge difference method (FCD) ==========" 2796 if (itranstype==1) write(*,*) "-1 Show information of some highest occupied MOs" 2797 if (itranstype==2) write(*,*) "-1 Show information of some lowest unoccupied MOs" 2798 write(*,*) "0 Start calculation" 2799 write(*,"(a,i4)") " 1 Define donor fragment, number of atoms:",ndonoratm 2800 write(*,"(a,i4)") " 2 Define bridge fragment (if any), number of atoms:",nbridgeatm 2801 write(*,"(a,i4)") " 3 Define acceptor fragment, number of atoms:",nacceptoratm 2802 write(*,"(a,f5.1,a)") " 4 Set composition threshold for determining attribution, current:",compthresFCD*100,"%" 2803 write(*,"(a,3i2)") " 5 The number of MOs attributed to donor, bridge and acceptor, current:",ndonorMO,nbridgeMO,nacceptorMO 2804 if (itranstype==1) write(*,*) "6 Switch type of electron transfer, current: Hole transfer" 2805 if (itranstype==2) write(*,*) "6 Switch type of electron transfer, current: Electron transfer" 2806! if (wfntype==1.and.iFCDorbtype==1) write(*,*) "7 Switch molecular orbital type, current: Alpha" 2807! if (wfntype==1.and.iFCDorbtype==2) write(*,*) "7 Switch molecular orbital type, current: Beta" 2808 if (imethod==1) write(*,*) "8 Switch method between FCD and GMH, current: FCD" 2809 if (imethod==2) write(*,*) "8 Switch method between FCD and GMH, current: GMH" 2810 read(*,*) isel 2811 if (isel==-1) then 2812 if (itranstype==1) write(*,*) "Show how many highest occupied MOs? e.g. 30" 2813 if (itranstype==2) write(*,*) "Show how many lowest unoccupied MOs? e.g. 30" 2814 write(*,*) "If press ENTER directly, 50 MOs will be shown" 2815 read(*,"(a)") c80tmp 2816 if (c80tmp==" ") then 2817 nshowMO=50 2818 else 2819 read(c80tmp,*) nshowMO 2820 end if 2821 exit 2822 else if (isel==0) then 2823 if (ndonoratm==0.or.nacceptoratm==0) then 2824 write(*,*) "Error: You should define donor and acceptor first!" 2825 cycle 2826 end if 2827 exit 2828 else if (isel==1) then 2829 write(*,*) "Input the atom indices constituting the donor" 2830 write(*,*) "e.g. 1,3-5,9" 2831 read(*,"(a)") c2000tmp 2832 if (allocated(donor)) deallocate(donor) 2833 call str2arr(c2000tmp,ndonoratm) 2834 allocate(donor(ndonoratm)) 2835 call str2arr(c2000tmp,ndonoratm,donor) 2836 write(*,*) "Done!" 2837 else if (isel==2) then 2838 write(*,*) "Input the atom indices constituting the bridge" 2839 write(*,*) "e.g. 1,3-5,9" 2840 write(*,*) "(If directly press ENTER then bridge will not be defined)" 2841 read(*,"(a)") c2000tmp 2842 if (c2000tmp==" ") then 2843 nbridgeatm=0 2844 cycle 2845 end if 2846 if (allocated(bridge)) deallocate(bridge) 2847 call str2arr(c2000tmp,nbridgeatm) 2848 allocate(bridge(nbridgeatm)) 2849 call str2arr(c2000tmp,nbridgeatm,bridge) 2850 write(*,*) "Done!" 2851 else if (isel==3) then 2852 write(*,*) "Input the atom indices constituting the acceptor" 2853 write(*,*) "e.g. 1,3-5,9" 2854 read(*,"(a)") c2000tmp 2855 if (allocated(acceptor)) deallocate(acceptor) 2856 call str2arr(c2000tmp,nacceptoratm) 2857 allocate(acceptor(nacceptoratm)) 2858 call str2arr(c2000tmp,nacceptoratm,acceptor) 2859 write(*,*) "Done!" 2860 else if (isel==4) then 2861 write(*,*) "Input composition threshold (%), e.g. 89" 2862 read(*,*) compthresFCD 2863 compthresFCD=compthresFCD/100 2864 else if (isel==5) then 2865 write(*,*) "Input the number of MOs attributed to donor, bridge and acceptor, respectively" 2866 write(*,*) "e.g. 2,2,1" 2867 read(*,*) ndonorMO,nbridgeMO,nacceptorMO 2868 else if (isel==6) then 2869 if (itranstype==1) then 2870 itranstype=2 2871 else if (itranstype==2) then 2872 itranstype=1 2873 end if 2874 else if (isel==7) then 2875 if (iFCDorbtype==1) then 2876 iFCDorbtype=2 2877 else if (iFCDorbtype==2) then 2878 iFCDorbtype=1 2879 end if 2880 else if (isel==8) then 2881 if (imethod==1) then 2882 if (igenDbas==0) then !Haven't calculated dipole moment integral matrix 2883 write(*,"(a)") " Warning: To use GMH method, you should set ""igenDbas"" in settings.ini to 1, & 2884 and then restart Multiwfn, so that dipole moment integral matrix could be generated when loading file" 2885 write(*,"(a)") " Alternatively, now you can press ENTER to reload input file and meantime generate dipole moment integral matrix" 2886 read(*,*) 2887 call dealloall 2888 write(*,*) "Reloading input file..." 2889 igenDbas=1 2890 call readinfile(firstfilename,1) 2891 end if 2892 imethod=2 2893 else if (imethod==2) then 2894 imethod=1 2895 end if 2896 end if 2897end do 2898 2899iHOMO=nelec/2 2900iLUMO=iHOMO+1 2901write(*,"(' HOMO is',i6)") iHOMO 2902write(*,"(' LUMO is',i6)") iLUMO 2903 2904!Determine attribution of complex MOs 2905if (isel==-1) then 2906 write(*,*) " MO Energy(eV) Donor(%) Bridge(%) Acceptor(%)" 2907 icount=0 2908else 2909 allocate(donorMO(ndonorMO),bridgeMO(nbridgeMO),acceptorMO(nacceptorMO)) 2910end if 2911 2912idonor=1 2913ibridge=1 2914iacceptor=1 2915if (itranstype==1) then !hole transfer, use HOMO and below MOs 2916 iMOstart=iHOMO 2917 iMOend=1 2918 iMOstep=-1 2919else !electron transfer, use LUMO and above MOs 2920 iMOstart=iLUMO 2921 iMOend=nmo 2922 iMOstep=1 2923end if 2924do imo=iMOstart,iMOend,iMOstep 2925 compdonor=0 2926 do iatmidx=1,ndonoratm 2927 iatm=donor(iatmidx) 2928 do ibas=basstart(iatm),basend(iatm) 2929 do jbas=1,nbasis 2930 compdonor=compdonor+CObasa(ibas,imo)*CObasa(jbas,imo)*Sbas(ibas,jbas) 2931 end do 2932 end do 2933 end do 2934 compbridge=0 2935 do iatmidx=1,nbridgeatm 2936 iatm=bridge(iatmidx) 2937 do ibas=basstart(iatm),basend(iatm) 2938 do jbas=1,nbasis 2939 compbridge=compbridge+CObasa(ibas,imo)*CObasa(jbas,imo)*Sbas(ibas,jbas) 2940 end do 2941 end do 2942 end do 2943 compacceptor=0 2944 do iatmidx=1,nacceptoratm 2945 iatm=acceptor(iatmidx) 2946 do ibas=basstart(iatm),basend(iatm) 2947 do jbas=1,nbasis 2948 compacceptor=compacceptor+CObasa(ibas,imo)*CObasa(jbas,imo)*Sbas(ibas,jbas) 2949 end do 2950 end do 2951 end do 2952 MOcomp(imo,1)=compdonor 2953 MOcomp(imo,2)=compbridge 2954 MOcomp(imo,3)=compacceptor 2955 if (isel==-1) then !Examine composition of some MOs 2956 write(*,"(i8,f10.3,3f12.2)") imo,MOene(imo)*au2eV,MOcomp(imo,:)*100 2957 icount=icount+1 2958 if (icount==nshowMO) exit !Show up to "nshowMO" MOs 2959 cycle 2960 end if 2961 if (idonor<=ndonorMO.and.compdonor>compthresFCD) then 2962 donorMO(idonor)=imo 2963 idonor=idonor+1 2964 write(*,"(' MO',i4,' (',f9.3,'eV) is attributed to donor with composition of',f5.1'%')") imo,MOene(imo)*au2eV,compdonor*100 2965 else if (iacceptor<=nacceptorMO.and.compacceptor>compthresFCD) then 2966 acceptorMO(iacceptor)=imo 2967 iacceptor=iacceptor+1 2968 write(*,"(' MO',i4,' (',f9.3,'eV) is attributed to acceptor with composition of',f5.1'%')") imo,MOene(imo)*au2eV,compacceptor*100 2969 else if (ibridge<=nbridgeMO.and.compbridge>compthresFCD) then 2970 bridgeMO(ibridge)=imo 2971 ibridge=ibridge+1 2972 write(*,"(' MO',i4,' (',f9.3,'eV) is attributed to bridge with composition of',f5.1'%')") imo,MOene(imo)*au2eV,compbridge*100 2973 end if 2974 if (idonor>ndonorMO.and.ibridge>nbridgeMO.and.iacceptor>nacceptorMO) exit 2975end do 2976 2977if (isel==-1) goto 130 2978 2979write(*,*) 2980write(*,"(' Donor atoms:')") 2981write(*,"(15i5)") donor 2982if (nbridgeatm>0.and.nbridgeMO>0) then 2983 write(*,"(' Bridge atoms:')") 2984 write(*,"(15i5)") bridge 2985end if 2986write(*,"(' Acceptor atoms:')") 2987write(*,"(15i5)") acceptor 2988 2989write(*,*) 2990write(*,*) 2991write(*,*) " ======== Results between donor and acceptor ========" 2992write(*,*) 2993rmsH_DA=0 2994avgE_DA=0 2995do idx=1,ndonorMO 2996 imo=donorMO(idx) 2997 E1=MOene(imo) 2998 q1D=MOcomp(imo,1) 2999 q1A=MOcomp(imo,3) 3000 dq1=q1D-q1A 3001 do jdx=1,nacceptorMO 3002 jmo=acceptorMO(jdx) 3003 E2=MOene(jmo) 3004 enediff=abs(E2-E1) !Always take positive energy difference 3005 q2D=MOcomp(jmo,1) 3006 q2A=MOcomp(jmo,3) 3007 dq2=q2D-q2A 3008 q12_1=0 3009 write(*,"(' **** MO',i4,' (donor) --- MO',i4,' (acceptor) ****')") imo,jmo 3010 write(*,"(' MO energy difference:',f12.4,' eV')") enediff*au2eV 3011 if (imethod==1) then 3012 do iatmidx=1,ndonoratm 3013 iatm=donor(iatmidx) 3014 do ibas=basstart(iatm),basend(iatm) 3015 do jbas=1,nbasis 3016 q12_1=q12_1+( CObasa(ibas,imo)*CObasa(jbas,jmo)+CObasa(ibas,jmo)*CObasa(jbas,imo) )*Sbas(ibas,jbas) 3017 end do 3018 end do 3019 end do 3020 q12_1=q12_1/2 3021 q12_2=0 3022 do iatmidx=1,nacceptoratm 3023 iatm=acceptor(iatmidx) 3024 do ibas=basstart(iatm),basend(iatm) 3025 do jbas=1,nbasis 3026 q12_2=q12_2+( CObasa(ibas,imo)*CObasa(jbas,jmo)+CObasa(ibas,jmo)*CObasa(jbas,imo) )*Sbas(ibas,jbas) 3027 end do 3028 end do 3029 end do 3030 q12_2=q12_2/2 3031 dq12=q12_1-q12_2 3032 H_DA=enediff*abs(dq12)/dsqrt((dq1-dq2)**2+4*dq12**2) 3033 write(*,"(' q1(donor):',f10.6,' q1(acceptor):',f10.6,' delta_q1:',f10.6)") q1D,q1A,dq1 3034 write(*,"(' q2(donor):',f10.6,' q2(acceptor):',f10.6,' delta_q2:',f10.6)") q2D,q2A,dq2 3035 write(*,"(' delta_q1 - delta_q2:',f10.6)") dq1-dq2 3036 write(*,"(' q12(donor):',f10.6,' q12(acceptor):',f10.6,' delta_q12:',f10.6)") q12_1,q12_2,dq12 3037 else if (imethod==2) then 3038 miu1=0 3039 miu2=0 3040 miu12=0 3041 do ibas=1,nbasis !Note that Dbas has already taken the negative sign of electron charge into account 3042 do jbas=1,nbasis 3043 miu1(:)=miu1(:)+CObasa(ibas,imo)*CObasa(jbas,imo)*Dbas(:,ibas,jbas) 3044 miu2(:)=miu2(:)+CObasa(ibas,jmo)*CObasa(jbas,jmo)*Dbas(:,ibas,jbas) 3045 miu12(:)=miu12(:)+( CObasa(ibas,imo)*CObasa(jbas,jmo)+CObasa(ibas,jmo)*CObasa(jbas,imo) )*Dbas(:,ibas,jbas) 3046 end do 3047 end do 3048 miu12=miu12/2 3049 valmiu1=dsqrt(sum(miu1**2)) 3050 valmiu2=dsqrt(sum(miu2**2)) 3051 valmiu12=dsqrt(sum(miu12**2)) 3052 write(*,"(' miu1 (X/Y/Z): ',3f11.6,' Norm:',f11.6,' a.u.')") miu1,valmiu1 3053 write(*,"(' miu2 (X/Y/Z): ',3f11.6,' Norm:',f11.6,' a.u.')") miu2,valmiu2 3054 write(*,"(' miu1-miu2 (X/Y/Z):',3f11.6,' Norm:',f11.6,' a.u.')") miu1-miu2,dsqrt(sum((miu1-miu2)**2)) 3055 write(*,"(' miu12 (X/Y/Z): ',3f11.6,' Norm:',f11.6,' a.u.')") miu12,valmiu12 3056 H_DA=enediff*valmiu12/dsqrt(sum((miu1-miu2)**2)+4*valmiu12**2) 3057 end if 3058 write(*,"(' H_DA:',f14.8,' eV')") H_DA*au2eV 3059 write(*,*) 3060 rmsH_DA=rmsH_DA+H_DA**2 3061 avgE_DA=avgE_DA+enediff 3062 end do 3063end do 3064rmsH_DA=dsqrt(rmsH_DA/ndonorMO/nacceptorMO) 3065avgE_DA=avgE_DA/(ndonorMO*nacceptorMO) 3066write(*,"(' Average MO energy difference E_DA:',f14.8,'eV')") avgE_DA*au2eV 3067write(*,"(' Root mean square H_DA:',f14.8,'eV')") rmsH_DA*au2eV 3068 3069if (nbridgeatm>0.and.nbridgeMO>0) then 3070 write(*,*) 3071 write(*,*) 3072 write(*,*) " ======== Results between donor and bridge ========" 3073 write(*,*) 3074 rmsH_DB=0 3075 avgE_DB=0 3076 do idx=1,ndonorMO 3077 imo=donorMO(idx) 3078 E1=MOene(imo) 3079 q1D=MOcomp(imo,1) 3080 q1B=MOcomp(imo,2) 3081 dq1=q1D-q1B 3082 do jdx=1,nbridgeMO 3083 jmo=bridgeMO(jdx) 3084 E2=MOene(jmo) 3085 enediff=abs(E2-E1) !Always take positive energy difference 3086 q2D=MOcomp(jmo,1) 3087 q2B=MOcomp(jmo,2) 3088 dq2=q2D-q2B 3089 q12_1=0 3090 write(*,"(' **** MO',i4,' (donor) --- MO',i4,' (bridge) ****')") imo,jmo 3091 write(*,"(' MO energy difference:',f12.4,' eV')") enediff*au2eV 3092 if (imethod==1) then 3093 do iatmidx=1,ndonoratm 3094 iatm=donor(iatmidx) 3095 do ibas=basstart(iatm),basend(iatm) 3096 do jbas=1,nbasis 3097 q12_1=q12_1+( CObasa(ibas,imo)*CObasa(jbas,jmo)+CObasa(ibas,jmo)*CObasa(jbas,imo) )*Sbas(ibas,jbas) 3098 end do 3099 end do 3100 end do 3101 q12_1=q12_1/2 3102 q12_2=0 3103 do iatmidx=1,nbridgeatm 3104 iatm=bridge(iatmidx) 3105 do ibas=basstart(iatm),basend(iatm) 3106 do jbas=1,nbasis 3107 q12_2=q12_2+( CObasa(ibas,imo)*CObasa(jbas,jmo)+CObasa(ibas,jmo)*CObasa(jbas,imo) )*Sbas(ibas,jbas) 3108 end do 3109 end do 3110 end do 3111 q12_2=q12_2/2 3112 dq12=q12_1-q12_2 3113 H_DB=enediff*abs(dq12)/dsqrt((dq1-dq2)**2+4*dq12**2) 3114 write(*,"(' q1(donor):',f10.6,' q1(bridge):',f10.6,' delta_q1:',f10.6)") q1D,q1B,dq1 3115 write(*,"(' q2(donor):',f10.6,' q2(bridge):',f10.6,' delta_q2:',f10.6)") q2D,q2B,dq2 3116 write(*,"(' delta_q1 - delta_q2:',f10.6)") dq1-dq2 3117 write(*,"(' q12(donor):',f10.6,' q12(bridge):',f10.6,' delta_q12:',f10.6)") q12_1,q12_2,dq12 3118 else if (imethod==2) then 3119 miu1=0 3120 miu2=0 3121 miu12=0 3122 do ibas=1,nbasis !Note that Dbas has already taken the negative sign of electron charge into account 3123 do jbas=1,nbasis 3124 miu1(:)=miu1(:)+CObasa(ibas,imo)*CObasa(jbas,imo)*Dbas(:,ibas,jbas) 3125 miu2(:)=miu2(:)+CObasa(ibas,jmo)*CObasa(jbas,jmo)*Dbas(:,ibas,jbas) 3126 miu12(:)=miu12(:)+( CObasa(ibas,imo)*CObasa(jbas,jmo)+CObasa(ibas,jmo)*CObasa(jbas,imo) )*Dbas(:,ibas,jbas) 3127 end do 3128 end do 3129 miu12=miu12/2 3130 valmiu1=dsqrt(sum(miu1**2)) 3131 valmiu2=dsqrt(sum(miu2**2)) 3132 valmiu12=dsqrt(sum(miu12**2)) 3133 write(*,"(' miu1 (X/Y/Z): ',3f11.6,' Norm:',f11.6,' a.u.')") miu1,valmiu1 3134 write(*,"(' miu2 (X/Y/Z): ',3f11.6,' Norm:',f11.6,' a.u.')") miu2,valmiu2 3135 write(*,"(' miu1-miu2 (X/Y/Z):',3f11.6,' Norm:',f11.6,' a.u.')") miu1-miu2,dsqrt(sum((miu1-miu2)**2)) 3136 write(*,"(' miu12 (X/Y/Z): ',3f11.6,' Norm:',f11.6,' a.u.')") miu12,valmiu12 3137 H_DB=enediff*valmiu12/dsqrt(sum((miu1-miu2)**2)+4*valmiu12**2) 3138 end if 3139 write(*,"(' H_DB:',f14.8,' eV')") H_DB*au2eV 3140 write(*,*) 3141 rmsH_DB=rmsH_DB+H_DB**2 3142 avgE_DB=avgE_DB+enediff 3143 end do 3144 end do 3145 rmsH_DB=dsqrt(rmsH_DB/ndonorMO/nbridgeMO) 3146 avgE_DB=avgE_DB/(ndonorMO*nbridgeMO) 3147 write(*,"(' Average MO energy difference E_DB:',f14.8,'eV')") avgE_DB*au2eV 3148 write(*,"(' Root mean square H_DB:',f14.8,'eV')") rmsH_DB*au2eV 3149 3150 write(*,*) 3151 write(*,*) 3152 write(*,*) " ======== Results between bridge and acceptor ========" 3153 write(*,*) 3154 rmsH_BA=0 3155 avgE_BA=0 3156 do idx=1,nbridgeMO 3157 imo=bridgeMO(idx) 3158 E1=MOene(imo) 3159 q1B=MOcomp(imo,2) 3160 q1A=MOcomp(imo,3) 3161 dq1=q1B-q1A 3162 do jdx=1,nacceptorMO 3163 jmo=acceptorMO(jdx) 3164 E2=MOene(jmo) 3165 enediff=abs(E2-E1) !Always take positive energy difference 3166 q2B=MOcomp(jmo,2) 3167 q2A=MOcomp(jmo,3) 3168 dq2=q2B-q2A 3169 q12_1=0 3170 write(*,"(' **** MO',i4,' (bridge) --- MO',i4,' (acceptor) ****')") imo,jmo 3171 write(*,"(' MO energy difference:',f12.4,' eV')") enediff*au2eV 3172 if (imethod==1) then 3173 do iatmidx=1,nbridgeatm 3174 iatm=bridge(iatmidx) 3175 do ibas=basstart(iatm),basend(iatm) 3176 do jbas=1,nbasis 3177 q12_1=q12_1+( CObasa(ibas,imo)*CObasa(jbas,jmo)+CObasa(ibas,jmo)*CObasa(jbas,imo) )*Sbas(ibas,jbas) 3178 end do 3179 end do 3180 end do 3181 q12_1=q12_1/2 3182 q12_2=0 3183 do iatmidx=1,nacceptoratm 3184 iatm=acceptor(iatmidx) 3185 do ibas=basstart(iatm),basend(iatm) 3186 do jbas=1,nbasis 3187 q12_2=q12_2+( CObasa(ibas,imo)*CObasa(jbas,jmo)+CObasa(ibas,jmo)*CObasa(jbas,imo) )*Sbas(ibas,jbas) 3188 end do 3189 end do 3190 end do 3191 q12_2=q12_2/2 3192 dq12=q12_1-q12_2 3193 H_BA=enediff*abs(dq12)/dsqrt((dq1-dq2)**2+4*dq12**2) 3194 write(*,"(' q1(bridge):',f10.6,' q1(acceptor):',f10.6,' delta_q1:',f10.6)") q1D,q1A,dq1 3195 write(*,"(' q2(bridge):',f10.6,' q2(acceptor):',f10.6,' delta_q2:',f10.6)") q2D,q2A,dq2 3196 write(*,"(' delta_q1 - delta_q2:',f10.6)") dq1-dq2 3197 write(*,"(' q12(bridge):',f10.6,' q12(acceptor):',f10.6,' delta_q12:',f10.6)") q12_1,q12_2,dq12 3198 else if (imethod==2) then 3199 miu1=0 3200 miu2=0 3201 miu12=0 3202 do ibas=1,nbasis !Note that Dbas has already taken the negative sign of electron charge into account 3203 do jbas=1,nbasis 3204 miu1(:)=miu1(:)+CObasa(ibas,imo)*CObasa(jbas,imo)*Dbas(:,ibas,jbas) 3205 miu2(:)=miu2(:)+CObasa(ibas,jmo)*CObasa(jbas,jmo)*Dbas(:,ibas,jbas) 3206 miu12(:)=miu12(:)+( CObasa(ibas,imo)*CObasa(jbas,jmo)+CObasa(ibas,jmo)*CObasa(jbas,imo) )*Dbas(:,ibas,jbas) 3207 end do 3208 end do 3209 miu12=miu12/2 3210 valmiu1=dsqrt(sum(miu1**2)) 3211 valmiu2=dsqrt(sum(miu2**2)) 3212 valmiu12=dsqrt(sum(miu12**2)) 3213 write(*,"(' miu1 (X/Y/Z): ',3f11.6,' Norm:',f11.6,' a.u.')") miu1,valmiu1 3214 write(*,"(' miu2 (X/Y/Z): ',3f11.6,' Norm:',f11.6,' a.u.')") miu2,valmiu2 3215 write(*,"(' miu1-miu2 (X/Y/Z):',3f11.6,' Norm:',f11.6,' a.u.')") miu1-miu2,dsqrt(sum((miu1-miu2)**2)) 3216 write(*,"(' miu12 (X/Y/Z): ',3f11.6,' Norm:',f11.6,' a.u.')") miu12,valmiu12 3217 H_BA=enediff*valmiu12/dsqrt(sum((miu1-miu2)**2)+4*valmiu12**2) 3218 end if 3219 write(*,"(' H_BA:',f14.8,' eV')") H_BA*au2eV 3220 write(*,*) 3221 rmsH_BA=rmsH_BA+H_BA**2 3222 avgE_BA=avgE_BA+enediff 3223 end do 3224 end do 3225 rmsH_BA=dsqrt(rmsH_BA/nbridgeMO/nacceptorMO) 3226 avgE_BA=avgE_BA/(nbridgeMO*nacceptorMO) 3227 write(*,"(' Average MO energy difference E_BA:',f14.8,'eV')") avgE_BA*au2eV 3228 write(*,"(' Root mean square H_BA:',f14.8,'eV')") rmsH_BA*au2eV 3229 write(*,*) 3230! write(*,*) " ======== Final result ========" 3231! write(*,*) 3232! H_SE=rmsH_DB*rmsH_BA/(avgE_DA-avgE_DB/2) 3233! write(*,"(' H_SE (Superexchange): ',f14.8,'eV')") H_SE 3234! H_DBA=rmsH_DA+H_SE 3235! write(*,"(' H_DBA (H_DA via bridge):',f14.8,'eV')") H_DBA 3236end if 3237 3238deallocate(donorMO,bridgeMO,acceptorMO) 3239 3240goto 130 3241end subroutine 3242 3243 3244 3245!---------- Perform Pipek-Mezey orbital localization 3246!The algorithm can be found in J. Comput. Chem., 14, 736 (1993) 3247!Performing E2 analysis based on LMO oribtal is not appropriate, the non-diagonal elements of Fock matrix are almost zero, & 3248!this phenonmenon can also be seen in NLMO Fock matrix ($NBO FNLMO $END). This is because we don't allow mixure between occupied and unoccupied MOs 3249! 3250!The final wavefunction can be exported as .fch, I don't select .molden because it doesn't record atomic charge, this will be problematic when ECP is used 3251subroutine pipek_mezey 3252use defvar 3253use util 3254implicit real*8 (a-h,o-z) 3255integer :: maxcyc=80,ireload=1,idoene=0,domark(4) 3256real*8 :: arrayi(nbasis),arrayj(nbasis),crit=1D-4,tmpbasarr(nbasis),tmpprimarr(nprims) 3257real*8,pointer :: Cmat(:,:) 3258real*8,allocatable :: FmatA(:,:),FmatB(:,:),FLMOA(:,:),FLMOB(:,:) 3259character c200tmp*200,typestr*4 3260if (wfntype==2.or.wfntype==3.or.wfntype==4) then 3261 write(*,*) "Error: This function only works for restricted or unrestricted SCF wavefunction!" 3262 write(*,*) "Press ENTER to return" 3263 read(*,*) 3264 return 3265end if 3266if (.not.allocated(CObasa)) then 3267 write(*,*) "Error: Basis function information was not provided by your input file!" 3268 write(*,*) "Press ENTER to return" 3269 read(*,*) 3270 return 3271end if 3272 3273do while(.true.) 3274 if (idoene==1) write(*,*) "-4 If calculate and print orbital energies, current: Yes" 3275 if (idoene==0) write(*,*) "-4 If calculate and print orbital energies, current: No" 3276 if (ireload==1) write(*,*) "-3 If finally reload newly generated .fch file, current: Yes" 3277 if (ireload==0) write(*,*) "-3 If finally reload newly generated .fch file, current: No" 3278 write(*,"(a,f12.8)") " -2 Set criterion of convergence, current:",crit 3279 write(*,"(a,i4)") " -1 Set maximum cycles, current:",maxcyc 3280 write(*,*) "0 Return" 3281 write(*,*) "1 Localizing occupied orbitals only" 3282 write(*,*) "2 Localizing both occupied and unoccupied orbitals separately" 3283 read(*,*) isel 3284 if (isel==0) then 3285 return 3286 else if (isel==-1) then 3287 write(*,*) "Input maximum cycles, e.g. 30" 3288 read(*,*) maxcyc 3289 else if (isel==-2) then 3290 write(*,*) "Input criterion of convergence, e.g. 0.0001" 3291 read(*,*) crit 3292 else if (isel==-3) then 3293 if (ireload==1) then 3294 ireload=0 3295 else if (ireload==0) then 3296 ireload=1 3297 end if 3298 else if (isel==-4) then 3299 if (idoene==1) then 3300 idoene=0 3301 else if (idoene==0) then 3302 write(*,"(a)") " Input the file recording Fock matrix in original basis functions in lower triangular form, e.g. C:\fock.txt" 3303 read(*,"(a)") c200tmp 3304 inquire(file=c200tmp,exist=alive) 3305 if (alive.eqv. .false.) then 3306 write(*,*) "Error: Unable to find this file!" 3307 cycle 3308 end if 3309 open(10,file=c200tmp,status="old") 3310 if (.not.allocated(FmatA)) allocate(FmatA(nbasis,nbasis)) 3311 read(10,*) ((FmatA(i,j),j=1,i),i=1,nbasis) 3312 do i=1,nbasis !Fill upper triangular part 3313 do j=i+1,nbasis 3314 FmatA(i,j)=FmatA(j,i) 3315 end do 3316 end do 3317 if (wfntype==1) then 3318 if (.not.allocated(FmatB)) allocate(FmatB(nbasis,nbasis)) 3319 read(10,*) ((FmatB(i,j),j=1,i),i=1,nbasis) 3320 do i=1,nbasis 3321 do j=i+1,nbasis 3322 FmatB(i,j)=FmatB(j,i) 3323 end do 3324 end do 3325 end if 3326 close(10) 3327 write(*,*) "Fock matrix loaded successfully!" 3328 write(*,*) 3329 idoene=1 3330 end if 3331 else if (isel==1.or.isel==2) then 3332 exit 3333 end if 3334end do 3335 3336call walltime(iwalltime1) 3337CALL CPU_TIME(time_begin) 3338 3339domark=0 3340if (wfntype==0) then 3341 domark(1)=1 3342 if (isel==2) domark(2)=1 3343else if (wfntype==1) then 3344 domark(1)=1 3345 domark(3)=1 3346 if (isel==2) domark=1 3347end if 3348 3349!Alpha-occ,Alpha-vir,Beta-occ,Beta-vir 3350do itime=1,4 3351 if (domark(itime)==0) cycle 3352 if (itime<=2) then 3353 Cmat=>CObasa 3354 else 3355 Cmat=>CObasb 3356 end if 3357 if (itime==1) then 3358 nmobeg=1 3359 nmoend=naelec 3360 else if (itime==2) then 3361 nmobeg=naelec+1 3362 nmoend=nbasis 3363 else if (itime==3) then 3364 nmobeg=1 3365 nmoend=nbelec 3366 else if (itime==4) then 3367 nmobeg=nbelec+1 3368 nmoend=nbasis 3369 end if 3370 3371 if (wfntype==0) then 3372 if (itime==1) write(*,"(/,a)") " Localizing occupied orbitals..." 3373 if (itime==2) write(*,"(/,a)") " Localizing unoccupied orbitals..." 3374 else if (wfntype==1) then 3375 if (itime==1) write(*,"(/,a)") " Localizing alpha occupied orbitals..." 3376 if (itime==2) write(*,"(/,a)") " Localizing alpha unoccupied orbitals..." 3377 if (itime==3) write(*,"(/,a)") " Localizing beta occupied orbitals..." 3378 if (itime==4) write(*,"(/,a)") " Localizing beta unoccupied orbitals..." 3379 end if 3380 3381 Pvalold=0 3382 do icyc=1,maxcyc 3383 do imo=nmobeg,nmoend-1 3384 do jmo=imo+1,nmoend 3385 Aval=0 3386 Bval=0 3387nthreads=getNThreads() 3388!$OMP parallel shared(Aval,Bval) private(Avalpriv,Bvalpriv,Qij,Qii,Qjj) num_threads(nthreads) 3389 Avalpriv=0 3390 Bvalpriv=0 3391!$OMP do schedule(DYNAMIC) 3392 do iatm=1,ncenter 3393 Qij=Qval(Cmat,imo,jmo,iatm) 3394 Qii=Qval(Cmat,imo,imo,iatm) 3395 Qjj=Qval(Cmat,jmo,jmo,iatm) 3396 Avalpriv=Avalpriv+( Qij**2-(Qii-Qjj)**2/4D0 ) 3397 Bvalpriv=Bvalpriv+( Qij*(Qii-Qjj) ) 3398 end do 3399!$OMP END DO 3400!$OMP CRITICAL 3401 Aval=Aval+Avalpriv 3402 Bval=Bval+Bvalpriv 3403!$OMP END CRITICAL 3404!$OMP END PARALLEL 3405 if (Aval==0.and.Bval==0) cycle 3406 gamma=sign(1D0,Bval)*acos(-Aval/dsqrt(Aval**2+Bval**2))/4D0 3407 arrayi=cos(gamma)*Cmat(:,imo)+sin(gamma)*Cmat(:,jmo) 3408 arrayj=-sin(gamma)*Cmat(:,imo)+cos(gamma)*Cmat(:,jmo) 3409 Cmat(:,imo)=arrayi 3410 Cmat(:,jmo)=arrayj 3411 ! write(*,*) imo,jmo,gamma,Aval,Bval 3412 ! read(*,*) 3413 end do 3414 end do 3415 Pval=0 3416 do imo=nmobeg,nmoend 3417 do iatm=1,ncenter 3418 Pval=Pval+Qval(Cmat,imo,imo,iatm)**2 3419 end do 3420 end do 3421 deltaPval=Pval-Pvalold 3422 write(*,"(' Cycle:',i5,' P:',f16.8,' Delta P:',f16.8)") icyc,Pval,deltaPval 3423 if (abs(deltaPval)<crit) exit 3424 Pvalold=Pval 3425 end do 3426 if (icyc==maxcyc+1) then 3427 write(*,*) "Warning: Convergence failed!" 3428 else 3429 write(*,"(a)") " Successfully converged!" 3430 end if 3431end do 3432 3433CALL CPU_TIME(time_end) 3434call walltime(iwalltime2) 3435write(*,"(/,' Calculation took up CPU time',f12.2,'s, wall clock time',i10,'s')") time_end-time_begin,iwalltime2-iwalltime1 3436 3437!Print orbital energies, sort orbital according to energies and do E2 analysis 3438if (idoene==1) then 3439 !Do Alpha part or closed-shell orbitals 3440 allocate(FLMOA(nbasis,nbasis)) 3441 FLMOA=matmul(matmul(transpose(CObasa),FmatA),CObasa) 3442 if (isel==1) nmoend=naelec 3443 if (isel==2) nmoend=nbasis 3444 do iorb=1,nmoend 3445 MOene(iorb)=FLMOA(iorb,iorb) 3446 end do 3447 do iorb=1,nmoend 3448 do jorb=iorb+1,nmoend 3449 if (MOene(iorb)>MOene(jorb)) then 3450 tmpbasarr=CObasa(:,iorb) 3451 CObasa(:,iorb)=CObasa(:,jorb) 3452 CObasa(:,jorb)=tmpbasarr 3453 tmpprimarr=CO(iorb,:) 3454 CO(iorb,:)=CO(jorb,:) 3455 CO(jorb,:)=tmpprimarr 3456 tmpene=MOene(iorb) 3457 MOene(iorb)=MOene(jorb) 3458 MOene(jorb)=tmpene 3459 end if 3460 end do 3461 end do 3462 !Do beta part 3463 if (wfntype==1) then 3464 allocate(FLMOB(nbasis,nbasis)) 3465 FLMOB=matmul(matmul(transpose(CObasb),FmatB),CObasb) 3466 if (isel==1) nmoend=nbelec 3467 if (isel==2) nmoend=nbasis 3468 do iorb=1,nmoend 3469 MOene(nbasis+iorb)=FLMOB(iorb,iorb) 3470 end do 3471 do iorb=1,nmoend 3472 do jorb=iorb+1,nmoend 3473 if (MOene(nbasis+iorb)>MOene(nbasis+jorb)) then 3474 tmpbasarr=CObasb(:,iorb) 3475 CObasb(:,iorb)=CObasb(:,jorb) 3476 CObasb(:,jorb)=tmpbasarr 3477 tmpprimarr=CO(nbasis+iorb,:) 3478 CO(nbasis+iorb,:)=CO(nbasis+jorb,:) 3479 CO(nbasis+jorb,:)=tmpprimarr 3480 tmpene=MOene(nbasis+iorb) 3481 MOene(nbasis+iorb)=MOene(nbasis+jorb) 3482 MOene(nbasis+jorb)=tmpene 3483 end if 3484 end do 3485 end do 3486 end if 3487 !Print energies, those not localized are not printed 3488 write(*,*) "Energies of localized orbitals:" 3489 do iorb=1,nmo 3490 if (isel==1) then 3491 if (wfntype==0) then 3492 if (iorb>nint(naelec)) cycle 3493 else if (wfntype==1) then 3494 if (MOtype(iorb)==1.and.iorb>nint(naelec)) cycle 3495 if (MOtype(iorb)==2.and.iorb>nbasis+nint(nbelec)) cycle 3496 end if 3497 end if 3498 if (MOtype(iorb)==0) typestr="A+B" 3499 if (MOtype(iorb)==1) typestr="A" 3500 if (MOtype(iorb)==2) typestr="B" 3501 write(*,"(i6,' Energy (a.u.):',f16.8,' Type: ',a,' Occ:',f4.1)") iorb,MOene(iorb),typestr,MOocc(iorb) 3502 end do 3503 if (isel==1) write(*,*) "Energies of unoccupied orbitals are not updated since they were not localized" 3504 3505 !Second-order perturbation analysis between occupied and virtual orbitals (like NBO E2), this only works when both of them have been localized 3506 !This part is commented since it don't print any useful result, because it is easy to proved that Fock element between occupied and virtual LMOs are exactly zero 3507! if (isel==1) then 3508! write(*,*) " Note: E(2) analysis is skipped since virtual orbitals were not localized" 3509! else if (isel==2) then 3510! write(*,*) "Second-order perturbation theory analysis of interaction energy:" 3511! !Regenerated Fock matrix in LMO, since they have been sorted 3512! FLMOA=matmul(matmul(transpose(CObasa),FmatA),CObasa) 3513! call showmatgau(FLMOA) 3514! read(*,*) 3515! coeff=2 3516! do iocc=1,naelec 3517! do ivir=naelec+1,nbasis 3518! E2val=-coeff* FLMOA(iocc,ivir)**2/(MOene(ivir)-MOene(iocc)) 3519! qCT=coeff* ( FLMOA(iocc,ivir)/(MOene(ivir)-MOene(iocc)) )**2 3520! if (abs(E2val*au2kcal)>0.2D0) then 3521! write(*,"(' Donor:',i5,' - Acceptor:',i5,' E(2):',f7.2,' kcal/mol q_CT:',f10.5)") iocc,ivir,E2val*au2kcal,qCT 3522! write(*,"(3f16.10)") FLMOA(iocc,ivir),MOene(ivir)-MOene(iocc) 3523! end if 3524! end do 3525! end do 3526! end if 3527end if 3528 3529write(*,*) 3530call outfch("new.fch",10) 3531write(*,*) "The localized orbitals have been exported as new.fch in current folder" 3532if (ireload==1) then 3533 call dealloall 3534 write(*,*) "Loading new.fch..." 3535 call readfch("new.fch",1) 3536 write(*,"(a)") " Loading finished, now you can use main function 0 to visualize them as isosurface" 3537end if 3538end subroutine 3539 3540!------ Calculate Q value used in Pipek-Mezey localization 3541real*8 function Qval(Cmat,imo,jmo,iatm) 3542use defvar 3543implicit real*8 (a-h,o-z) 3544real*8 Cmat(nbasis,nbasis) 3545Qval=0 3546do ibas=basstart(iatm),basend(iatm) 3547 Qval=Qval+sum( (Cmat(:,imo)*Cmat(ibas,jmo)+Cmat(ibas,imo)*Cmat(:,jmo)) *Sbas(ibas,:) ) 3548end do 3549Qval=Qval/2D0 3550end function 3551 3552 3553 3554 3555 3556!!------------ Integrate any real space function within isosurface of a real space function 3557!I use the same data structure as basin analysis to illustrate definition of isosurfaces 3558subroutine intisosurface 3559use defvar 3560use util 3561use function 3562use basinintmod !Use its vec26 array 3563implicit real*8 (a-h,o-z) 3564integer :: ifunciso=13,ifuncint=1 3565integer,allocatable :: mergelist(:),grididx(:,:,:),dogrid(:,:) 3566character :: defdomain*20="<0.5",c1000tmp*1000 3567 3568if (allocated(gridxyz)) deallocate(gridxyz) 3569if (allocated(domainsize)) deallocate(domainsize) 3570if (allocated(domaingrid)) deallocate(domaingrid) 3571do while(.true.) 3572 write(*,*) 3573 if (allocated(cubmat)) write(*,*) "-1 Directly use the grid data in memory to yield domains" 3574 write(*,*) "0 Return" 3575 write(*,*) "1 Calculate grid data and yield domains" 3576 write(*,"(a,i5)") " 2 Select real space function to be calculated, current:",ifunciso 3577 write(*,"(a,a)") " 3 Set criterion for defining domain, current: ",trim(defdomain) 3578 read(*,*) isel 3579 if (isel==0) then 3580 return 3581 else if (isel==1.or.isel==-1) then 3582 exit 3583 else if (isel==2) then 3584 call selfunc_interface(ifunciso) 3585 else if (isel==3) then 3586 write(*,*) "Input the definition, e.g. <0.05" 3587 write(*,*) "Note: The first character must be < or >" 3588 read(*,"(a)") defdomain 3589 end if 3590end do 3591 3592!Set grid and generate grid data 3593if (isel==1) then 3594 call setgridfixspc 3595 dvol=dx*dy*dz 3596 if (allocated(cubmat)) deallocate(cubmat) 3597 allocate(cubmat(nx,ny,nz)) 3598 call savecubmat(ifunciso,0,iorbsel) 3599end if 3600 3601!Count the number of grids satisfying the criterion 3602read(defdomain(2:),*) valiso 3603if (defdomain(1:1)=='<') then 3604 ngrid=count(cubmat<valiso) 3605else if (defdomain(1:1)=='>') then 3606 ngrid=count(cubmat>valiso) 3607else 3608 write(*,*) "Error: Parsing of the domain definition failed!" 3609 return 3610end if 3611write(*,"(/,' The number of grids satisfied the criterion:',i10)") ngrid 3612 3613!Clustering grids that satisfied criterion to domain 3614!The idea is very clever: I assign each grid with a different index, and frequently perform iteration, in each cycle all grids (that meets isovalue criterion) && 3615!are run over, and index of each grid is set to that of one of the nearest 26 grids as long as its index is larger than current grid. Finally, & 3616!grids in each domian will have identical indices 3617write(*,*) "Clustering domains..." 3618call walltime(iwalltime1) 3619CALL CPU_TIME(time_begin) 3620!dogrid: XYZ index of useful grids 3621!gridxyz: XYZ coordinate of useful grids 3622!grididx: currently records initial index of all gris 3623allocate(dogrid(ngrid,3),grididx(nx,ny,nz),gridxyz(ngrid,3)) 3624grididx=-1 !Irrelevant grids have very small value 3625idx=0 3626do iz=1,nz 3627 do iy=1,ny 3628 do ix=1,nx 3629 if ((defdomain(1:1)=='<'.and.cubmat(ix,iy,iz)<valiso).or.(defdomain(1:1)=='>'.and.cubmat(ix,iy,iz)>valiso)) then 3630 idx=idx+1 3631 dogrid(idx,1)=ix 3632 dogrid(idx,2)=iy 3633 dogrid(idx,3)=iz 3634 grididx(ix,iy,iz)=idx 3635 gridxyz(idx,1)=orgx+(ix-1)*dx 3636 gridxyz(idx,2)=orgy+(iy-1)*dy 3637 gridxyz(idx,3)=orgz+(iz-1)*dz 3638 end if 3639 end do 3640 end do 3641end do 3642 3643call setupmovevec 3644icyc=0 3645do while(.true.) 3646 iupdate=0 3647 icyc=icyc+1 3648! write(*,*) icyc 3649 do itmp=1,ngrid 3650 ix=dogrid(itmp,1) 3651 iy=dogrid(itmp,2) 3652 iz=dogrid(itmp,3) 3653 do imove=1,26 3654 ixtmp=ix+vec26x(imove) 3655 iytmp=iy+vec26y(imove) 3656 iztmp=iz+vec26z(imove) 3657 if (ixtmp<1.or.ixtmp>nx.or.iytmp<1.or.iytmp>ny.or.iztmp<1.or.iztmp>nz) then 3658 write(*,"(a)") " Error: Some domains may be truncated! You should enlarge extension size or completely re-defined grid to avoid this circumstance!" 3659 write(*,*) "Press Enter to exit" 3660 read(*,*) 3661 return 3662 end if 3663 if (grididx(ix,iy,iz)<grididx(ixtmp,iytmp,iztmp)) then 3664 grididx(ix,iy,iz)=grididx(ixtmp,iytmp,iztmp) 3665 iupdate=1 3666 exit 3667 end if 3668 end do 3669 end do 3670 if (iupdate==0) exit 3671end do 3672! do itmp=1,ngrid 3673! ix=dogrid(itmp,1) 3674! iy=dogrid(itmp,2) 3675! iz=dogrid(itmp,3) 3676! write(*,*) itmp,grididx(ix,iy,iz) 3677! end do 3678! read(*,*) 3679ndone=0 3680ndomain=0 3681do while(.true.) 3682 ndomain=ndomain+1 3683 mintmp=1000000 3684 do itmp=1,ngrid 3685 idxtmp=grididx(dogrid(itmp,1),dogrid(itmp,2),dogrid(itmp,3)) 3686 if (idxtmp<ndomain) cycle 3687 if (idxtmp<mintmp) mintmp=idxtmp 3688 end do 3689 do itmp=1,ngrid 3690 idxtmp=grididx(dogrid(itmp,1),dogrid(itmp,2),dogrid(itmp,3)) 3691 if (idxtmp==mintmp) then 3692 grididx(dogrid(itmp,1),dogrid(itmp,2),dogrid(itmp,3))=ndomain 3693 ndone=ndone+1 3694 end if 3695 end do 3696 if (ndone==ngrid) exit 3697! write(*,*) ndomain,ndone 3698end do 3699!Generate domainsize and domaingrid (grid index that contained in each domain) 3700allocate(domainsize(ndomain),domaingrid(ndomain,ngrid)) 3701do idom=1,ndomain 3702 j=0 3703 do itmp=1,ngrid 3704 if (grididx(dogrid(itmp,1),dogrid(itmp,2),dogrid(itmp,3))==idom) then 3705 j=j+1 3706 domaingrid(idom,j)=itmp 3707 end if 3708 end do 3709 domainsize(idom)=j 3710end do 3711 3712write(*,*) "Clustering domains finished!" 3713CALL CPU_TIME(time_end) 3714call walltime(iwalltime2) 3715write(*,"(' Clustering took up CPU time',f12.2,'s, wall clock time',i10,'s')") time_end-time_begin,iwalltime2-iwalltime1 3716 3717do idom=1,ndomain 3718 write(*,"(' Domain:',i6,' Grids:',i8)") idom,domainsize(idom) 3719end do 3720 3721do while(.true.) 3722 write(*,*) 3723 write(*,*) "-1 Merge specific domains" 3724 write(*,*) "0 Exit" 3725 write(*,*) "1 Perform integration for a domain" 3726 write(*,*) "2 Perform integration for all domains" 3727 write(*,*) "3 Visualize domains" 3728 write(*,"(a,i5)") " 4 Select the real space function to be integrated, current:",ifuncint 3729 write(*,*) "5 Calculate q_bind index for a domain" 3730 read(*,*) isel2 3731 if (isel2==0) then 3732 return 3733 else if (isel2==-1) then 3734 if (ndomain<2) then 3735 write(*,*) "Error: At least two domains must be presented!" 3736 cycle 3737 end if 3738 write(*,*) "Input indices of the domains you want to merge, e.g. 4,5,8-10" 3739 read(*,"(a)") c1000tmp 3740 call str2arr(c1000tmp,nmerge) !Find how many terms 3741 allocate(mergelist(nmerge)) 3742 call str2arr(c1000tmp,nmerge,mergelist) 3743 call sorti4(mergelist) 3744 idom=mergelist(1) 3745 do jdx=nmerge,2,-1 3746 jdom=mergelist(jdx) 3747 nsizei=domainsize(idom) 3748 nsizej=domainsize(jdom) 3749 domainsize(idom)=nsizei+nsizej 3750 domaingrid(idom,nsizei+1:nsizei+nsizej)=domaingrid(jdom,1:nsizej) 3751 ndomain=ndomain-1 3752 domainsize(jdom:ndomain)=domainsize(jdom+1:ndomain+1) 3753 domaingrid(jdom:ndomain,:)=domaingrid(jdom+1:ndomain+1,:) 3754 end do 3755 deallocate(mergelist) 3756 write(*,"(a,i6)") " Done! The domains you selected have been merged as domain",idom 3757 write(*,*) "Size of current domains:" 3758 do idom=1,ndomain 3759 write(*,"(' Domain:',i6,' Grids:',i8)") idom,domainsize(idom) 3760 end do 3761 else if (isel2==1) then 3762 write(*,*) "Input the index of the domain to be integrated, e.g. 3" 3763 read(*,*) intdom 3764 if (intdom<1.or.intdom>ndomain) then 3765 write(*,"(a)") " Error: The index of the domain to be integrated is incorrect" 3766 cycle 3767 end if 3768 valint=0 3769 volint=0 3770 valmin=1D100 3771 valmax=-1D100 3772 do igrd=1,domainsize(intdom) 3773 idx=domaingrid(intdom,igrd) 3774 tmpval=calcfuncall(ifuncint,gridxyz(idx,1),gridxyz(idx,2),gridxyz(idx,3)) 3775 valint=valint+tmpval 3776 volint=volint+1 3777 if (tmpval<valmin) valmin=tmpval 3778 if (tmpval>valmax) valmax=tmpval 3779 end do 3780 avgval=valint/domainsize(intdom) 3781 valint=valint*dvol 3782 volint=volint*dvol 3783 write(*,"(' Integration result:',E20.10,' a.u.')") valint 3784 write(*,"(' Volume:',f12.6,' Bohr^3 (',f12.6,' Angstrom^3 )')") volint,volint*b2a**3 3785 write(*,"(' Average:',E20.10)") avgval 3786 write(*,"(' Maximum:',E20.10,' Minimum:',E20.10)") valmax,valmin 3787 else if (isel2==2) then 3788 write(*,*) "Domain Integral (a.u.) Volume (Bohr^3) Average" 3789 valinttot=0 3790 volinttot=0 3791 do intdom=1,ndomain 3792 valint=0 3793 volint=0 3794 do igrd=1,domainsize(intdom) 3795 idx=domaingrid(intdom,igrd) 3796 valint=valint+calcfuncall(ifuncint,gridxyz(idx,1),gridxyz(idx,2),gridxyz(idx,3)) 3797 volint=volint+1 3798 end do 3799 avgval=valint/domainsize(intdom) 3800 valint=valint*dvol 3801 volint=volint*dvol 3802 write(*,"(i6,E20.10,f17.6,E20.10)") intdom,valint,volint,avgval 3803 valinttot=valinttot+valint 3804 volinttot=volinttot+volint 3805 end do 3806 write(*,"(' Integration result of all domains:',E20.10,' a.u.')") valinttot 3807 write(*,"(' Volume of all domains:',f13.6,' Bohr^3 ',f13.6,' Angstrom^3')") volinttot,volinttot*b2a**3 3808 else if (isel2==3) then 3809 idrawdomain=1 3810 aug3Dold=aug3D 3811 if (aug3D<3) aug3D=3 !Often we set extension distance to zero, e.g. RDG, in this case the molecule will be truncated, therefore here temporarily augment it 3812 aug3D=aug3Dold 3813 idrawdomain=0 3814 else if (isel2==4) then 3815 call selfunc_interface(ifuncint) 3816 else if (isel2==5) then 3817 write(*,*) "Input the index of the domain to be integrated, e.g. 3" 3818 read(*,*) intdom 3819 if (intdom<1.or.intdom>ndomain) then 3820 write(*,"(a)") " Error: The index of the domain to be integrated is incorrect" 3821 cycle 3822 end if 3823 qatt=0 3824 qrep=0 3825 expfac=4D0/3D0 3826 volneg=0 3827 volpos=0 3828 do igrd=1,domainsize(intdom) 3829 idx=domaingrid(intdom,igrd) 3830 call signlambda2rho_RDG(gridxyz(idx,1),gridxyz(idx,2),gridxyz(idx,3),sl2r,RDG,rho) 3831 if (sl2r<0) then 3832 qatt=qatt+rho**expfac 3833 volneg=volneg+1 3834 else 3835 qrep=qrep+rho**expfac 3836 volpos=volpos+1 3837 end if 3838 end do 3839 qatt=qatt*dvol 3840 qrep=qrep*dvol 3841 qbind=-(qatt-qrep) 3842 volneg=volneg*dvol 3843 volpos=volpos*dvol 3844 write(*,"(' q_att: 'f16.8,' a.u.')") qatt 3845 write(*,"(' q_rep: 'f16.8,' a.u.')") qrep 3846 write(*,"(' q_bind:'f16.8,' a.u.')") qbind 3847 write(*,"(' Volume (lambda2<0):',f13.6,' Bohr^3 ',f13.6,' Angstrom^3')") volneg 3848 write(*,"(' Volume (lambda2>0):',f13.6,' Bohr^3 ',f13.6,' Angstrom^3')") volpos 3849 write(*,"(' Volume (Total): ',f13.6,' Bohr^3 ',f13.6,' Angstrom^3')") volneg+volpos 3850 end if 3851end do 3852end subroutine 3853 3854 3855!------ Calculate electron correlation index proposed by Matito et al. 3856subroutine elecorridx 3857use defvar 3858real*8 occ(nmo),I_ND,I_D,I_T 3859write(*,*) "Citation: Phys. Chem. Chem. Phys., 18, 24015 (2016)" 3860I_ND=0 3861I_D=0 3862if (wfntype==3) then 3863 occ=MOocc/2 3864 where(occ>1) occ=1 !Remove unphysical occupation number larger than unity 3865 where(occ<0) occ=0 !Remove unphysical negative occupation number 3866 do i=1,nmo 3867 I_D=I_D+ dsqrt(occ(i)*(1-occ(i))) - 2*occ(i)*(1-occ(i)) 3868 end do 3869 I_D=I_D/4 3870 do i=1,nmo 3871 I_ND=I_ND+ occ(i)*(1-occ(i)) 3872 end do 3873 I_ND=I_ND/2 3874 !Above we only consider half part, another part is identical to that, so double the result 3875 I_D=I_D*2 3876 I_ND=I_ND*2 3877else if (wfntype==4) then 3878 occ=MOocc 3879 where(occ>1) occ=1 3880 where(occ<0) occ=0 3881 do i=1,nmo 3882 I_D=I_D+ dsqrt(occ(i)*(1-occ(i))) - 2*occ(i)*(1-occ(i)) 3883 end do 3884 I_D=I_D/4 3885 do i=1,nmo 3886 I_ND=I_ND+ occ(i)*(1-occ(i)) 3887 end do 3888 I_ND=I_ND/2 3889end if 3890I_T=I_ND+I_D 3891write(*,"(' Nondynamic correlation index:',f12.8)") I_ND 3892write(*,"(' Dynamic correlation index: ',f12.8)") I_D 3893write(*,"(' Total correlation index: ',f12.8)") I_T 3894end subroutine 3895 3896 3897 3898!!------ Generate natural orbitals based on the density matrix in .fch/.fchk file 3899subroutine fch_gennatorb 3900use util 3901use defvar 3902implicit real*8 (a-h,o-z) 3903real*8,allocatable :: tmparr(:),Pspin(:,:) 3904real*8 Xmat(nbasis,nbasis) 3905character selectyn,denstype*10,locstr*40 3906if (ifiletype/=1) then 3907 write(*,*) "Error: .fch/.fchk should be used as input file for this function" 3908 write(*,*) "Press ENTER to return" 3909 read(*,*) 3910 return 3911end if 3912write(*,*) "Input type of density, e.g. SCF, MP2, CI, CC, MP4..." 3913write(*,*) "e.g. If the .fch was produced under MP2, you may input ""SCF"" or ""MP2""" 3914read(*,"(a)") denstype 3915write(locstr,"('Total ',a,' Density')") trim(denstype) 3916open(10,file=filename,status="old") 3917call loclabel(10,trim(locstr),ifoundDM) 3918if (ifoundDM==0) then 3919 write(*,"(' Error: Unable to find ""',a,'"" from the input file')") trim(locstr) 3920 write(*,*) "Press ENTER to return" 3921 read(*,*) 3922 return 3923end if 3924iNOtype=1 3925if (wfntype==1.or.wfntype==2.or.wfntype==4) then 3926 write(*,*) "Select natural orbitals you want to obtain" 3927 write(*,*) "1 Spatial natural orbitals (diagonalization of total density matrix)" 3928 write(*,*) "2 Alpha and beta natural orbitals (diagonalization of respective DM)" 3929 write(*,*) "3 Spin natural orbitals (diagonalization of spin density matrix)" 3930 read(*,*) iNOtype 3931end if 3932 3933write(*,*) "Loading density matrix..." 3934!Load total density matrix 3935Ptot=0D0 3936call loclabel(10,trim(locstr)) 3937read(10,*) 3938read(10,"(5(1PE16.8))") ((Ptot(i,j),j=1,i),i=1,nbasis) 3939Ptot=Ptot+transpose(Ptot) 3940do i=1,nbasis 3941 Ptot(i,i)=Ptot(i,i)/2D0 3942end do 3943!Load spin density matrix 3944if (iNOtype>1) then 3945 allocate(Pspin(nbasis,nbasis)) 3946 Pspin=0 3947 read(10,*) 3948 read(10,"(5(1PE16.8))") ((Pspin(i,j),j=1,i),i=1,nbasis) 3949 Pspin=Pspin+transpose(Pspin) 3950 do i=1,nbasis 3951 Pspin(i,i)=Pspin(i,i)/2D0 3952 end do 3953 Palpha=(Ptot+Pspin)/2D0 3954 Pbeta=(Ptot-Pspin)/2D0 3955end if 3956close(10) 3957write(*,*) "Density matrix was loaded from .fch/.fchk file" 3958 3959!To produce natural orbitals, we need to convert P to orthogonalized basis and then diagonalize it 3960allocate(tmparr(nbasis)) 3961write(*,*) 3962if (iNOtype==1.or.iNOtype==3) then 3963 if (iNOtype==1) write(*,*) "Generating NOs, please wait..." 3964 if (iNOtype==3) write(*,*) "Generating SNOs, please wait..." 3965 call symmorthomat(nbasis,Sbas,Xmat,0) !Xmat=S^0.5 3966 if (iNOtype==1) call diagsymat(matmul(matmul(transpose(Xmat),Ptot),Xmat),CObasa,MOocc,ierror) !CObasa now is NOs in orthogonalized basis 3967 if (iNOtype==3) call diagsymat(matmul(matmul(transpose(Xmat),Pspin),Xmat),CObasa,MOocc,ierror) !CObasa now is SNOs in orthogonalized basis 3968 MOene=0 3969 call symmorthomat(nbasis,Sbas,Xmat,1) !Xmat=S^-0.5 3970 CObasa=matmul(Xmat,CObasa) !Convert CObasa to original basis 3971 !Sort NOs according to occupation number 3972 do i=1,nbasis 3973 do j=i+1,nbasis 3974 if (MOocc(i)<MOocc(j)) then 3975 tmpocc=MOocc(i) 3976 MOocc(i)=MOocc(j) 3977 MOocc(j)=tmpocc 3978 tmparr=CObasa(:,i) 3979 CObasa(:,i)=CObasa(:,j) 3980 CObasa(:,j)=tmparr 3981 end if 3982 end do 3983 end do 3984 if (wfntype==1.or.wfntype==4) then !Then wfntype will be 3, deallocate useless arrays and resize arrays 3985 deallocate(CObasb,Palpha,Pbeta,MOene,tmparr) 3986 allocate(MOene(nbasis)) 3987 MOene=0 3988 allocate(tmparr(nmo)) 3989 tmparr=MOocc 3990 deallocate(MOocc) 3991 allocate(MOocc(nbasis)) 3992 MOocc=tmparr(1:nbasis) 3993 nmo=nbasis 3994 end if 3995 write(*,*) "Occupation numbers:" 3996 write(*,"(6f12.6)") MOocc 3997 wfntype=3 3998else 3999 write(*,*) "Generating alpha and beta NOs, please wait..." 4000 call symmorthomat(nbasis,Sbas,Xmat,0)!Xmat=S^0.5 4001 call diagsymat(matmul(matmul(transpose(Xmat),Palpha),Xmat),CObasa,MOocc(1:nbasis),ierror) 4002 call symmorthomat(nbasis,Sbas,Xmat,1) !Xmat=S^-0.5 4003 CObasa=matmul(Xmat,CObasa) !Convert CObasa to original basis, as what we did in HF calculation 4004 MOene(1:nbasis)=0 4005 do i=1,nbasis 4006 do j=i+1,nbasis 4007 if (MOocc(i)<MOocc(j)) then 4008 tmpocc=MOocc(i) 4009 MOocc(i)=MOocc(j) 4010 MOocc(j)=tmpocc 4011 tmparr=CObasa(:,i) 4012 CObasa(:,i)=CObasa(:,j) 4013 CObasa(:,j)=tmparr 4014 end if 4015 end do 4016 end do 4017 write(*,*) "Occupation numbers of Alpha NOs:" 4018 write(*,"(6f12.6)") MOocc(1:nbasis) 4019 write(*,*) 4020 call symmorthomat(nbasis,Sbas,Xmat,0) !Xmat=S^0.5 4021 call diagsymat(matmul(matmul(transpose(Xmat),Pbeta),Xmat),CObasb,MOocc(nbasis+1:nmo),ierror) 4022 call symmorthomat(nbasis,Sbas,Xmat,1) !Xmat=S^-0.5 4023 CObasb=matmul(Xmat,CObasb) 4024 MOene(nbasis+1:nmo)=0 4025 do i=1,nbasis 4026 ii=i+nbasis 4027 do j=i+1,nbasis 4028 jj=j+nbasis 4029 if (MOocc(ii)<MOocc(jj)) then 4030 tmpocc=MOocc(ii) 4031 MOocc(ii)=MOocc(jj) 4032 MOocc(jj)=tmpocc 4033 tmparr=CObasb(:,i) 4034 CObasb(:,i)=CObasb(:,j) 4035 CObasb(:,j)=tmparr 4036 end if 4037 end do 4038 end do 4039 write(*,*) "Occupation numbers of Beta NOs:" 4040 write(*,"(6f12.6)") MOocc(nbasis+1:nmo) 4041 wfntype=4 4042end if 4043write(*,*) "Done! Basis function information now correspond to natural orbitals" 4044 4045write(*,"(/,a)") " If next you intend to analyze real space functions based on the NOs, you should export .molden file & 4046and then reload it, so that GTF information will also correspond to NOs" 4047write(*,*) "Would you like to do this immediately? (y/n)" 4048read(*,*) selectyn 4049if (selectyn=='y') then 4050 call outmolden("new.molden",10) 4051 write(*,*) "The NOs have been exported to NO.molden in current folder" 4052 call dealloall 4053 write(*,*) "Loading NO.molden..." 4054 call readmolden("new.molden",1) 4055 write(*,"(a)") " Loading finished, now you can use main function 0 to visualize NOs as isosurfaces" 4056end if 4057end subroutine 4058 4059