******************************************************************************** * * * 2DHF version 1-2003 * * Copyright (C) 1996 Jacek Kobus, Leif Laaksonen, Dage Sundholm * * * * This software may be used and distributed according to the terms * * of the GNU General Public License, see README and COPYING. * * * ******************************************************************************** c ### inigauss ### c c Molecular orbitals are initialized as a linear combination of c gaussian orbitals. Parameters of the orbitals and the coeeficients c are provided by the GAUSSIAN program. c c Coulomb (HF) potentials are initialized as a linear combination of c -ez1/r1 and -ez2/r2. For the initialization of exchange potentials c see routine tfpot. subroutine inigauss (psi,pot,excp,f2,f4,wgt2,wk0) implicit integer*4 (i-n) implicit real*8 (a-h,o-z) include 'commons8.inc' c set the maximum value of an argument of the exp function parameter (exp_max=700.d0) parameter (thnorm2=1.d-2, thnorm3=1.d-2) dimension psi(*),pot(*),excp(*),f2(*),f4(*),wgt2(*),wk0(*) c evaluate overlap matrix over gaussian basis functions igauss=0 if (idbg(552).ne.0) igauss=1 if (igauss.ne.0) then call gauss_ovlap (psi,pot,excp,f2,f4,wgt2,wk0) stop endif c Initialization of molecular orbitals if (ini.eq.4.) then norbt=1 else norbt=norb endif write(*,1114) icen1=1 icen2=2 icen3=3 iorb=0 iorbnorm=1 900 continue c Gaussian normal output c In the case of homonuclear molecules Gaussian program reverses c the order of the centres (with respect to the order defined c by the input data (sigma g/u and nonsigma orbitals are effected). c That is why the other ordering is tried when the norm of a numerical c orbital is off from unity by thnorm2. if ((ini.eq.2.or.ini.eq.3).and.idbg(551).eq.0) then icen1=1 icen2=2 icen3=3 xnorm_prv=0.d0 elseif ((ini.eq.2.or.ini.eq.3).and.idbg(551).eq.1) then icen1=3 icen2=2 icen3=1 endif iorb=iorb+1 if (iorb.gt.norbt) goto 910 if (iorbnorm.eq.1) then if(ini.eq.3.and.mm(iorb).ne.0) then iprint=0 else iprint=1 endif if (ini.eq.2) then iprint=0 endif endif psi00=0.d0 psi01=0.d0 if (idbg(550).ne.0) then i1b(iorb)=1 ngrid= i1si(iorb) do i=1,ngrid psi(i)=0.d0 enddo endif ishift=i1b(iorb)-1 ngrid= i1si(iorb) igp=ishift do i=1,ngrid psi(ishift+i)=0.d0 enddo c loop over basis set functions ipnzero=0 nc2=0 c1sum=0.d0 do ipb=1,npbasis l1 =lprim(ipb) m1 =iabs(mprim(ipb)) if (m1.ne.mm(iorb)) goto 999 c Gaussian long (customized) output. c In most cases PX components are used to construct numerical c molecular orbital. Sometimes PY components are used instead to c get the norm of the orbital correct. If the norm is off by more c than thnorm3 than idbg(551) is set to 1 and the orbital is c initialized afresh. if (ini.eq.3.and.idbg(551).eq.0) then if (mprim(ipb).lt.0) goto 999 elseif (ini.eq.3.and.idbg(551).eq.1) then if (mprim(ipb).gt.0) goto 999 endif ic1 =icgau(ipb) d1 =primexp(ipb) c sometimes the gaussian long output contains unprintable exponents c which must be singled out if (ini.eq.3) then if (d1.gt.99999999.d0) goto 999 endif fn1 =fngau2(ipb) shn1=shngau(ipb) fnorm=fn1*shn1 c1 =primcoef(iorb,ipb) if (abs(c1).gt.0.d0) then ipnzero=ipnzero+1 c write(*,1111) iorb,ipb,ic1,l1,m1,d1,c1,fn1,shn1 c write(*,1112) iorb,ipb,ic1,l1,mprim(ipb),c1,d1 1111 format(5i4,f20.9,3d15.5) 1112 format(5i4,2f20.9) c ic1=1 basis function at centre Z1 c ic1=2 basis function at bond centre c ic1=3 basis function at centre Z2 c loop over grid points do imu=1,mxnmu inioff=(imu-1)*nni do in=1,nni igp=ishift+inioff+in c for each grid point, i.e. for (vmu(imu),vni(ini)) determine c its distance |_r1| from the nuclei Z_1 and Z_2 and cosine c of the polar angles costh1 and costh2 betp ween z axis c and the vectors _r1 and _r2 c rr=(r/2.d0)*sqrt(vxisq(imu)+vetasq(in)-1.d0) z=(r/2.d0)*vxi(imu)*veta(in) if (ic1.eq.icen1) then c gaussians centred on Z1 r1=(r/2.d0)*(vxi(imu)+veta(in)) if (r1.lt.precis) then costh1=0.d0 else costh1=(z+r/2.d0)/r1 endif elseif (ic1.eq.icen3) then c gaussians centred on Z2 r1=(r/2.d0)*(vxi(imu)-veta(in)) if (r1.lt.precis) then costh1=0.d0 else costh1=(z-r/2.d0)/r1 endif elseif (ic1.eq.icen2) then c gaussians centred on the bond centre r1=sqrt(vxisq(imu)+vetasq(in)-1.d0)*r/2.d0 if (r1.lt.precis) then costh1=0.d0 else costh1= (z)/r1 endif endif c calculate the gaussian function centered on one of the nuclei d1r=d1*r1*r1 if (d1r.gt.exp_max) then expt=0.d0 else expt=1.d0/exp(d1r) endif if (r1.lt.precis.and.l1.eq.0) then psi1=fnorm*expt*plegendg(l1,m1,costh1) elseif (r1.lt.precis.and.l1.gt.0) then psi1=0.d0 else psi1=fnorm*r1**l1*expt*plegendg(l1,m1,costh1) endif 99 continue psi(igp)=psi(igp)+c1*psi1 enddo enddo c c1.ne.0 endif c 999 continue enddo call norm94 (iorb,psi,f4,wgt2,wk0,xnorm) write (*,1115) iorn(iorb),bond(iorb),gut(iorb),xnorm if (iprint.eq.1) write (*,1115) iorn(iorb),bond(iorb), & gut(iorb),xnorm if (iorbnorm.eq.1) then c if(ini.eq.2) then if(ini.eq.2.and. & abs(xnorm-1.d0).gt.thnorm2) then c Try reversed ordering of atomic centres. This might lead to a better orbital idbg(551)=1 iorb=iorb-1 iorbnorm=2 xnorm_prev=xnorm endif c if(ini.eq.3.and.mm(iorb).ne.0) then if(ini.eq.3.and.mm(iorb).ne.0.and. & abs(xnorm-1.d0).gt.thnorm3) then idbg(551)=1 iorb=iorb-1 iorbnorm=2 xnorm_prev=xnorm endif elseif (iorbnorm.eq.2) then if (abs(xnorm-1.d0).le.abs(xnorm_prev-1.d0)) then iorbnorm=3 idbg(551)=0 iprint=1 else iorbnorm=3 xnorm=xnorm_prev iorb=iorb-1 idbg(551)=0 iprint=1 endif endif goto 900 910 continue write(*,*) if (idbg(550).ne.0) stop 'inigauss' c initialize Coulomb and exchange potentials call tfpot(psi,pot,excp,f2,f4,wk0) 1114 format(/1x,' orbital norm ') 1115 format(1x,i3,1x,a8,a1,d20.12) return end