From bb091d247569c5276006e9225c1282a4a7a2cc75 Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" Date: Sat, 29 Oct 2022 12:38:01 +0900 Subject: [PATCH] Merge orbderv into one routine --- .../Multiwfn_3.8_dev_src_Linux/define.f90 | 13 +- .../Multiwfn_3.8_dev_src_Linux/function.f90 | 862 ++---------------- 2 files changed, 102 insertions(+), 773 deletions(-) diff --git a/MultiWFN/Multiwfn_3.8_dev_src_Linux/define.f90 b/MultiWFN/Multiwfn_3.8_dev_src_Linux/define.f90 index ee7e6f5..b8917c7 100644 --- a/MultiWFN/Multiwfn_3.8_dev_src_Linux/define.f90 +++ b/MultiWFN/Multiwfn_3.8_dev_src_Linux/define.f90 @@ -257,23 +257,26 @@ real*8 :: totenergy=0,virialratio=2,nelec=0,naelec=0,nbelec=0 integer :: loadmulti=-99,loadcharge=-99 !Spin multiplicity and net charge, loaded directly from input file (e.g. from .gjf or .xyz), only utilized in rare cases. -99 means unloaded !-------- Variables for nuclei & GTF & Orbitals. Note: Row and column of CO(:,:) correspond to orbital and GTF, respectively, in contrary to convention type(atomtype),allocatable :: a(:),a_org(:),a_tmp(:) !a_tmp is only used in local temporary operation, should be destoried immediatedly after using -type(primtype),allocatable :: b(:),b_org(:),b_tmp(:) +type(primtype),allocatable,target :: b(:) +type(primtype),allocatable :: b_org(:),b_tmp(:) real*8,allocatable :: MOocc(:),MOocc_org(:),MOene(:),MOene_org(:) !Occupation number & energy of orbital integer,allocatable :: MOtype(:) !The type of orbitals, (alpha&beta)=0/alpha=1/beta=2, not read from .wfn directly character(len=10) :: orbtypename(0:2)=(/ character(len=10) :: "Alpha&Beta","Alpha","Beta" /) character(len=4),allocatable :: MOsym(:) !The symmetry of orbitals, meaningful when .mwfn/molden/gms is used -real*8,allocatable :: CO(:,:),CO_org(:,:),CO_tmp(:,:) !Coefficient matrix of primitive basis functions, including both normalization and contraction coefficients +real*8,allocatable, target :: CO(:,:) +real*8,allocatable :: CO_org(:,:),CO_tmp(:,:) !Coefficient matrix of primitive basis functions, including both normalization and contraction coefficients real*8,allocatable :: COtr(:,:) !Transposed CO matrix, which is used in some routines for faster calculation than using CO. Must be deallocated after using !Unique GTFs (the GTFs with identical center, type and exponent are combined together and leave only one). Can be activated after running gen_GTFuniq integer :: nprims_uniq=0 !0 means uninitialized -type(primtype),allocatable :: b_uniq(:) !b of unique GTFs -real*8,allocatable :: CO_uniq(:,:) !CO of b_uniq. The coefficients of duplicated GTFs are summed together +type(primtype),allocatable,target :: b_uniq(:) !b of unique GTFs +real*8,allocatable,target :: CO_uniq(:,:) !CO of b_uniq. The coefficients of duplicated GTFs are summed together !-------- Describe inner electron density in EDF section type(primtype),allocatable :: b_EDF(:) real*8,allocatable :: CO_EDF(:) integer :: nEDFprims=0,nEDFelec=0 !Electrons represented by EDF !-------- Promolecular wavefunction. Share same "a" and "b" -real*8,allocatable :: CO_pmol(:,:),MOocc_pmol(:),MOene_pmol(:) +real*8,allocatable,target :: CO_pmol(:,:) +real*8,allocatable :: MOocc_pmol(:),MOene_pmol(:) integer,allocatable :: MOtype_pmol(:) integer :: nmo_pmol=0 !-------- Variables when basis functions are basis rather than primitive function as basis diff --git a/MultiWFN/Multiwfn_3.8_dev_src_Linux/function.f90 b/MultiWFN/Multiwfn_3.8_dev_src_Linux/function.f90 index a6d448d..0fc0863 100644 --- a/MultiWFN/Multiwfn_3.8_dev_src_Linux/function.f90 +++ b/MultiWFN/Multiwfn_3.8_dev_src_Linux/function.f90 @@ -392,244 +392,96 @@ end function !! 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 =3: value+dxx/yy/zz(diagonal of hess) =4: value+dx/y/z+Hessian !! =5: value+dx/y/z+hess+3-order derivative tensor -subroutine orbderv(runtype,istart,iend,x,y,z,wfnval,grad,hess,tens3) +!! if promol == .true., compute for promolecular wavefunction +subroutine orbderv(runtype,istart,iend,x,y,z,wfnval,grad,hess,tens3,promol) real*8 x,y,z,wfnval(nmo) real*8,optional :: grad(3,nmo),hess(3,3,nmo),tens3(3,3,3,nmo) integer runtype,istart,iend -!real*8 GTFexpterm(nprims) - -if (ifPBC>0) then !Consider PBC - call orbderv_PBC(runtype,istart,iend,x,y,z,wfnval,grad,hess,tens3) - return -end if +logical,optional :: promol +logical :: promol_ +real*8 tvec(3), xmove,ymove,zmove +type(primtype), pointer :: b_p(:) +real*8, pointer :: CO_p(:,:) +integer :: nprims_tmp +real*8 :: cutoff_tmp wfnval=0D0 if (present(grad)) grad=0D0 if (present(hess)) hess=0D0 if (present(tens3)) tens3=0D0 + +promol_ = .false. +if (present(promol)) promol_ = promol + lastcen=-1 !Arbitrary value -!If the center/exp of current GTF is the same as previous, then we do not need to recalculate them -if (nprims_uniq==0) then !Unique GTF information is not available - do j=1,nprims - jcen=b(j)%center - if (jcen/=lastcen) then - sftx=x-a(jcen)%x - sfty=y-a(jcen)%y - sftz=z-a(jcen)%z - sftx2=sftx*sftx - sfty2=sfty*sfty - sftz2=sftz*sftz - rr=sftx2+sfty2+sftz2 - end if - - ep=b(j)%exp - tmpval=-ep*rr - lastcen=jcen - if (tmpval>expcutoff.or.expcutoff>0) then - expterm=exp(tmpval) - else - cycle - end if - - !Calculate value for current GTF - jtype=b(j)%type - ix=type2ix(jtype) - iy=type2iy(jtype) - iz=type2iz(jtype) - if (jtype==1) then !Some functype use manually optimized formula for cutting down computational time - GTFval=expterm - else if (jtype==2) then - GTFval=sftx*expterm - else if (jtype==3) then - GTFval=sfty*expterm - else if (jtype==4) then - GTFval=sftz*expterm - else if (jtype==5) then - GTFval=sftx2*expterm - else if (jtype==6) then - GTFval=sfty2*expterm - else if (jtype==7) then - GTFval=sftz2*expterm - else if (jtype==8) then - GTFval=sftx*sfty*expterm - else if (jtype==9) then - GTFval=sftx*sftz*expterm - else if (jtype==10) then - GTFval=sfty*sftz*expterm - else !If above conditions are not satisfied (angular moment higher than f), the function will be calculated explicitly - GTFval=sftx**ix *sfty**iy *sftz**iz *expterm - end if - !Calculate orbital wavefunction value. This is bottle neck of cost of this function - 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=0D0 - ty=0D0 - tz=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 +tvec = 0D0 +xmove=tvec(1) +ymove=tvec(2) +zmove=tvec(3) - if (runtype>=3) then - !Calculate 2-order derivative for current GTF - txx=0D0 - tyy=0D0 - tzz=0D0 - if (ix>=2) txx=ix*(ix-1)*sftx**(ix-2) - if (iy>=2) tyy=iy*(iy-1)*sfty**(iy-2) - if (iz>=2) tzz=iz*(iz-1)*sftz**(iz-2) - GTFdxx=sfty**iy *sftz**iz *expterm*( txx + 2*ep*sftx**ix*(-2*ix+2*ep*sftx2-1) ) - GTFdyy=sftx**ix *sftz**iz *expterm*( tyy + 2*ep*sfty**iy*(-2*iy+2*ep*sfty2-1) ) - GTFdzz=sftx**ix *sfty**iy *expterm*( tzz + 2*ep*sftz**iz*(-2*iz+2*ep*sftz2-1) ) - ttx=tx-2*ep*sftx**(ix+1) - tty=ty-2*ep*sfty**(iy+1) - ttz=tz-2*ep*sftz**(iz+1) - GTFdxy=sftz**iz *expterm*ttx*tty - GTFdyz=sftx**ix *expterm*tty*ttz - GTFdxz=sfty**iy *expterm*ttx*ttz - !Calculate diagonal Hessian elements for orbitals - do imo=istart,iend - hess(1,1,imo)=hess(1,1,imo)+CO(imo,j)*GTFdxx !dxx - hess(2,2,imo)=hess(2,2,imo)+CO(imo,j)*GTFdyy !dyy - hess(3,3,imo)=hess(3,3,imo)+CO(imo,j)*GTFdzz !dzz - end do - if (runtype>=4) then !Also process nondiagonal elements - do imo=istart,iend - hess(1,2,imo)=hess(1,2,imo)+CO(imo,j)*GTFdxy !dxy - hess(2,3,imo)=hess(2,3,imo)+CO(imo,j)*GTFdyz !dyz - hess(1,3,imo)=hess(1,3,imo)+CO(imo,j)*GTFdxz !dxz - end do - hess(2,1,:)=hess(1,2,:) - hess(3,2,:)=hess(2,3,:) - hess(3,1,:)=hess(1,3,:) - end if - - if (runtype>=5) then - !Calculate 3-order derivative for current GTF - ep2=ep*2D0 - ep4=ep*4D0 - epep4=ep2*ep2 - epep8=epep4*2D0 - !dxyz - a1=0D0 - b1=0D0 - c1=0D0 - if (ix>=1) a1=ix*sftx**(ix-1) - if (iy>=1) b1=iy*sfty**(iy-1) - if (iz>=1) c1=iz*sftz**(iz-1) - a2=-ep2*sftx**(ix+1) - b2=-ep2*sfty**(iy+1) - c2=-ep2*sftz**(iz+1) - GTFdxyz=(a1+a2)*(b1+b2)*(c1+c2)*expterm - !dxyy,dxxy,dxxz,dxzz,dyzz,dyyz - atmp=0D0 - btmp=0D0 - ctmp=0D0 - if (ix>=2) atmp=ix*(ix-1)*sftx**(ix-2) - if (iy>=2) btmp=iy*(iy-1)*sfty**(iy-2) - if (iz>=2) ctmp=iz*(iz-1)*sftz**(iz-2) - GTFdxyy=(a1+a2)*sftz**iz *expterm*(-ep4*iy*sfty**iy+epep4*sfty**(iy+2)+btmp-ep2*sfty**iy) - GTFdxxy=(b1+b2)*sftz**iz *expterm*(-ep4*ix*sftx**ix+epep4*sftx**(ix+2)+atmp-ep2*sftx**ix) !=dyxx - GTFdxxz=(c1+c2)*sfty**iy *expterm*(-ep4*ix*sftx**ix+epep4*sftx**(ix+2)+atmp-ep2*sftx**ix) !=dzxx - GTFdxzz=(a1+a2)*sfty**iy *expterm*(-ep4*iz*sftz**iz+epep4*sftz**(iz+2)+ctmp-ep2*sftz**iz) - GTFdyzz=(b1+b2)*sftx**ix *expterm*(-ep4*iz*sftz**iz+epep4*sftz**(iz+2)+ctmp-ep2*sftz**iz) - GTFdyyz=(c1+c2)*sftx**ix *expterm*(-ep4*iy*sfty**iy+epep4*sfty**(iy+2)+btmp-ep2*sfty**iy) !=dzyy - !dxxx,dyyy,dzzz - aatmp1=0D0 - bbtmp1=0D0 - cctmp1=0D0 - if (ix>=1) aatmp1=ep2*ix*sftx**(ix-1) - if (iy>=1) bbtmp1=ep2*iy*sfty**(iy-1) - if (iz>=1) cctmp1=ep2*iz*sftz**(iz-1) - aatmp2=0D0 - bbtmp2=0D0 - cctmp2=0D0 - if (ix>=2) aatmp2=ep2*ix*(ix-1)*sftx**(ix-1) - if (iy>=2) bbtmp2=ep2*iy*(iy-1)*sfty**(iy-1) - if (iz>=2) cctmp2=ep2*iz*(iz-1)*sftz**(iz-1) - aatmp3=0D0 - bbtmp3=0D0 - cctmp3=0D0 - if (ix>=3) aatmp3=ix*(ix-1)*(ix-2)*sftx**(ix-3) - if (iy>=3) bbtmp3=iy*(iy-1)*(iy-2)*sfty**(iy-3) - if (iz>=3) cctmp3=iz*(iz-1)*(iz-2)*sftz**(iz-3) - GTFdxxx=sfty**iy*sftz**iz*expterm*( (-2*ix+ep2*sftx2-1)*(-epep4*sftx**(ix+1) + aatmp1) - aatmp2 + epep8*sftx**(ix+1) + aatmp3 ) - GTFdyyy=sftx**ix*sftz**iz*expterm*( (-2*iy+ep2*sfty2-1)*(-epep4*sfty**(iy+1) + bbtmp1) - bbtmp2 + epep8*sfty**(iy+1) + bbtmp3 ) - GTFdzzz=sfty**iy*sftx**ix*expterm*( (-2*iz+ep2*sftz2-1)*(-epep4*sftz**(iz+1) + cctmp1) - cctmp2 + epep8*sftz**(iz+1) + cctmp3 ) - - !Calculate 3-order derivative tensor for orbital wavefunction - do imo=istart,iend - tens3(1,1,1,imo)=tens3(1,1,1,imo)+CO(imo,j)*GTFdxxx !dxxx - tens3(2,2,2,imo)=tens3(2,2,2,imo)+CO(imo,j)*GTFdyyy !dyyy - tens3(3,3,3,imo)=tens3(3,3,3,imo)+CO(imo,j)*GTFdzzz !dzzz - tens3(1,2,2,imo)=tens3(1,2,2,imo)+CO(imo,j)*GTFdxyy !dxyy - tens3(1,1,2,imo)=tens3(1,1,2,imo)+CO(imo,j)*GTFdxxy !dxxy - tens3(1,1,3,imo)=tens3(1,1,3,imo)+CO(imo,j)*GTFdxxz !dxxz - tens3(1,3,3,imo)=tens3(1,3,3,imo)+CO(imo,j)*GTFdxzz !dxzz - tens3(2,3,3,imo)=tens3(2,3,3,imo)+CO(imo,j)*GTFdyzz !dyzz - tens3(2,2,3,imo)=tens3(2,2,3,imo)+CO(imo,j)*GTFdyyz !dyyz - tens3(1,2,3,imo)=tens3(1,2,3,imo)+CO(imo,j)*GTFdxyz !dxyz - end do - tens3(1,2,1,:)=tens3(1,1,2,:) !dxyx=dxxy - tens3(1,3,1,:)=tens3(1,1,3,:) !dxzx=dxxz - tens3(1,3,2,:)=tens3(1,2,3,:) !dxzy=dxyz - tens3(2,1,1,:)=tens3(1,1,2,:) !dyxx=dxxy - tens3(2,1,2,:)=tens3(1,2,2,:) !dyxy=dxyy - tens3(2,1,3,:)=tens3(1,2,3,:) !dyxz=dxyz - tens3(2,2,1,:)=tens3(1,2,2,:) !dyyx=dxyy - tens3(2,3,1,:)=tens3(1,2,3,:) !dyzx=dxyz - tens3(2,3,2,:)=tens3(2,2,3,:) !dyzy=dyyz - tens3(3,1,1,:)=tens3(1,1,3,:) !dzxx=dxxz - tens3(3,1,2,:)=tens3(1,2,3,:) !dzxy=dxyz - tens3(3,1,3,:)=tens3(1,3,3,:) !dzxz=dxzz - tens3(3,2,1,:)=tens3(1,2,3,:) !dzyx=dxyz - tens3(3,2,2,:)=tens3(2,2,3,:) !dzyy=dyyz - tens3(3,2,3,:)=tens3(2,3,3,:) !dzyz=dyzz - tens3(3,3,1,:)=tens3(1,3,3,:) !dzzx=dxzz - tens3(3,3,2,:)=tens3(2,3,3,:) !dzzy=dyzz - end if !end runtype>=5 - - end if !end runtype>=3 - end if !end runtype>=2 - end do +!calculate orbderv for promolecular wavefunction +if (promol_) then + nprims_tmp = nprims + b_p => b + CO_p => CO_pmol +!If the center/exp of current GTF is the same as previous, then we do not need to recalculate them +else if (nprims_uniq==0) then !Unique GTF information is not available + nprims_tmp = nprims + b_p => b + CO_p => CO else !Unique GTF information has been constructed by gen_GTFuniq - do j=1,nprims_uniq - jcen=b_uniq(j)%center + nprims_tmp = nprims_uniq + b_p => b_uniq + CO_p => CO_uniq +end if + +if (ifPBC>0) then !Consider PBC + cutoff_tmp = expcutoff_PBC + call getpointcell(x,y,z,ic,jc,kc) + do icell=ic-PBCnx,ic+PBCnx + do jcell=jc-PBCny,jc+PBCny + do kcell=kc-PBCnz,kc+PBCnz + lastcen=-1 !Arbitrary value + call tvec_PBC(icell,jcell,kcell,tvec) + xmove=tvec(1) + ymove=tvec(2) + zmove=tvec(3) + call orbderv_eval() + end do + end do + end do +else + cutoff_tmp = expcutoff + call orbderv_eval() +end if + +contains + subroutine orbderv_eval() + do j=1,nprims_tmp + jcen=b_p(j)%center if (jcen/=lastcen) then - sftx=x-a(jcen)%x - sfty=y-a(jcen)%y - sftz=z-a(jcen)%z + sftx=x-(a(jcen)%x+xmove) + sfty=y-(a(jcen)%y+ymove) + sftz=z-(a(jcen)%z+zmove) sftx2=sftx*sftx sfty2=sfty*sfty sftz2=sftz*sftz rr=sftx2+sfty2+sftz2 end if - - ep=b_uniq(j)%exp + + ep=b_p(j)%exp tmpval=-ep*rr lastcen=jcen - if (tmpval>expcutoff.or.expcutoff>0) then + if (tmpval>cutoff_tmp.or.cutoff_tmp>0) then expterm=exp(tmpval) else cycle end if - + !Calculate value for current GTF - jtype=b_uniq(j)%type + jtype=b_p(j)%type ix=type2ix(jtype) iy=type2iy(jtype) iz=type2iz(jtype) @@ -658,9 +510,9 @@ else !Unique GTF information has been constructed by gen_GTFuniq end if !Calculate orbital wavefunction value. This is bottle neck of cost of this function do imo=istart,iend - wfnval(imo)=wfnval(imo)+CO_uniq(imo,j)*GTFval + wfnval(imo)=wfnval(imo)+CO_p(imo,j)*GTFval end do - + if (runtype>=2) then !Calculate 1-order derivative for current GTF tx=0D0 @@ -674,9 +526,9 @@ else !Unique GTF information has been constructed by gen_GTFuniq 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_uniq(imo,j)*GTFdx - grad(2,imo)=grad(2,imo)+CO_uniq(imo,j)*GTFdy - grad(3,imo)=grad(3,imo)+CO_uniq(imo,j)*GTFdz + grad(1,imo)=grad(1,imo)+CO_p(imo,j)*GTFdx + grad(2,imo)=grad(2,imo)+CO_p(imo,j)*GTFdy + grad(3,imo)=grad(3,imo)+CO_p(imo,j)*GTFdz end do if (runtype>=3) then @@ -698,21 +550,22 @@ else !Unique GTF information has been constructed by gen_GTFuniq GTFdxz=sfty**iy *expterm*ttx*ttz !Calculate diagonal Hessian elements for orbitals do imo=istart,iend - hess(1,1,imo)=hess(1,1,imo)+CO_uniq(imo,j)*GTFdxx !dxx - hess(2,2,imo)=hess(2,2,imo)+CO_uniq(imo,j)*GTFdyy !dyy - hess(3,3,imo)=hess(3,3,imo)+CO_uniq(imo,j)*GTFdzz !dzz + hess(1,1,imo)=hess(1,1,imo)+CO_p(imo,j)*GTFdxx !dxx + hess(2,2,imo)=hess(2,2,imo)+CO_p(imo,j)*GTFdyy !dyy + hess(3,3,imo)=hess(3,3,imo)+CO_p(imo,j)*GTFdzz !dzz end do + if (runtype>=4) then !Also process nondiagonal elements do imo=istart,iend - hess(1,2,imo)=hess(1,2,imo)+CO_uniq(imo,j)*GTFdxy !dxy - hess(2,3,imo)=hess(2,3,imo)+CO_uniq(imo,j)*GTFdyz !dyz - hess(1,3,imo)=hess(1,3,imo)+CO_uniq(imo,j)*GTFdxz !dxz + hess(1,2,imo)=hess(1,2,imo)+CO_p(imo,j)*GTFdxy !dxy + hess(2,3,imo)=hess(2,3,imo)+CO_p(imo,j)*GTFdyz !dyz + hess(1,3,imo)=hess(1,3,imo)+CO_p(imo,j)*GTFdxz !dxz end do hess(2,1,:)=hess(1,2,:) hess(3,2,:)=hess(2,3,:) hess(3,1,:)=hess(1,3,:) end if - + if (runtype>=5) then !Calculate 3-order derivative for current GTF ep2=ep*2D0 @@ -765,20 +618,21 @@ else !Unique GTF information has been constructed by gen_GTFuniq GTFdxxx=sfty**iy*sftz**iz*expterm*( (-2*ix+ep2*sftx2-1)*(-epep4*sftx**(ix+1) + aatmp1) - aatmp2 + epep8*sftx**(ix+1) + aatmp3 ) GTFdyyy=sftx**ix*sftz**iz*expterm*( (-2*iy+ep2*sfty2-1)*(-epep4*sfty**(iy+1) + bbtmp1) - bbtmp2 + epep8*sfty**(iy+1) + bbtmp3 ) GTFdzzz=sfty**iy*sftx**ix*expterm*( (-2*iz+ep2*sftz2-1)*(-epep4*sftz**(iz+1) + cctmp1) - cctmp2 + epep8*sftz**(iz+1) + cctmp3 ) - + !Calculate 3-order derivative tensor for orbital wavefunction do imo=istart,iend - tens3(1,1,1,imo)=tens3(1,1,1,imo)+CO_uniq(imo,j)*GTFdxxx !dxxx - tens3(2,2,2,imo)=tens3(2,2,2,imo)+CO_uniq(imo,j)*GTFdyyy !dyyy - tens3(3,3,3,imo)=tens3(3,3,3,imo)+CO_uniq(imo,j)*GTFdzzz !dzzz - tens3(1,2,2,imo)=tens3(1,2,2,imo)+CO_uniq(imo,j)*GTFdxyy !dxyy - tens3(1,1,2,imo)=tens3(1,1,2,imo)+CO_uniq(imo,j)*GTFdxxy !dxxy - tens3(1,1,3,imo)=tens3(1,1,3,imo)+CO_uniq(imo,j)*GTFdxxz !dxxz - tens3(1,3,3,imo)=tens3(1,3,3,imo)+CO_uniq(imo,j)*GTFdxzz !dxzz - tens3(2,3,3,imo)=tens3(2,3,3,imo)+CO_uniq(imo,j)*GTFdyzz !dyzz - tens3(2,2,3,imo)=tens3(2,2,3,imo)+CO_uniq(imo,j)*GTFdyyz !dyyz - tens3(1,2,3,imo)=tens3(1,2,3,imo)+CO_uniq(imo,j)*GTFdxyz !dxyz + tens3(1,1,1,imo)=tens3(1,1,1,imo)+CO_p(imo,j)*GTFdxxx !dxxx + tens3(2,2,2,imo)=tens3(2,2,2,imo)+CO_p(imo,j)*GTFdyyy !dyyy + tens3(3,3,3,imo)=tens3(3,3,3,imo)+CO_p(imo,j)*GTFdzzz !dzzz + tens3(1,2,2,imo)=tens3(1,2,2,imo)+CO_p(imo,j)*GTFdxyy !dxyy + tens3(1,1,2,imo)=tens3(1,1,2,imo)+CO_p(imo,j)*GTFdxxy !dxxy + tens3(1,1,3,imo)=tens3(1,1,3,imo)+CO_p(imo,j)*GTFdxxz !dxxz + tens3(1,3,3,imo)=tens3(1,3,3,imo)+CO_p(imo,j)*GTFdxzz !dxzz + tens3(2,3,3,imo)=tens3(2,3,3,imo)+CO_p(imo,j)*GTFdyzz !dyzz + tens3(2,2,3,imo)=tens3(2,2,3,imo)+CO_p(imo,j)*GTFdyyz !dyyz + tens3(1,2,3,imo)=tens3(1,2,3,imo)+CO_p(imo,j)*GTFdxyz !dxyz end do + tens3(1,2,1,:)=tens3(1,1,2,:) !dxyx=dxxy tens3(1,3,1,:)=tens3(1,1,3,:) !dxzx=dxxz tens3(1,3,2,:)=tens3(1,2,3,:) !dxzy=dxyz @@ -797,539 +651,11 @@ else !Unique GTF information has been constructed by gen_GTFuniq tens3(3,3,1,:)=tens3(1,3,3,:) !dzzx=dxzz tens3(3,3,2,:)=tens3(2,3,3,:) !dzzy=dyzz end if !end runtype>=5 - + end if !end runtype>=3 end if !end runtype>=2 end do -end if -end subroutine - - - - -!!--------- The same as orbderv, but explicitly consider PBC -!! Calculate wavefunction value of a range of orbitals and their derivatives at a given point, up to third-order -!! 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 =3: value+dxx/yy/zz(diagonal of hess) =4: value+dx/y/z+Hessian -!! =5: value+dx/y/z+hess+3-order derivative tensor -subroutine orbderv_PBC(runtype,istart,iend,x,y,z,wfnval,grad,hess,tens3) -real*8 x,y,z,wfnval(nmo),tvec(3) -real*8,optional :: grad(3,nmo),hess(3,3,nmo),tens3(3,3,3,nmo) -integer runtype,istart,iend -real*8 GTFexpterm(nprims) - -wfnval=0D0 -if (present(grad)) grad=0D0 -if (present(hess)) hess=0D0 -if (present(tens3)) tens3=0D0 - -call getpointcell(x,y,z,ic,jc,kc) -do icell=ic-PBCnx,ic+PBCnx - do jcell=jc-PBCny,jc+PBCny - do kcell=kc-PBCnz,kc+PBCnz - lastcen=-1 !Arbitrary value - call tvec_PBC(icell,jcell,kcell,tvec) - xmove=tvec(1) - ymove=tvec(2) - zmove=tvec(3) - - if (nprims_uniq==0) then !Unique GTF information is not available - do j=1,nprims - jcen=b(j)%center - if (jcen/=lastcen) then - sftx=x-(a(jcen)%x+xmove) - sfty=y-(a(jcen)%y+ymove) - sftz=z-(a(jcen)%z+zmove) - sftx2=sftx*sftx - sfty2=sfty*sfty - sftz2=sftz*sftz - rr=sftx2+sfty2+sftz2 - end if - ep=b(j)%exp - tmpval=-ep*rr - lastcen=jcen - if (tmpval>expcutoff_PBC.or.expcutoff_PBC>0) then - expterm=exp(tmpval) - else - cycle - end if - - !Calculate value for current GTF - jtype=b(j)%type - ix=type2ix(jtype) - iy=type2iy(jtype) - iz=type2iz(jtype) - if (jtype==1) then - GTFval=expterm - else if (jtype==2) then - GTFval=sftx*expterm - else if (jtype==3) then - GTFval=sfty*expterm - else if (jtype==4) then - GTFval=sftz*expterm - else if (jtype==5) then - GTFval=sftx2*expterm - else if (jtype==6) then - GTFval=sfty2*expterm - else if (jtype==7) then - GTFval=sftz2*expterm - else if (jtype==8) then - GTFval=sftx*sfty*expterm - else if (jtype==9) then - GTFval=sftx*sftz*expterm - else if (jtype==10) then - GTFval=sfty*sftz*expterm - else - GTFval=sftx**ix *sfty**iy *sftz**iz *expterm - end if - !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=0D0 - ty=0D0 - tz=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 - - if (runtype>=3) then - !Calculate 2-order derivative for current GTF - txx=0D0 - tyy=0D0 - tzz=0D0 - if (ix>=2) txx=ix*(ix-1)*sftx**(ix-2) - if (iy>=2) tyy=iy*(iy-1)*sfty**(iy-2) - if (iz>=2) tzz=iz*(iz-1)*sftz**(iz-2) - GTFdxx=sfty**iy *sftz**iz *expterm*( txx + 2*ep*sftx**ix*(-2*ix+2*ep*sftx2-1) ) - GTFdyy=sftx**ix *sftz**iz *expterm*( tyy + 2*ep*sfty**iy*(-2*iy+2*ep*sfty2-1) ) - GTFdzz=sftx**ix *sfty**iy *expterm*( tzz + 2*ep*sftz**iz*(-2*iz+2*ep*sftz2-1) ) - ttx=tx-2*ep*sftx**(ix+1) - tty=ty-2*ep*sfty**(iy+1) - ttz=tz-2*ep*sftz**(iz+1) - GTFdxy=sftz**iz *expterm*ttx*tty - GTFdyz=sftx**ix *expterm*tty*ttz - GTFdxz=sfty**iy *expterm*ttx*ttz - !Calculate diagonal Hessian elements for orbitals - do imo=istart,iend - hess(1,1,imo)=hess(1,1,imo)+CO(imo,j)*GTFdxx !dxx - hess(2,2,imo)=hess(2,2,imo)+CO(imo,j)*GTFdyy !dyy - hess(3,3,imo)=hess(3,3,imo)+CO(imo,j)*GTFdzz !dzz - end do - if (runtype>=4) then !Also process nondiagonal elements - do imo=istart,iend - hess(1,2,imo)=hess(1,2,imo)+CO(imo,j)*GTFdxy !dxy - hess(2,3,imo)=hess(2,3,imo)+CO(imo,j)*GTFdyz !dyz - hess(1,3,imo)=hess(1,3,imo)+CO(imo,j)*GTFdxz !dxz - end do - hess(2,1,:)=hess(1,2,:) - hess(3,2,:)=hess(2,3,:) - hess(3,1,:)=hess(1,3,:) - end if - - if (runtype>=5) then - !Calculate 3-order derivative for current GTF - ep2=ep*2D0 - ep4=ep*4D0 - epep4=ep2*ep2 - epep8=epep4*2D0 - !dxyz - a1=0D0 - b1=0D0 - c1=0D0 - if (ix>=1) a1=ix*sftx**(ix-1) - if (iy>=1) b1=iy*sfty**(iy-1) - if (iz>=1) c1=iz*sftz**(iz-1) - a2=-ep2*sftx**(ix+1) - b2=-ep2*sfty**(iy+1) - c2=-ep2*sftz**(iz+1) - GTFdxyz=(a1+a2)*(b1+b2)*(c1+c2)*expterm - !dxyy,dxxy,dxxz,dxzz,dyzz,dyyz - atmp=0D0 - btmp=0D0 - ctmp=0D0 - if (ix>=2) atmp=ix*(ix-1)*sftx**(ix-2) - if (iy>=2) btmp=iy*(iy-1)*sfty**(iy-2) - if (iz>=2) ctmp=iz*(iz-1)*sftz**(iz-2) - GTFdxyy=(a1+a2)*sftz**iz *expterm*(-ep4*iy*sfty**iy+epep4*sfty**(iy+2)+btmp-ep2*sfty**iy) - GTFdxxy=(b1+b2)*sftz**iz *expterm*(-ep4*ix*sftx**ix+epep4*sftx**(ix+2)+atmp-ep2*sftx**ix) !=dyxx - GTFdxxz=(c1+c2)*sfty**iy *expterm*(-ep4*ix*sftx**ix+epep4*sftx**(ix+2)+atmp-ep2*sftx**ix) !=dzxx - GTFdxzz=(a1+a2)*sfty**iy *expterm*(-ep4*iz*sftz**iz+epep4*sftz**(iz+2)+ctmp-ep2*sftz**iz) - GTFdyzz=(b1+b2)*sftx**ix *expterm*(-ep4*iz*sftz**iz+epep4*sftz**(iz+2)+ctmp-ep2*sftz**iz) - GTFdyyz=(c1+c2)*sftx**ix *expterm*(-ep4*iy*sfty**iy+epep4*sfty**(iy+2)+btmp-ep2*sfty**iy) !=dzyy - !dxxx,dyyy,dzzz - aatmp1=0D0 - bbtmp1=0D0 - cctmp1=0D0 - if (ix>=1) aatmp1=ep2*ix*sftx**(ix-1) - if (iy>=1) bbtmp1=ep2*iy*sfty**(iy-1) - if (iz>=1) cctmp1=ep2*iz*sftz**(iz-1) - aatmp2=0D0 - bbtmp2=0D0 - cctmp2=0D0 - if (ix>=2) aatmp2=ep2*ix*(ix-1)*sftx**(ix-1) - if (iy>=2) bbtmp2=ep2*iy*(iy-1)*sfty**(iy-1) - if (iz>=2) cctmp2=ep2*iz*(iz-1)*sftz**(iz-1) - aatmp3=0D0 - bbtmp3=0D0 - cctmp3=0D0 - if (ix>=3) aatmp3=ix*(ix-1)*(ix-2)*sftx**(ix-3) - if (iy>=3) bbtmp3=iy*(iy-1)*(iy-2)*sfty**(iy-3) - if (iz>=3) cctmp3=iz*(iz-1)*(iz-2)*sftz**(iz-3) - GTFdxxx=sfty**iy*sftz**iz*expterm*( (-2*ix+ep2*sftx2-1)*(-epep4*sftx**(ix+1) + aatmp1) - aatmp2 + epep8*sftx**(ix+1) + aatmp3 ) - GTFdyyy=sftx**ix*sftz**iz*expterm*( (-2*iy+ep2*sfty2-1)*(-epep4*sfty**(iy+1) + bbtmp1) - bbtmp2 + epep8*sfty**(iy+1) + bbtmp3 ) - GTFdzzz=sfty**iy*sftx**ix*expterm*( (-2*iz+ep2*sftz2-1)*(-epep4*sftz**(iz+1) + cctmp1) - cctmp2 + epep8*sftz**(iz+1) + cctmp3 ) - - !Calculate 3-order derivative tensor for orbital wavefunction - do imo=istart,iend - tens3(1,1,1,imo)=tens3(1,1,1,imo)+CO(imo,j)*GTFdxxx !dxxx - tens3(2,2,2,imo)=tens3(2,2,2,imo)+CO(imo,j)*GTFdyyy !dyyy - tens3(3,3,3,imo)=tens3(3,3,3,imo)+CO(imo,j)*GTFdzzz !dzzz - tens3(1,2,2,imo)=tens3(1,2,2,imo)+CO(imo,j)*GTFdxyy !dxyy - tens3(1,1,2,imo)=tens3(1,1,2,imo)+CO(imo,j)*GTFdxxy !dxxy - tens3(1,1,3,imo)=tens3(1,1,3,imo)+CO(imo,j)*GTFdxxz !dxxz - tens3(1,3,3,imo)=tens3(1,3,3,imo)+CO(imo,j)*GTFdxzz !dxzz - tens3(2,3,3,imo)=tens3(2,3,3,imo)+CO(imo,j)*GTFdyzz !dyzz - tens3(2,2,3,imo)=tens3(2,2,3,imo)+CO(imo,j)*GTFdyyz !dyyz - tens3(1,2,3,imo)=tens3(1,2,3,imo)+CO(imo,j)*GTFdxyz !dxyz - end do - tens3(1,2,1,:)=tens3(1,1,2,:) !dxyx=dxxy - tens3(1,3,1,:)=tens3(1,1,3,:) !dxzx=dxxz - tens3(1,3,2,:)=tens3(1,2,3,:) !dxzy=dxyz - tens3(2,1,1,:)=tens3(1,1,2,:) !dyxx=dxxy - tens3(2,1,2,:)=tens3(1,2,2,:) !dyxy=dxyy - tens3(2,1,3,:)=tens3(1,2,3,:) !dyxz=dxyz - tens3(2,2,1,:)=tens3(1,2,2,:) !dyyx=dxyy - tens3(2,3,1,:)=tens3(1,2,3,:) !dyzx=dxyz - tens3(2,3,2,:)=tens3(2,2,3,:) !dyzy=dyyz - tens3(3,1,1,:)=tens3(1,1,3,:) !dzxx=dxxz - tens3(3,1,2,:)=tens3(1,2,3,:) !dzxy=dxyz - tens3(3,1,3,:)=tens3(1,3,3,:) !dzxz=dxzz - tens3(3,2,1,:)=tens3(1,2,3,:) !dzyx=dxyz - tens3(3,2,2,:)=tens3(2,2,3,:) !dzyy=dyyz - tens3(3,2,3,:)=tens3(2,3,3,:) !dzyz=dyzz - tens3(3,3,1,:)=tens3(1,3,3,:) !dzzx=dxzz - tens3(3,3,2,:)=tens3(2,3,3,:) !dzzy=dyzz - end if !end runtype>=5 - - end if !end runtype>=3 - end if !end runtype>=2 - end do - - else !Unique GTF information has been constructed by gen_GTFuniq - do j=1,nprims_uniq - jcen=b_uniq(j)%center - if (jcen/=lastcen) then - sftx=x-(a(jcen)%x+xmove) - sfty=y-(a(jcen)%y+ymove) - sftz=z-(a(jcen)%z+zmove) - sftx2=sftx*sftx - sfty2=sfty*sfty - sftz2=sftz*sftz - rr=sftx2+sfty2+sftz2 - end if - ep=b_uniq(j)%exp - tmpval=-ep*rr - lastcen=jcen - if (tmpval>expcutoff_PBC.or.expcutoff_PBC>0) then - expterm=exp(tmpval) - else - cycle - end if - - !Calculate value for current GTF - jtype=b_uniq(j)%type - ix=type2ix(jtype) - iy=type2iy(jtype) - iz=type2iz(jtype) - if (jtype==1) then - GTFval=expterm - else if (jtype==2) then - GTFval=sftx*expterm - else if (jtype==3) then - GTFval=sfty*expterm - else if (jtype==4) then - GTFval=sftz*expterm - else if (jtype==5) then - GTFval=sftx2*expterm - else if (jtype==6) then - GTFval=sfty2*expterm - else if (jtype==7) then - GTFval=sftz2*expterm - else if (jtype==8) then - GTFval=sftx*sfty*expterm - else if (jtype==9) then - GTFval=sftx*sftz*expterm - else if (jtype==10) then - GTFval=sfty*sftz*expterm - else - GTFval=sftx**ix *sfty**iy *sftz**iz *expterm - end if - !Calculate orbital wavefunction value - do imo=istart,iend - wfnval(imo)=wfnval(imo)+CO_uniq(imo,j)*GTFval - end do - - if (runtype>=2) then - !Calculate 1-order derivative for current GTF - tx=0D0 - ty=0D0 - tz=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_uniq(imo,j)*GTFdx - grad(2,imo)=grad(2,imo)+CO_uniq(imo,j)*GTFdy - grad(3,imo)=grad(3,imo)+CO_uniq(imo,j)*GTFdz - end do - - if (runtype>=3) then - !Calculate 2-order derivative for current GTF - txx=0D0 - tyy=0D0 - tzz=0D0 - if (ix>=2) txx=ix*(ix-1)*sftx**(ix-2) - if (iy>=2) tyy=iy*(iy-1)*sfty**(iy-2) - if (iz>=2) tzz=iz*(iz-1)*sftz**(iz-2) - GTFdxx=sfty**iy *sftz**iz *expterm*( txx + 2*ep*sftx**ix*(-2*ix+2*ep*sftx2-1) ) - GTFdyy=sftx**ix *sftz**iz *expterm*( tyy + 2*ep*sfty**iy*(-2*iy+2*ep*sfty2-1) ) - GTFdzz=sftx**ix *sfty**iy *expterm*( tzz + 2*ep*sftz**iz*(-2*iz+2*ep*sftz2-1) ) - ttx=tx-2*ep*sftx**(ix+1) - tty=ty-2*ep*sfty**(iy+1) - ttz=tz-2*ep*sftz**(iz+1) - GTFdxy=sftz**iz *expterm*ttx*tty - GTFdyz=sftx**ix *expterm*tty*ttz - GTFdxz=sfty**iy *expterm*ttx*ttz - !Calculate diagonal Hessian elements for orbitals - do imo=istart,iend - hess(1,1,imo)=hess(1,1,imo)+CO_uniq(imo,j)*GTFdxx !dxx - hess(2,2,imo)=hess(2,2,imo)+CO_uniq(imo,j)*GTFdyy !dyy - hess(3,3,imo)=hess(3,3,imo)+CO_uniq(imo,j)*GTFdzz !dzz - end do - if (runtype>=4) then !Also process nondiagonal elements - do imo=istart,iend - hess(1,2,imo)=hess(1,2,imo)+CO_uniq(imo,j)*GTFdxy !dxy - hess(2,3,imo)=hess(2,3,imo)+CO_uniq(imo,j)*GTFdyz !dyz - hess(1,3,imo)=hess(1,3,imo)+CO_uniq(imo,j)*GTFdxz !dxz - end do - hess(2,1,:)=hess(1,2,:) - hess(3,2,:)=hess(2,3,:) - hess(3,1,:)=hess(1,3,:) - end if - - if (runtype>=5) then - !Calculate 3-order derivative for current GTF - ep2=ep*2D0 - ep4=ep*4D0 - epep4=ep2*ep2 - epep8=epep4*2D0 - !dxyz - a1=0D0 - b1=0D0 - c1=0D0 - if (ix>=1) a1=ix*sftx**(ix-1) - if (iy>=1) b1=iy*sfty**(iy-1) - if (iz>=1) c1=iz*sftz**(iz-1) - a2=-ep2*sftx**(ix+1) - b2=-ep2*sfty**(iy+1) - c2=-ep2*sftz**(iz+1) - GTFdxyz=(a1+a2)*(b1+b2)*(c1+c2)*expterm - !dxyy,dxxy,dxxz,dxzz,dyzz,dyyz - atmp=0D0 - btmp=0D0 - ctmp=0D0 - if (ix>=2) atmp=ix*(ix-1)*sftx**(ix-2) - if (iy>=2) btmp=iy*(iy-1)*sfty**(iy-2) - if (iz>=2) ctmp=iz*(iz-1)*sftz**(iz-2) - GTFdxyy=(a1+a2)*sftz**iz *expterm*(-ep4*iy*sfty**iy+epep4*sfty**(iy+2)+btmp-ep2*sfty**iy) - GTFdxxy=(b1+b2)*sftz**iz *expterm*(-ep4*ix*sftx**ix+epep4*sftx**(ix+2)+atmp-ep2*sftx**ix) !=dyxx - GTFdxxz=(c1+c2)*sfty**iy *expterm*(-ep4*ix*sftx**ix+epep4*sftx**(ix+2)+atmp-ep2*sftx**ix) !=dzxx - GTFdxzz=(a1+a2)*sfty**iy *expterm*(-ep4*iz*sftz**iz+epep4*sftz**(iz+2)+ctmp-ep2*sftz**iz) - GTFdyzz=(b1+b2)*sftx**ix *expterm*(-ep4*iz*sftz**iz+epep4*sftz**(iz+2)+ctmp-ep2*sftz**iz) - GTFdyyz=(c1+c2)*sftx**ix *expterm*(-ep4*iy*sfty**iy+epep4*sfty**(iy+2)+btmp-ep2*sfty**iy) !=dzyy - !dxxx,dyyy,dzzz - aatmp1=0D0 - bbtmp1=0D0 - cctmp1=0D0 - if (ix>=1) aatmp1=ep2*ix*sftx**(ix-1) - if (iy>=1) bbtmp1=ep2*iy*sfty**(iy-1) - if (iz>=1) cctmp1=ep2*iz*sftz**(iz-1) - aatmp2=0D0 - bbtmp2=0D0 - cctmp2=0D0 - if (ix>=2) aatmp2=ep2*ix*(ix-1)*sftx**(ix-1) - if (iy>=2) bbtmp2=ep2*iy*(iy-1)*sfty**(iy-1) - if (iz>=2) cctmp2=ep2*iz*(iz-1)*sftz**(iz-1) - aatmp3=0D0 - bbtmp3=0D0 - cctmp3=0D0 - if (ix>=3) aatmp3=ix*(ix-1)*(ix-2)*sftx**(ix-3) - if (iy>=3) bbtmp3=iy*(iy-1)*(iy-2)*sfty**(iy-3) - if (iz>=3) cctmp3=iz*(iz-1)*(iz-2)*sftz**(iz-3) - GTFdxxx=sfty**iy*sftz**iz*expterm*( (-2*ix+ep2*sftx2-1)*(-epep4*sftx**(ix+1) + aatmp1) - aatmp2 + epep8*sftx**(ix+1) + aatmp3 ) - GTFdyyy=sftx**ix*sftz**iz*expterm*( (-2*iy+ep2*sfty2-1)*(-epep4*sfty**(iy+1) + bbtmp1) - bbtmp2 + epep8*sfty**(iy+1) + bbtmp3 ) - GTFdzzz=sfty**iy*sftx**ix*expterm*( (-2*iz+ep2*sftz2-1)*(-epep4*sftz**(iz+1) + cctmp1) - cctmp2 + epep8*sftz**(iz+1) + cctmp3 ) - - !Calculate 3-order derivative tensor for orbital wavefunction - do imo=istart,iend - tens3(1,1,1,imo)=tens3(1,1,1,imo)+CO_uniq(imo,j)*GTFdxxx !dxxx - tens3(2,2,2,imo)=tens3(2,2,2,imo)+CO_uniq(imo,j)*GTFdyyy !dyyy - tens3(3,3,3,imo)=tens3(3,3,3,imo)+CO_uniq(imo,j)*GTFdzzz !dzzz - tens3(1,2,2,imo)=tens3(1,2,2,imo)+CO_uniq(imo,j)*GTFdxyy !dxyy - tens3(1,1,2,imo)=tens3(1,1,2,imo)+CO_uniq(imo,j)*GTFdxxy !dxxy - tens3(1,1,3,imo)=tens3(1,1,3,imo)+CO_uniq(imo,j)*GTFdxxz !dxxz - tens3(1,3,3,imo)=tens3(1,3,3,imo)+CO_uniq(imo,j)*GTFdxzz !dxzz - tens3(2,3,3,imo)=tens3(2,3,3,imo)+CO_uniq(imo,j)*GTFdyzz !dyzz - tens3(2,2,3,imo)=tens3(2,2,3,imo)+CO_uniq(imo,j)*GTFdyyz !dyyz - tens3(1,2,3,imo)=tens3(1,2,3,imo)+CO_uniq(imo,j)*GTFdxyz !dxyz - end do - tens3(1,2,1,:)=tens3(1,1,2,:) !dxyx=dxxy - tens3(1,3,1,:)=tens3(1,1,3,:) !dxzx=dxxz - tens3(1,3,2,:)=tens3(1,2,3,:) !dxzy=dxyz - tens3(2,1,1,:)=tens3(1,1,2,:) !dyxx=dxxy - tens3(2,1,2,:)=tens3(1,2,2,:) !dyxy=dxyy - tens3(2,1,3,:)=tens3(1,2,3,:) !dyxz=dxyz - tens3(2,2,1,:)=tens3(1,2,2,:) !dyyx=dxyy - tens3(2,3,1,:)=tens3(1,2,3,:) !dyzx=dxyz - tens3(2,3,2,:)=tens3(2,2,3,:) !dyzy=dyyz - tens3(3,1,1,:)=tens3(1,1,3,:) !dzxx=dxxz - tens3(3,1,2,:)=tens3(1,2,3,:) !dzxy=dxyz - tens3(3,1,3,:)=tens3(1,3,3,:) !dzxz=dxzz - tens3(3,2,1,:)=tens3(1,2,3,:) !dzyx=dxyz - tens3(3,2,2,:)=tens3(2,2,3,:) !dzyy=dyyz - tens3(3,2,3,:)=tens3(2,3,3,:) !dzyz=dyzz - tens3(3,3,1,:)=tens3(1,3,3,:) !dzzx=dxzz - tens3(3,3,2,:)=tens3(2,3,3,:) !dzzy=dyzz - end if !end runtype>=5 - - end if !end runtype>=3 - end if !end runtype>=2 - end do - end if - - end do - end do -end do -end subroutine - - - - -!!--------- The same as subroutine orbderv, but calculate for promolecular wavefunction (CO_pmol ...), up to Hessian -!! istart and iend is the range of the orbitals will be calculated, to calculate all orbitals, use 1,nmo_pmol -!! runtype=1: value =2: value+dx/y/z =3: value+dxx/yy/zz(diagonal of hess) =4: value+dx/y/z+Hessian -subroutine orbderv_pmol(runtype,istart,iend,x,y,z,wfnval,grad,hess) -real*8 x,y,z,wfnval(nmo_pmol) -real*8,optional :: grad(3,nmo_pmol),hess(3,3,nmo_pmol) -integer runtype,istart,iend - -wfnval=0D0 -if (present(grad)) grad=0D0 -if (present(hess)) hess=0D0 -lastcen=-1 !Arbitrary value - -!If the center/exp of current GTF is the same as previous, then we do not need to recalculate them -do j=1,nprims - ix=type2ix(b(j)%type) - iy=type2iy(b(j)%type) - iz=type2iz(b(j)%type) - ep=b(j)%exp - - if (b(j)%center/=lastcen) then - sftx=x-a(b(j)%center)%x - sfty=y-a(b(j)%center)%y - sftz=z-a(b(j)%center)%z - sftx2=sftx*sftx - sfty2=sfty*sfty - sftz2=sftz*sftz - rr=sftx2+sfty2+sftz2 - end if - if (expcutoff>0.or.-ep*rr>expcutoff) then - expterm=exp(-ep*rr) - else - expterm=0D0 - end if - lastcen=b(j)%center - if (expterm==0D0) cycle - - !Calculate value for current GTF - GTFval=sftx**ix *sfty**iy *sftz**iz *expterm - !Calculate orbital wavefunction value - do imo=istart,iend - wfnval(imo)=wfnval(imo)+CO_pmol(imo,j)*GTFval - end do - - if (runtype>=2) then - !Calculate 1-order derivative for current GTF - tx=0D0 - ty=0D0 - tz=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_pmol(imo,j)*GTFdx - grad(2,imo)=grad(2,imo)+CO_pmol(imo,j)*GTFdy - grad(3,imo)=grad(3,imo)+CO_pmol(imo,j)*GTFdz - end do - - if (runtype>=3) then - !Calculate 2-order derivative for current GTF - txx=0D0 - tyy=0D0 - tzz=0D0 - if (ix>=2) txx=ix*(ix-1)*sftx**(ix-2) - if (iy>=2) tyy=iy*(iy-1)*sfty**(iy-2) - if (iz>=2) tzz=iz*(iz-1)*sftz**(iz-2) - GTFdxx=sfty**iy *sftz**iz *expterm*( txx + 2*ep*sftx**ix*(-2*ix+2*ep*sftx2-1) ) - GTFdyy=sftx**ix *sftz**iz *expterm*( tyy + 2*ep*sfty**iy*(-2*iy+2*ep*sfty2-1) ) - GTFdzz=sftx**ix *sfty**iy *expterm*( tzz + 2*ep*sftz**iz*(-2*iz+2*ep*sftz2-1) ) - ttx=tx-2*ep*sftx**(ix+1) - tty=ty-2*ep*sfty**(iy+1) - ttz=tz-2*ep*sftz**(iz+1) - GTFdxy=sftz**iz *expterm*ttx*tty - GTFdyz=sftx**ix *expterm*tty*ttz - GTFdxz=sfty**iy *expterm*ttx*ttz - !Calculate diagonal Hessian elements for orbitals - do imo=istart,iend - hess(1,1,imo)=hess(1,1,imo)+CO_pmol(imo,j)*GTFdxx !dxx - hess(2,2,imo)=hess(2,2,imo)+CO_pmol(imo,j)*GTFdyy !dyy - hess(3,3,imo)=hess(3,3,imo)+CO_pmol(imo,j)*GTFdzz !dzz - end do - if (runtype>=4) then !Also process nondiagonal elements - do imo=istart,iend - hess(1,2,imo)=hess(1,2,imo)+CO_pmol(imo,j)*GTFdxy !dxy - hess(2,3,imo)=hess(2,3,imo)+CO_pmol(imo,j)*GTFdyz !dyz - hess(1,3,imo)=hess(1,3,imo)+CO_pmol(imo,j)*GTFdxz !dxz - end do - hess(2,1,:)=hess(1,2,:) - hess(3,2,:)=hess(2,3,:) - hess(3,1,:)=hess(1,3,:) - end if - end if !end runtype>=3 - end if !end runtype>=2 -end do + end subroutine orbderv_eval end subroutine @@ -2082,7 +1408,7 @@ end subroutine subroutine calchessmat_dens_promol(itype,x,y,z,elerho,elegrad,elehess) real*8 x,y,z,elerho,elegrad(3),elehess(3,3),wfnval(nmo_pmol),wfnderv(3,nmo_pmol),wfnhess(3,3,nmo_pmol) integer itype -call orbderv_pmol(4,1,nmo_pmol,x,y,z,wfnval,wfnderv,wfnhess) +call orbderv(4,1,nmo_pmol,x,y,z,wfnval,wfnderv,wfnhess,promol=.true.) elerho=sum( MOocc_pmol(1:nmo_pmol)*wfnval(1:nmo_pmol)*wfnval(1:nmo_pmol) ) do itmp=1,3 elegrad(itmp)=2*sum( MOocc_pmol(1:nmo_pmol)*wfnval(1:nmo_pmol)*wfnderv(itmp,1:nmo_pmol) ) @@ -5003,7 +4329,7 @@ real*8 x,y,z,wfnval_pmol(nmo_pmol) rho=fdens(x,y,z) -call orbderv_pmol(1,1,nmo_pmol,x,y,z,wfnval_pmol) +call orbderv(1,1,nmo_pmol,x,y,z,wfnval_pmol,promol=.true.) rho_0=sum( MOocc_pmol(1:nmo_pmol)*wfnval_pmol(1:nmo_pmol)**2 ) relShannon=rho*log(rho/rho_0) @@ -6782,4 +6108,4 @@ else if (ifunc==21) then end do write(*,*) end if -end subroutine \ No newline at end of file +end subroutine -- 2.34.1