************************************************************************
C                            MGAUSS_GRSOL.F
************************************************************************
C Calculate the mean (rbar) and the covariance (a) of a vector r, as 
C described in CP, sec 3.3.  The matrix a's inverse and determinant are 
C calculated.  
C Periodic BC exist only in the (xy) plane.  
C Uses gasdev,ratint,locate.
************************************************************************
      subroutine mgauss(x,y,z,xbar0,ybar0,zbar0,r,
     x     dutfix,ddutfix,iatom,idum,natp,sig2,nr,ipc,
     x     islice,nslm,klev,move,t,det)
      include 'pimc_sqr3sc.dcl'
 
      dimension x(2*m,nat),y(2*m,nat),z(2*m,nat)
      dimension xbar0(4),ybar0(4),zbar0(4),xbar(4),ybar(4),zbar(4)
      dimension r(nr1)
      dimension dutfix(lev,nr1),ddutfix(lev,nr1)
 
      dimension sumkk(ndim,ndim),a(4*ndim,4*ndim),sinv(4*ndim,4*ndim)
      dimension ainv(4*ndim,4*ndim)
      dimension p(4*ndim),sumchi(4*ndim),delr(4*ndim)
      dimension r0(ndim)
      dimension ydum1(nait),ydum2(nait)
      
      integer iatom(nperm,4)

      told=0.0
      tnew=0.0
      alpha=0.0
      r0min=2.
      rmax=8.  !(nr)
      idiagnostic=0
      iargm=islice+nslm
      rstep=r(2)-r(1)
 
      ducut=0.0
      d2ucut=0.0
      rlow=(rmax-r(1))/rstep
      ilow=int(rlow)+1
c      k=min(max(ilow-(nait-1)/2,1),nr+1-nait)
      
      xi=rlow-ilow+1.                    ! 3-point interpolation
      duk=dutfix(klev,ilow)
      duk1=dutfix(klev,ilow+1)
      duk2=dutfix(klev,ilow+2)
      t1=duk+(duk1-duk)*xi
      t2=duk1+(duk2-duk1)*(xi-1.)
      ducut=t1+(t2-t1)*xi*.5

      xi=rlow-ilow+1.                    ! 3-point interpolation
      d2uk=ddutfix(klev,ilow)
      d2uk1=ddutfix(klev,ilow+1)
      d2uk2=ddutfix(klev,ilow+2)
      t1=d2uk+(d2uk1-d2uk)*xi
      t2=d2uk1+(d2uk2-d2uk1)*(xi-1.)
      d2ucut=t1+(t2-t1)*xi*.5
      
      do i=1,natp*ndim
         p(i)=0.0
         sumchi(i)=0.0
         delr(i)=0.0
         do j=1,natp*ndim
            a(i,j)=0.0
            ainv(i,j)=0.0
            sinv(i,j)=0.0
         end do
      end do
      DO i=1,natp
         ii=iatom(ipc,i)
         sumx=0.0
         sumy=0.0
         sumz=0.0
         DO k1=1,ndim         
            DO k2=1,ndim
               sumkk(k1,k2)=0.0
            END DO
         END DO

         DO j=1,nat               ! Mobile helium contribution
            iflag=0
            isave=0
            IF(j.eq.ii) GOTO 5
            
            call dstpbc(xbar0(i),x(iargm,j),xdim,t1,sgnx)
            call dstpbc(ybar0(i),y(iargm,j),ydim,t2,sgny)
            r0(1)=sgnx*t1
            r0(2)=sgny*t2
            r0(3)=zbar0(i)-z(iargm,j)
            r0ij=sqrt(r0(1)**2+r0(2)**2+r0(3)**2)
            IF(r0ij.gt.rmax)THEN
               du=0.0
               d2u=0.0
            ELSEIF(r0ij.lt.r0min)THEN
               rlow=(r0min-r(1))/rstep
               ilow=int(rlow)+1
               
               xi=rlow-ilow+1.
               duk=dutfix(klev,ilow)
               duk1=dutfix(klev,ilow+1)
               duk2=dutfix(klev,ilow+2)
               t1=duk+(duk1-duk)*xi
               t2=duk1+(duk2-duk1)*(xi-1.)
               du=t1+(t2-t1)*xi*.5
               
               xi=rlow-ilow+1.
               d2uk=ddutfix(klev,ilow)
               d2uk1=ddutfix(klev,ilow+1)
               d2uk2=ddutfix(klev,ilow+2)
               t1=d2uk+(d2uk1-d2uk)*xi
               t2=d2uk1+(d2uk2-d2uk1)*(xi-1.)
               d2u=t1+(t2-t1)*xi*.5
               
               du=du+du*(r0min-r0ij)
            ELSE
               rlow=(r0ij-r(1))/rstep
               ilow=int(rlow)+1
               
               xi=rlow-ilow+1.
               duk=dutfix(klev,ilow)
               duk1=dutfix(klev,ilow+1)
               duk2=dutfix(klev,ilow+2)
               t1=duk+(duk1-duk)*xi
               t2=duk1+(duk2-duk1)*(xi-1.)
               du=t1+(t2-t1)*xi*.5
               
               xi=rlow-ilow+1.
               d2uk=ddutfix(klev,ilow)
               d2uk1=ddutfix(klev,ilow+1)
               d2uk2=ddutfix(klev,ilow+2)
               t1=d2uk+(d2uk1-d2uk)*xi
               t2=d2uk1+(d2uk2-d2uk1)*(xi-1.)
               d2u=t1+(t2-t1)*xi*.5
               
            END IF

            sumx=sumx+du*r0(1)/r0ij
            sumy=sumy+du*r0(2)/r0ij
            sumz=sumz+du*r0(3)/r0ij
            DO k1=1,ndim
               r0ij2=r0ij*r0ij
               sumkk(k1,k1)=sumkk(k1,k1)
     x              +d2u*r0(k1)*r0(k1)/r0ij2
     x              +du*(r0ij*r0ij-r0(k1)*r0(k1))/r0ij**3
               DO k2=k1+1,ndim
                  r0ij2=r0ij*r0ij
                  sumkk(k1,k2)=sumkk(k1,k2)
     x                 +d2u*r0(k1)*r0(k2)/r0ij2
     x                 -du*r0(k1)*r0(k2)/r0ij/r0ij/r0ij
               END DO
            END DO   
 5       END DO

C Calculate elements of the mean

         xbar(i)=xbar0(i)-sig2*sumx
         t1=xbar(i)
         if(t1.gt.xdim) xbar(i)=t1-xdim      ! PBC
         if(t1.lt.0.) xbar(i)=t1+xdim

         ybar(i)=ybar0(i)-sig2*sumy          
         t1=ybar(i)
         if(t1.gt.ydim) ybar(i)=t1-ydim      ! PBC
         if(t1.lt.0.) ybar(i)=t1+ydim

         zbar(i)=zbar0(i)-sig2*sumz

C Calculate diagonal elements of the submatrices A(xx), A(xy), etc.
         
         DO k1=1,ndim
            lll=(k1-1)*natp
            a(i+lll,i+lll)=sig2-sig2*sig2*sumkk(k1,k1)
            DO k2=k1+1,ndim
               mmm=(k2-1)*natp
               a(i+lll,i+mmm)=-sig2*sig2*sumkk(k1,k2)
               a(i+mmm,i+lll)=a(i+lll,i+mmm)
            END DO
         END DO

      END DO
      
C Cholesky decompose A into S.  Find SINV, S's inverse.

      nmat=ndim*natp

 200  format(1x,12(f6.2,1x))
      
      cholflg=0.0
      asum=0.0
      det=1.
      DO i=1,nmat
         DO j=i,nmat
            asum=a(i,j)
c            if(move.eq.0)print*,i,j,a(i,j)
            DO k=i-1,1,-1
               asum=asum-a(i,k)*a(j,k)
            END DO
            IF (i.eq.j)then
               IF(asum.le.0.) THEN  ! Choldc failed.  
c                  print*,'choldc failed'
                  move=2
                  t=0.0
                  det=0.0
                  goto 750
               ENDIF
               p(i)=sqrt(asum)
               det=det*p(i)
            ELSE
               a(j,i)=asum/p(i)
               sinv(j,i)=a(j,i)
            END IF
         END DO
      END DO

 500  continue  
      
      DO i=1,nmat              ! Calc SINV, inverse of matrix S
         sinv(i,i)=1./p(i)
         DO j=i+1,nmat
            sums=0.0
            DO k1=i,j-1
               sums=sums-sinv(j,k1)*sinv(k1,i)
            END DO
            sinv(j,i)=sums/p(j)
         END DO
         sumchi(i)=0.0          ! Convenient zeroing of sumchi
      END DO

C Calculate proposed new particle positions. 
         
      IF(move.eq.1) THEN
         DO j=1,nmat
            chi=gasdev(idum)
            sumchi(j)=sumchi(j)+p(j)*chi
            DO i=j+1,nmat
               sumchi(i)=sumchi(i)+a(i,j)*chi
            END DO
         END DO
         
         DO i=1,natp
c            print*,x(iargm,ii),y(iargm,ii),z(iargm,ii)
            ii=iatom(ipc,i)
            x(iargm,ii)=sumchi(i)+xbar(i)
            y(iargm,ii)=sumchi(i+natp)+ybar(i)
            z(iargm,ii)=sumchi(i+2*natp)+zbar(i)
            t1=x(iargm,ii)
            t2=y(iargm,ii)
            t3=z(iargm,ii)
            IF(t1.gt.xdim) x(iargm,ii)=t1-xdim
            IF(t2.gt.ydim) y(iargm,ii)=t2-ydim
            If(t1.lt.0.) x(iargm,ii)=xdim+t1
            IF(t2.lt.0.) y(iargm,ii)=ydim+t2
c            print*,x(iargm,ii),y(iargm,ii),z(iargm,ii)
c            print*
         END DO
      END IF

C Calculate transition amplitude

      DO i=1,natp
         ii=iatom(ipc,i)
         call dstpbc(x(iargm,ii),xbar(i),xdim,t1,sgn)
         delr(i)=sgn*t1
         call dstpbc(y(iargm,ii),ybar(i),ydim,t1,sgn)
         delr(i+natp)=sgn*t1
         t1=z(iargm,ii)-zbar(i)
         delr(i+2*natp)=t1
      END DO
      DO i=1,natp*ndim
         DO j=i,natp*ndim
            sainv=0.0
            DO k=j,natp*ndim
               sainv=sainv+sinv(k,j)*sinv(k,i)
            END DO
            ainv(j,i)=sainv
            ainv(i,j)=ainv(j,i)
         END DO
      END DO
      sumfinal=0.0 
      DO i=1,natp*ndim
c         if(move.eq.0) print*,i,delr(i)
         DO j=1,natp*ndim
c            if(move.eq.0) print*,i,j,ainv(i,j)
            sumfinal=sumfinal+delr(i)*ainv(i,j)*delr(j)/2.
         END DO   
      END DO
      t=-sumfinal
 750  continue
 1000 RETURN
      END
************************************************************************
      subroutine locate(r,nr,r0ij,num,k)
c      include 'pimc.dcl'
      implicit real*4 (a-h,o-z)
      dimension r(nr)
      klo=1
      khi=nr
 1    if(khi-klo.gt.1) then
	 kdum=(khi+klo)/2
	 if(r(kdum).gt.r0ij) then
	    khi=kdum
         else
            klo=kdum
         endif
         goto 1
      endif
      k=min(max(klo-(num-1)/2,1),nr+1-num)
      return
      end
