******************************************************************************** * * * 2DHF version 9-2002 * * 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 ### prepfix ### c c Initializes arrays within common blocks and prepares grids. c subroutine prepfix implicit real*8 (a-h,o-z) character*8 sigma,pi,delta,phi character*8 spac,ger,uger include 'commons8.inc' data sigma/'sigma'/,pi/'pi'/,delta/'delta'/,phi/'phi'/ data spac/' '/,ger/'g'/,uger/'u'/ c c homonuclear case is treated as a heteronuclear one; if ihomon=2, i.e c if HOMO label is used, the g/u symmetry is imposed explicitely c ihomon=0 -- heteronuclear case c ihomon=1 -- homonuclear case |z1-z2|<1.d-6 c ihomon=2 -- homonuclear case with forced symmetry c ihomon=2 is equivalent to idbg(90).ne.0 c since in this version 7 point integration formula is used the number c of mesh points in the variable ni and the total number of points c in mu must be equal to 6*n+1 c Another restriction has to be imposed on the number of points in the c mu variable if the mcsor scheme is to be used: nmu must be c of the form kN+k+1, where k=iorder/2+1 for each grid (in this program c iorder, i.e. the order of descritization is 8) c checking and possibly adjusting the grid size so that nmu is c of the form kN+k+1, where k=iorder/2+1 c check first if ngrids does not exceed maxgrids if (ngrids.gt.3) then write(iout6,*) '===== prepfix =====' write(iout6,*) 'number of grids can not exceed',maxgrids stop 'prepfix' endif k1=5 k2=6 do ig=1,ngrids 20 continue ichk1=0 ichk2=0 11 nord=(nmu(ig)-1)/k2 if (k2*nord.ne.nmu(ig)-1) then nmu(ig)=nmu(ig)-1 ichk2=ichk2+1 goto 11 endif if (ichk1.ne.0.or.ichk2.ne.0) goto 20 enddo 12 nord=(nni-1)/k2 if (k2*nord.ne.nni-1) then nni=nni-1 goto 12 endif c mxnmu -- maximum no. of grid points in mu variable c mxmnu should be of the form 6N+1 c When counting the maximum number of points in mu variable do not forget c that the first point of the next grid is the last on of the previous grid. c Arrays ibmu and iemu contain the first and last addresses of c particular grid within the combined grid ibmu(1)=1 iemu(1)=nmu(1) mxnmu=nmu(1) ioffs(1)=0 do ig=2,ngrids mxnmu=mxnmu+nmu(ig)-1 ibmu(ig)=iemu(ig-1) iemu(ig)=mxnmu ioffs(ig)=nni*(iemu(ig-1)-1) enddo c determine the size of each grid do ig=1,ngrids ngsize(ig)=nni*nmu(ig) enddo c total no of grid points mxsize=nni*mxnmu if ( norb.gt.60 ) then write(iout6,*) '===== prepfix =====' write(*,*) 'number of orbitals can not exceed 60' stop endif nel=0 do i=1,norb nn(i)=mgx(1,i) ll(i)=mgx(3,i) mm(i)=mgx(3,i) iocc(i)=occ(i)+.000001d0 nel=nel+iocc(i) enddo c hni - step in ni variable c hmu - step in mu variable c determine step size in ni variable hni=pii/dble(nni-1) c determine step size in mu variable for each grid c when ngrids=1 hmu can be calculated from the practical infinity c variable if its value is provided in the input; if (ngrids.eq.1) then xmi0=2.d0*rgrid(1)/r xmi0=log(xmi0+sqrt(xmi0*xmi0-1.d0)) hmu(1)=xmi0/dble(mxnmu-1) rgrid(1)=hmu(1) else do ig=1,ngrids if(ig.eq.1) then hmu(ig)=rgrid(ig) else hmu(ig)=hmu(1)*rgrid(ig) endif enddo endif c initialize mu and xi arryas vmu(1) =0.d0 vxi(1) =cosh(vmu(1)) vxisq(1)=vxi(1)*vxi(1) vxi1(1)=sqrt(vxi(1)*vxi(1)-1.d0) vxi2(1)=0.0d0 ib=2 do ig=1,ngrids ie=iemu(ig) do i=ib,ie vmu(i)=vmu(i-1)+hmu(ig) vxi(i)=cosh(vmu(i)) vxisq(i)=vxi(i)*vxi(i) vxi1(i)=sqrt(vxi(i)*vxi(i)-1.d0) vxi2(i)=vxi(i)/vxi1(i) c enddo over i enddo ib=ie+1 c enddo over ig enddo c initialize ni and eta arrays do i=1,nni vni(i)=dble((i-1))*hni veta(i)=cos(vni(i)) vetasq(i)=veta(i)*veta(i) c sqrt(1-veta*veta) and veta/sqrt(1-veta*veta) c veta1(i)=sqrt(1.d0-veta(i)*veta(i)) if (veta1(i).lt.precis) then veta2(i)=0.0 else veta2(i)=veta(i)/veta1(i) endif enddo c determine rinf if (idbg(155).ne.0) then write(iout6,*) write(iout6,*) 'subgrid nni nmu hni hmu rb ', & ' heta hxib hxie' endif rinf=r*vxi(iemu(ngrids))/2.d0 heta=abs(veta(nni)-veta(nni-1)) ib=1 do ig=1,ngrids rinfig=r*vxi(iemu(ig))/2.d0 ie=iemu(ig) hxibeg=vxi(ib+1)-vxi(ib) hxiend=vxi(ie)-vxi(ie-1) if (idbg(155).ne.0) write(iout6,1000) ig,nni,nmu(ig), & hni,hmu(ig),rinfig,heta,hxibeg,hxiend ib=ie enddo 01000 format(4x,i2,3x,i4,2x,i4,2f9.5,f8.3,3f9.5) c ige(i)=1 ungerade c ige(i)=2 gerade or heteronuclear case c if break is on if (ibreak.eq.1) then do i=1,norb gut(i)=spac enddo endif c give a number within a symmetry isu=0 isi=0 ipu=0 ipi=0 idu=0 idi=0 ipsu=0 ips=0 do ibo=1,norb iba=norb+1-ibo if (bond(iba).eq.sigma) then if (gut(iba).eq.ger) then isi=isi+1 iorn(iba)=isi endif if (gut(iba).eq.spac) then isi=isi+1 iorn(iba)=isi endif if (gut(iba).eq.uger) then isu=isu+1 iorn(iba)=isu endif goto 100 endif if (bond(iba).eq.pi) then if (gut(iba).eq.ger) then ipi=ipi+1 iorn(iba)=ipi endif if (gut(iba).eq.spac) then ipi=ipi+1 iorn(iba)=ipi endif if (gut(iba).eq.uger) then ipu=ipu+1 iorn(iba)=ipu endif goto 100 endif if (bond(iba).eq.delta) then if (gut(iba).eq.ger) then idi=idi+1 iorn(iba)=idi endif if (gut(iba).eq.spac) then idi=idi+1 iorn(iba)=idi endif if (gut(iba).eq.uger) then idu=idu+1 iorn(iba)=idu endif goto 100 endif if (bond(iba).eq.phi) then if (gut(iba).eq.ger) then ips=ips+1 iorn(iba)=ips endif if (gut(iba).eq.spac) then ips=ips+1 iorn(iba)=ips endif if (gut(iba).eq.uger) then ipsu=ipsu+1 iorn(iba)=ipsu endif endif 100 continue enddo c ihsym = 1 - symmetry g c ihsym =-1 - symmetry u do i=1,norb ihomo(i)=1 c if (gut(i).eq.'u') ihomo(i)=-1 if (gut(i).eq.'u'.and.orbsym(i).eq.sigma) ihomo(i)=-1 if (gut(i).eq.'u'.and.orbsym(i).eq.delta) ihomo(i)=-1 if (gut(i).eq.'g'.and.orbsym(i).eq.pi) ihomo(i)=-1 if (gut(i).eq.'g'.and.orbsym(i).eq.phi) ihomo(i)=-1 enddo do ial=1,norb ige(ial)=2 if (gut(ial).eq.'u') ige(ial)=1 enddo c calculate the weights of the two-electron contributions to the total energy c expression for an open shell scf case call oshell c set optimal value of omega for potentials if ovfcoul(1)< 0 if (ovfcoul(1).lt.0.0) then ovfcoul(1)=setomega() ovfexch(1)=ovfcoul(1) endif c initialize array calp of coefficients of associated Legendre c functions needed by mulex and multi routines c data c/1.25d-01, 5.d-01, 5.59016994d-01, c 4.33012702d-01, 1.2247448710d0, 3.95284708d-01, c 1.369306394d0, 6.12372436d-01, 1.479019946d0, c 5.59016994d-01, 5.22912517d-01, 7.07106781d-01/ calp( 1)=0.125d0 calp( 2)=0.5d0 calp( 3)=sqrt( 5.d0/ 16.d0) calp( 4)=sqrt( 3.d0/ 16.d0) calp( 5)=sqrt( 3.d0/ 2.d0) calp( 6)=sqrt( 5.d0/ 32.d0) calp( 7)=sqrt(15.d0/ 8.d0) calp( 8)=sqrt( 3.d0/ 8.d0) calp( 9)=sqrt(35.d0/ 16.d0) calp(10)=sqrt( 5.d0/ 16.d0) calp(11)=sqrt(35.d0/128.d0) calp(12)=sqrt( 1.d0/ 2.d0) return end