********************************************************************************
*                                                                              *
*  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 ### prepvar ###
c
c  Initialize arrays of variable length
c
      subroutine prepvar (borb,bpot,d,e,f0,f1,f2,f3,f4,g,
     &		          wjac1,wjac2,wgt1,wgt2)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension borb(nni,mxnmu),bpot(nni,mxnmu),
     &          d(nni,mxnmu),e(nni,mxnmu),
     &		f0(nni,mxnmu),f1(nni,mxnmu),f2(nni,mxnmu),
     &          f3(nni,mxnmu),f4(nni,mxnmu),g(nni,mxnmu),
     &          wjac1(nni,mxnmu),wjac2(nni,mxnmu),
     &          wgt1(mxsize),wgt2(mxsize)

c     iaddr4(1)=borb
c     iaddr4(2)=bpot
c     iaddr4(3)=d
c     iaddr4(4)=e
c     iaddr4(5)=f0
c     iaddr4(6)=f1
c     iaddr4(7)=f2
c     iaddr4(8)=f3
c     iaddr4(9)=f4
c     iaddr4(10)=g
c     iaddr4(11)=wjac1
c     iaddr4(12)=wjac2
c     iaddr4(13)=wgt1
c     iaddr4(14)=wgt2

c     wjac1 - jacobian for rayleigh quotient
c     wjac2 - jacobian for two-electron integrals

c     initialize arrays needed for the poisson equations

c     f4 - array of factors multiplying off-diagonal lagrange multipliers
c	   in the hf equation
c     off-diagonal l.m. must be multiplied by the same factor as
c     diagonal one i.e. f1 instead of f4

      if (ifefield.ne.0) then
         izz1=nint(z1)
         izz2=nint(z2)
         atw1=atweight(izz1)
         atw2=atweight(izz2)
         zcm=(-atw1+atw2)/(atw1+atw2)*r/2.d0
c         write(iout6,1000) element(izz1), atweight(izz1), 
c     &        element(izz2), atweight(izz2),zcm
c1000    format(5x,'atomic weight: ',a2,' --',f14.8,/,
c     &        5x,'               ',a2,' --',f14.8,/,
c     &        5x,'    z_cm (bohr)','     ',f14.8)
      endif
      
      wj1 =-r*pii/2.d0
      wj2 = r*r*pii/2.d0 
      xrr=r/2.d0
      ib=ibmu(1)
      do ig=1,ngrids
         ie=iemu(ig)
         diag(ig)=-14350.d0/(5040.d0)*
     &        ( 1.d0/(hmu(ig)*hmu(ig))+1.d0/(hni*hni) )
         do i=ib,ie
            w1=vxi1(i)
            w2=vxi(i)*vxi(i)
            w9=xrr*vxi(i)

c           ssf is included in differentiation array dmu1
c           w4=vxi2(i)/ssf

            w4=vxi2(i)
            if(i.eq.1) then
               w5=0.d0
               w6=0.d0
            else
               w5= vxi(i)/vxi1(i)-2.0*vxi1(i)/vxi(i)
               w6= 1.d0/(vxi1(i)*vxi1(i))
            endif
            
            w8=-r*r*r*pii*vxi(i)/2.d0
            
            do j=1,nni
               w10=veta1(j)
               borb(j,i)=w4
               bpot(j,i)=w5
               d   (j,i)=veta2(j)
               if (w10.lt.precis) then
                  e(j,i)= 0.d0
               else
                  e(j,i)=-w6 - 1.d0/(w10*w10)
               endif

               if (ifermi.eq.0) then
                  f0(j,i)=r*(vxi(i)*(z1+z2)+veta(j)*(z2-z1))
               elseif (ifermi.eq.1) then
                  z1t=zz1(i,j)
                  z2t=zz2(i,j)
                  f0(j,i)=r*(vxi(i)*(z1t+z2t)+veta(j)*(z2t-z1t))
c                 write(*,*) i,j,z1t,z2t,f0(j,i)
               elseif (ifermi.eq.2) then
                  z1t=zz1g(i,j)
                  z2t=zz2g(i,j)
                  f0(j,i)=r*(vxi(i)*(z1t+z2t)+veta(j)*(z2t-z1t))
c                 write(*,*) i,j,z1t,z2t,f0(j,i)
               endif

               w7=w2-veta(j)*veta(j)
               f1(j,i)= r*r*w7/2.d0
               
               if (ifefield.eq.1) then

c                 when external electric field is present (z-zcm)*E 
c                 has to be added (times -f1)

                  z=(r/2.d0)*vxi(i)*veta(j)
                  if (abs(z).lt.zcutoff) then
                     f0(j,i)=f0(j,i)-ffield*(z-zcm)*f1(j,i)
                  else
                     f0(j,i)=f0(j,i)-ffield*(z-zcm)*
     &                    exp(zcutoff-abs(z))*f1(j,i)
                  endif 
                  
c                 stop 'prepvar: ifefield'
               elseif (ifefield.eq.2) then

c                 if electric field gradient is nonzero the interaction 
c                 energy is -(Q_z/4)*(dE/dz) 

c                 if electric field gradient is nonzero the interaction 
c                 energy is  -R^2/3 P(2,0)*(dE/dz) (Greiner, Electrodynamics, p.91)

                  imu=i
                  in=j
                  z=(r/2.d0)*vxi(imu)*veta(in)-zcm
                  xxplusyy=((r/2.d0)*vxi1(imu)*veta1(in))**2
                  rr=sqrt(xxplusyy+z*z)
                  if (abs(rr).lt.precis) then
                     costh=0.d0
                  else
                     costh=z/rr
                  endif
                  
                  wktmp=rr**2*plegendg(2,0,costh)
                  z=(r/2.d0)*vxi(i)*veta(j)
                  if (abs(z).lt.zcutoff) then
                     f0(j,i)=f0(j,i)-fgrad*wktmp*f1(j,i)/3.d0
                  else
                     f0(j,i)=f0(j,i)-fgrad*wktmp*f1(j,i)/3.d0
     &                    *exp(zcutoff-abs(z))
                  endif
               endif
               
               if(i.eq.1) then
                  f2(j,i)=0.d0
               else
                  f2(j,i)=-r*w7/vxi(i)
               endif
               f3(j,i)=-2.d0/w2
               g (j,i)= w8*w7
               f4(j,i)= w9
               
c              Integration needs hmu*hni but the integration routine works
c              for hni*hni. That is why ssf must be inserted here.

c              wjac1(j,i)=wj1*w1*w5/ssf(ig)
c              wjac2(j,i)=wj2*w1*w5*(w3+w6)/(ssf(ig)*w0)

               wjac1(j,i)=wj1*w1*w10
               wjac2(j,i)=wj2*w1*w10*(w1*w1+w10*w10)/vxi(i)
            enddo
         enddo
         ib=ie+1
c        enddo over ig
      enddo

c     initialize vector of integration weights
c     (7-point formula	abramovic 25.4.16.)

      wmup=0.d0
      do ig=1,ngrids
         ib=ibmu(ig)
         ie=iemu(ig)
         w1 =hmu(ig)/140.d0
         wmu(ib)=wmup+w1*41.d0
         do i=ib+1,ie,6
            wmu(i  )=w1*216.d0
            wmu(i+1)=w1* 27.d0
            wmu(i+2)=w1*272.d0
            wmu(i+3)=w1* 27.d0
            wmu(i+4)=w1*216.d0
            if (i+5.eq.ie) then
               wmu(i+5)=w1*41.d0
            else
               wmu(i+5)=2.d0*w1*41.d0
            endif
         enddo
         wmup=wmu(ie)
      enddo

      w1    =hni/140.d0
      wni(1)=w1*41.d0
      do i=2,nni,6
         wni(i  )=w1*216.d0
         wni(i+1)=w1* 27.d0
         wni(i+2)=w1*272.d0
         wni(i+3)=w1* 27.d0
         wni(i+4)=w1*216.d0
         if (i+5.eq.nni) then
            wni(i+5)=w1*41.d0
         else
            wni(i+5)=2.d0*w1*41.d0
         endif
c        enddo over i
      enddo

c     multiply jacobians for one and two-electron integrals by elements
c     of wmu array

      do i=1,mxnmu
         w1=wmu(i)
         do j=1,nni
            wjac1(j,i)=w1*wjac1(j,i)
            wjac2(j,i)=w1*wjac2(j,i)
         enddo
      enddo

c     prepare superarrays with the integration weights over ni variable
c     and with jacobians included

      do i=1,mxnmu
         ii=(i-1)*nni
         do j=1,nni
            wgt1(ii+j)=wni(j)*wjac1(j,i)
            wgt2(ii+j)=wni(j)*wjac2(j,i)
         enddo
      enddo

c     initialize extrapolation coefficients for even functions
c     bondary values for odd functions must be all zero

      exeven(1)= 210.d0/126.d0
      exeven(2)=-120.d0/126.d0
      exeven(3)=  45.d0/126.d0
      exeven(4)= -10.d0/126.d0
      exeven(5)=   1.d0/126.d0

      return
      end