!================================================================= ! Below code is developed by Sobereva, 2013-Apr-11, adapted from Multiwfn 3.0 program ! Users of this code are asked to include this reference in their publications: Tian Lu, Feiwu Chen, J. Comp. Chem. 33, 580-592 (2012) !================================================================= module defvar type atomtype character*2 name !name of atom integer index !The index in periodic table, if ECP was used, charge will smaller than this value real*8 x,y,z,charge !Coordinate(Bohr) and charge of atoms end type type primtype integer center,functype !The No. of nuclei that the basis function centered on and function type real*8 exp !Exponent end type character*2 :: name2ind(0:109)=(/ "Bq","H ","He", & !X(number O) is ghost atom "Li","Be","B ","C ","N ","O ","F ","Ne", & !3~10 "Na","Mg","Al","Si","P ","S ","Cl","Ar", & !11~18 "K ","Ca","Sc","Ti","V ","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr", & !19~36 "Rb","Sr","Y ","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I ","Xe", & !37~54 "Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu", & !55~71 "Hf","Ta","W ","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn", & !72~86 "Fr","Ra","Ac","Th","Pa","U ","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr", & !87~103 "Rf","Db","Sg","Bh","Hs","Mt" /) !104~109 integer :: type2ix(35)=(/ 0,1,0,0, 2,0,0,1,1,0, 3,0,0,2,2,0,1,1,0,1, 0,0,0,0,0,1,1,1,1,2,2,2,3,3,4 /) integer :: type2iy(35)=(/ 0,0,1,0, 0,2,0,1,0,1, 0,3,0,1,0,2,2,0,1,1, 0,1,2,3,4,0,1,2,3,0,1,2,0,1,0 /) integer :: type2iz(35)=(/ 0,0,0,1, 0,0,2,0,1,1, 0,0,3,0,1,1,0,2,2,1, 4,3,2,1,0,3,2,1,0,2,1,0,1,0,0 /) integer :: nmo=0,nprims=0,ncenter=0 !Number of orbital/primitive basis function/nuclei integer wfntype !0/1/2/3/4 means R/U/ROHF /R/U-Post-HF wavefunction real*8 :: nelec=0,naelec=0,nbelec=0 type(atomtype),allocatable :: a(:) type(primtype),allocatable :: b(:) real*8,allocatable :: MOocc(:),MOene(:) !Occupation number & energy of orbital integer,allocatable :: MOtype(:) !The type of orbital, (alpha&beta)=0/alpha=1/beta=2, not read from .wfn directly real*8,allocatable :: CO(:,:) !Coefficient matrix of primitive basis functions, including both normalization and contraction coefficients. Row/column of CO denote MO/GTF respectively, in contrary to convention end module module function use defvar implicit real*8 (a-h,o-z) contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Calculate wavefunction value of all orbitals and its derivatives at given point !! istart and iend is the range of the orbitals will be calculated, to calculate all orbitals, use 1,nmo !! runtype=1: value =2: value+dx/y/z subroutine orbderv(runtype,istart,iend,x,y,z,wfnval,grad) real*8 x,y,z,wfnval(nmo) real*8,optional :: grad(3,nmo) integer runtype,istart,iend wfnval=0D0 if (present(grad)) grad=0D0 do j=1,nprims ix=type2ix(b(j)%functype) iy=type2iy(b(j)%functype) iz=type2iz(b(j)%functype) ep=b(j)%exp sftx=x-a(b(j)%center)%x sfty=y-a(b(j)%center)%y sftz=z-a(b(j)%center)%z rr=sftx**2+sfty**2+sftz**2 expterm=exp(-ep*rr) GTFval=sftx**ix *sfty**iy *sftz**iz *expterm !Calculate orbital wavefunction value do imo=istart,iend wfnval(imo)=wfnval(imo)+co(imo,j)*GTFval end do if (runtype==2) then !Calculate 1-order derivative for current GTF tx=0.0D0 ty=0.0D0 tz=0.0D0 if (ix/=0) tx=ix*sftx**(ix-1) if (iy/=0) ty=iy*sfty**(iy-1) if (iz/=0) tz=iz*sftz**(iz-1) GTFdx=sfty**iy *sftz**iz *expterm*(tx-2*ep*sftx**(ix+1)) GTFdy=sftx**ix *sftz**iz *expterm*(ty-2*ep*sfty**(iy+1)) GTFdz=sftx**ix *sfty**iy *expterm*(tz-2*ep*sftz**(iz+1)) !Calculate 1-order derivative for orbitals do imo=istart,iend grad(1,imo)=grad(1,imo)+co(imo,j)*GTFdx grad(2,imo)=grad(2,imo)+co(imo,j)*GTFdy grad(3,imo)=grad(3,imo)+co(imo,j)*GTFdz end do end if !end runtype>=2 end do end subroutine !--------Output electron density at a point real*8 function fdens(x,y,z) real*8 x,y,z,wfnval(nmo) call orbderv(1,1,nmo,x,y,z,wfnval) fdens=0.0D0 do i=1,nmo fdens=fdens+MOocc(i)*wfnval(i)**2 end do end function end module !!========== Main program =========== program wfnprop use function implicit real*8 (a-h,o-z) character c80*80 write(*,*) "Input .wfn file name" read(*,"(a)") c80 call readwfn(c80) do while(.true.) write(*,*) "Input the coordinate in Bohr e.g. 1.0,1.2,-3.5" read(*,*) x,y,z write(*,"(' The density is',f16.10,' e/Bohr^3',/)") fdens(x,y,z) end do end program !!!------------------------- Read .wfn file subroutine readwfn(name) use defvar implicit real*8 (a-h,o-z) CHARACTER(LEN=*) name character*80 wfntitle,c80tmp*80 open(10,file=name,access="sequential",status="old") read(10,"(a80)") wfntitle read(10,"(8x,i15,13x,i7,11x,i9)") nmo,nprims,ncenter allocate(a(ncenter)) allocate(b(nprims)) allocate(co(nmo,nprims)) allocate(MOocc(nmo)) allocate(MOene(nmo)) allocate(MOtype(nmo)) do i=1,ncenter read(10,"(a24,3f12.8,10x,f5.1)") c80tmp,a(i)%x,a(i)%y,a(i)%z,a(i)%charge read(c80tmp,*) a(i)%name call lc2uc(a(i)%name(1:1)) call uc2lc(a(i)%name(2:2)) do j=0,109 if (a(i)%name==name2ind(j)) then a(i)%index=j exit end if end do end do read(10,"(20x,20i3)") (b(i)%center,i=1,nprims) read(10,"(20x,20i3)") (b(i)%functype,i=1,nprims) read(10,"(10x,5D14.7)") (b(i)%exp,i=1,nprims) !Read orbitals do i=1,nmo read(10,"(35X,f12.7,15x,f12.6)") MOocc(i),MOene(i) read(10,"(5D16.8)") (co(i,j),j=1,nprims) ! note: row/column of CO denote MO/basis function respectively, in contrary to convention end do read(10,*) close (10) ! Determine type of wavefunction if (sum(MOocc)==2*nmo) then wfntype=0 !This is restricted wavefunction MOtype=0 else if (sum(MOocc)==nmo) then wfntype=1 !This is unrestricted wavefunction MOtype=1 !Set all MO is alpha do i=2,nmo !if nmo=1, i will be set to 2, and no errors will appear if (MOene(i) MOocc(i-1)) then MOtype(1:i-1)=1 MOtype(i:nmo)=2 exit end if end do end if else wfntype=2 !This is RO wavefunction MOtype=0 do i=1,nmo if (MOocc(i)==1) MOtype(i)=1 !alpha end do end if !Count electrons nelec=0 naelec=0 nbelec=0 do imo=1,nmo if (MOtype(imo)==0) then naelec=naelec+MOocc(imo)/2D0 nbelec=nbelec+MOocc(imo)/2D0 else if (MOtype(imo)==1) then naelec=naelec+MOocc(imo) else if (MOtype(imo)==2) then nbelec=nbelec+MOocc(imo) end if end do nelec=naelec+nbelec !Summary write(*,*) write(*,"('Total/Alpha/Beta electrons:',3f12.4)") nelec,naelec,nbelec write(*,"('Net charge:',f12.5,' Expected multiplicity:',i5)") sum(a(:)%charge)-nelec,nint(naelec-nbelec)+1 write(*,"('The number of orbitals:',i6,', Atoms:',i7,', GTFs:',i7)") nmo,ncenter,nprims if (wfntype==0) write(*,"('This is restricted close-shell single-determinant wavefunction')") if (wfntype==1) write(*,"('This is unrestricted single-determinant wavefunction')") if (wfntype==2) write(*,"('This is restricted open-shell wavefunction')") if (wfntype==3) write(*,"('This is close-shell post-HF wavefunction')") if (wfntype==4) write(*,"('This is open-shell post-HF wavefunction')") if (wfntype==1.or.wfntype==4) then do i=1,nmo if (MOtype(i)==2) exit end do write(*,"('Orbitals from 1 to',i6,' are alpha type, from',i6,' to',i6,' are beta type')") i-1,i,nmo end if write(*,"(/,'Title line of this file:',/,a79,/)") wfntitle end subroutine !!--------- Convert a character to lower case subroutine uc2lc(inc) character*1 inc itmp=ichar(inc) if (itmp>=65.and.itmp<=90) itmp=itmp+32 inc=char(itmp) end subroutine !!--------- Convert a character to upper case subroutine lc2uc(inc) character*1 inc itmp=ichar(inc) if (itmp>=97.and.itmp<=122) itmp=itmp-32 inc=char(itmp) end subroutine

map