module commvar
      integer,parameter :: n=2,n1=3*n*n*n*8,n2=n1+1,n3=8*n*n*n
      real*8 :: realn,mini
      real*8 :: pos(n3,3),dpos(n3,3),r0(3,3),o0,pi,dv(n1,n1),small
      integer :: neigh(n3,4),type(n3)
      logical :: flag(n3,n3)

!	Variables used
!
!	n = length of side of lattice in 8 atom blocks
!	n1,n2,n3 = sizes of arrays
!	pos = default positions of atoms
!	dpos = changes in positions of atoms
!     opos = positions loaded in
!	neigh(i,j) = j th neighbour of ith atom
!	r0(i,j)= perfect lattice length between atoms of type=i and type=j
!	o0= 109'
!	pi 
!	wrap = boundary conditions function 
!	type = ith atom is type i 1=Si 2=Ge 3=C
!

      contains

      subroutine values
      implicit none
      open(unit=101,file='mini2',status='old')
      read(101,*) mini
      close(101)
      r0(1,1)=dble(2.35156)
      r0(2,2)=dble(2.447)
      r0(3,3)=dble(1.545)
      r0(1,2)=(r0(1,1)+r0(2,2))/dble(2)
      r0(2,1)=r0(1,2)
      r0(2,3)=(r0(2,2)+r0(3,3))/dble(2)
      r0(3,2)=r0(2,3)
      r0(1,3)=dble(1.888)
      r0(3,1)=dble(1.888)
      realn=dble(4*n)*r0(1,1)/sqrt(dble(3))
      o0=acos(dble(-1)/dble(3))
      pi=acos(dble(-1))
      small=dble(.0001)
      end subroutine values

      end module commvar

!	MAIN PROGRAM

      program phd4
      use commvar
      use para
      implicit none
      call values
      call parameters
      call setup
      call leastsquares
      call writeout(2) 
      end program phd4

      subroutine lsqfun(iflag,m1,no,x,fvecc,iw,liw,w,lw)
      use commvar
      use para
      implicit none
      integer :: iflag,m1,no,liw,iw(liw),lw,i,j,a,b
      real*8 :: x(no),fvecc(m1),w(lw),sum
      sum=dble(0)
      bfr(1,2)=x(1)
      bfr(2,1)=x(1)
      call writeout(2)
      call relax
      call makedv
      call ek(sum)
      write(10,*) 'LVM=',sum,sum*dble(519)
      write(10,*)
      sum=sum-(dble(386)/dble(519))
      if (sum**2<mini) then
        mini=sum**2
        write(101,*) mini
        open(unit=102,file='mini2',status='replace')
        write(102,*) mini
        close(102)
        call writeout(3)
      end if
      fvecc(1)=sum
      return
      end subroutine
 
      subroutine lsqmon(m1,no,xc,fvecc,fjacc,ljc,s,igrade,niter,nf,&
                        iw,liw,w,lw)
      use commvar
      implicit none
      integer :: m1,no,ljc,igrade,niter,nf,liw,iw(liw),lw,i
      real*8 :: xc(no),fvecc(m1),fjacc(ljc,no),s(no),w(lw),sum
      write(10,*)
      write(10,*) 'niter=',niter,nf
      write(10,*)
      sum=dble(0)
      do i=1,m1
        sum=sum+fvecc(i)**2
      end do
      write(10,*) sum
      write(10,*)
      call writeout(2)
      return
      end subroutine lsqmon

      subroutine leastsquares
      use commvar
      use para
      implicit none
      integer :: m1,no,iprint,maxcal,lj,lv,niter,nf,liw,iw(10),lw,ifail
      real*8 :: eta,xtol,stepmx,x(1),fsumsq,fvec(1),fjac(1,1)
      real*8 :: s(1),v(1,1),w(10)
      integer :: i,j,k,l,m,a,b
      external lsqfun,lsqmon
      ifail=1
      m1=1
      no=1
      liw=10
      lj=m1
      lv=no
      lw=7+3*m1
      iprint=1
      maxcal=20000
      fsumsq=dble(0)
      eta=dble(0.0)
      xtol=dble(1.0e-9)
      stepmx=dble(.1)
      do i=1,3
        do j=1,3
          bfr(i,j)=sqrt(bfr(i,j))
        end do
      end do
      do i=1,n3
        type(i)=2
      end do
      type(1)=1      
      x(1)=bfr(1,2)
      call e04fcf(m1,no,lsqfun,lsqmon,iprint,maxcal,eta,xtol,stepmx,x,fsumsq,&
                  fvec,fjac,lj,s,v,lv,niter,nf,iw,liw,w,lw,ifail)
      bfr(1,2)=x(1)
      bfr(2,1)=x(1)
      call writeout(2)
      return
      end subroutine leastsquares

      subroutine setup
      use commvar
      implicit none
      integer :: count,set=n*n*n,alter,i,j,k,l,m
      real*8 :: match1,match2,match3,wrap
!	Set variables to zero
      do i=1,n3
            do j=1,3
                 pos(i,j)=dble(0)
                 dpos(i,j)=dble(0)
                 neigh(i,j)=0
            end do
            neigh(i,4)=0
      end do
      count=1
            do k=0,n-1
                 do j=0,n-1
                      do i=0,n-1
      do l=0,1
! fcc corners
      pos(count+l*4*set,1)=dble(i)+l*dble(0.25)
      pos(count+l*4*set,2)=dble(j)+l*dble(0.25)
      pos(count+l*4*set,3)=dble(k)+l*dble(0.25)
! centre of floor
      pos(count+(4*l+1)*set,1)=dble(i)+dble(0.5)+l*dble(0.25)
      pos(count+(4*l+1)*set,2)=dble(j)+dble(0.5)+l*dble(0.25)
      pos(count+(4*l+1)*set,3)=dble(k)+l*dble(0.25)
! front face
      pos(count+(4*l+2)*set,1)=dble(i)+dble(0.5)+l*dble(0.25)
      pos(count+(4*l+2)*set,2)=dble(j)+l*dble(0.25)
      pos(count+(4*l+2)*set,3)=dble(k)+dble(0.5)+l*dble(0.25)
! side face
      pos(count+(4*l+3)*set,1)=dble(i)+l*dble(0.25)
      pos(count+(4*l+3)*set,2)=dble(j)+dble(0.5)+l*dble(0.25)
      pos(count+(4*l+3)*set,3)=dble(k)+dble(0.5)+l*dble(0.25)
      end do
      count=count+1
                       end do
                 end do
           end do
      do i=1,n3/2
        match1=wrap(pos(i,1)+dble(0.25))
        match2=wrap(pos(i,2)+dble(0.25))
        match3=wrap(pos(i,3)+dble(0.25))
        call matched(1,i,match1,match2,match3)
        match1=wrap(pos(i,1)-dble(0.25))
        match2=wrap(pos(i,2)-dble(0.25))
        match3=wrap(pos(i,3)+dble(0.25))
        call matched(2,i,match1,match2,match3)
        match1=wrap(pos(i,1)-dble(0.25))
        match2=wrap(pos(i,2)+dble(0.25))
        match3=wrap(pos(i,3)-dble(0.25))
        call matched(3,i,match1,match2,match3)
        match1=wrap(pos(i,1)+dble(0.25))
        match2=wrap(pos(i,2)-dble(0.25))
        match3=wrap(pos(i,3)-dble(0.25))
        call matched(4,i,match1,match2,match3)
      end do
      do i=1,n3
        do j=1,3
           pos(i,j)=dble(4)*pos(i,j)*r0(1,1)/dsqrt(dble(3))
        end do
      end do
      do i=1,n3
       do j=1,n3
         if(i==j) flag(i,j)=.true.
       end do
      do j=1,4
      do k=1,4
       if (j/=k) then
       flag(i,neigh(neigh(i,j),k))=.true.
      do l=1,4
      if(l/=j.and.l/=k) then
            flag(i,neigh(neigh(neigh(i,k),j),l))=.true.
            flag(i,neigh(neigh(i,k),l))=.true.
                 do m=1,4
                 if(m/=j.and.m/=k.and.m/=l) then
              flag(i,neigh(neigh(neigh(i,k),l),m))=.true.
              flag(i,neigh(neigh(i,l),m))=.true.
                 end if
                 end do
      end if
      end do
      end if
      end do
      end do
      end do
      return
      end subroutine setup

      real*8 function wrap(a)
      use commvar
      implicit none
      real*8 ::a
      wrap=a
      if (a==dble(-0.25))wrap=dble(n)-dble(0.25)
      if (a==dble(n)) wrap=dble(0)
      end function wrap

      subroutine matched(a,b,c1,c2,c3)
      use commvar
      implicit none
      integer :: a,b,c0
      real*8 :: c1,c2,c3
      c0=n3/2
10    if (pos(c0,1)/=c1.or.pos(c0,2)/=c2.or.pos(c0,3)/=c3) then
        c0=c0+1
        goto 10
      end if
      neigh(b,a)=c0
      neigh(c0,a)=b
      flag(b,c0)=.true.
      flag(c0,b)=.true.
      return
      end subroutine matched


      subroutine relax
      use commvar
      implicit none
      integer,parameter :: liw=n2+2,lw=12*n2+n2*(n2-1)/2
      real*8 ::bl(n2),ul(n2),x(n2),f,w(lw)
      integer :: iw(liw),ifail,ibound
      integer :: i,j
      do i=1,n3
        do j=1,3
          x(3*(i-1)+j)=dpos(i,j)
        end do
      end do
      do i=1,n2
        bl(i)=dble(-4)
        ul(i)=dble(4)
      end do
      bl(n2)=realn-dble(4)
      x(n2)=realn
      ul(n2)=realn+dble(4)
      ifail=1
      ibound=0
      call e04jaf(n2,ibound,bl,ul,x,f,iw,liw,w,lw,ifail)
      write(10,*) 'Relax=',f
      realn=x(n2)
      write(10,*) 'realn=',realn
      do i=1,n3
        do j=1,3
          dpos(i,j)=x(3*(i-1)+j)
          if(abs(dpos(i,j))==dble(4)) write(10,*) 'Maxed out',i,j
        end do
      end do
      return
      end subroutine relax

      subroutine funct1(no,x,energy)
      use commvar
      implicit none
      integer ::no,i,j
      real*8 :: x(no),energy
      do i=1,n3
        do j=1,3
          dpos(i,j)=x(3*(i-1)+j)
        end do
      end do
      realn=x(no)
      call pot(energy)
      return
      end subroutine funct1

      subroutine pot(energy)
      use commvar
      implicit none
      real*8 :: ebit,energy
      integer ::z
      energy=dble(0)
      do z=1,n3
        ebit=dble(0)
        call poti(ebit,z)
        energy=energy+ebit
      end do
      return
      end subroutine pot
      
      subroutine poti(energy,a)
      use commvar
      use para
      implicit none
      real*8 :: energy
      integer ::a,as(4),ta,tas(4)
      integer :: i,j,k,l,m
      real*8 :: r(4),o(4,4),ot,angle,length
      real*8 :: r3,r41,r42,r51,r52
      energy=dble(0)
      ta=type(a)
      do j=1,4
           as(j)=neigh(a,j)
           tas(j)=type(neigh(a,j))
      end do
      do j=1,4
        do k=1,4
         o(j,k)=dble(0)
         if (j/=k)  o(j,k)=angle(as(j),a,as(k))-o0
        end do
      end do
      do i=1,4
           r(i)=length(a,as(i))-r0(ta,tas(i))
      end do
!		fr term
      do i=1,4
      energy=energy+dble(.25)*(bfr(ta,tas(i))**2)*(r(i)**2)
      end do
!               fo term
      do i=1,3
      do j=i+1,4
      energy=energy+dble(.5)*bfo(tas(i),ta,tas(j))*((r3(as(i),a,as(j))*&
                 o(i,j))**2)
      end do
      end do
!		frr term
      do i=1,3
      do j=i+1,4
      energy=energy+bfrr(tas(i),ta,tas(j))*r(i)*r(j)
      end do
      end do
!		fro term
      do i=1,3
      do j=i+1,4
      energy=energy+r3(as(i),a,as(j))*&
                      bfro(tas(i),ta,tas(j))*o(i,j)*(r(i)+r(j))
      end do
      end do
!		foo term
      do i=1,2
      do j=i+1,3
      do k=j+1,4
      energy=energy+(r41(a,as(i),as(j),as(k))**2)*(&
              bfoo(ta,tas(i),tas(j),tas(k))*o(i,j)*o(j,k)+&
              bfoo(ta,tas(j),tas(k),tas(i))*o(j,k)*o(k,i)+&
              bfoo(ta,tas(k),tas(i),tas(j))*o(k,i)*o(i,j))
      end do
      end do
      end do
!		foo1 term
      do i=1,4
      do j=1,4
      if (i/=j) then
         ot=angle(a,as(j),neigh(as(j),i))-o0
         energy=energy+(r42(as(i),a,as(j),neigh(as(j),i))**2)*o(i,j)*&
         bfoo1(tas(i),ta,tas(j),type(neigh(as(j),i)))*ot*dble(.5)
      end if
      end do
      end do 
!               foo4 term
       do i=1,4
       do j=1,4
       do k=1,4
       if (i/=j.and.i/=k.and.j/=k) then
       ot=angle(as(j),neigh(as(j),i),neigh(neigh(as(j),i),k))-o0
       energy=energy+(r51(as(i),a,as(j),neigh(as(j),i),&
       neigh(neigh(as(j),i),k))**2)*o(i,j)*&
       bfoo4(tas(i),ta,tas(j),type(neigh(as(j),i)),&
             type(neigh(neigh(as(j),i),k)))*ot
       end if
       end do
       end do
       end do
!               foo5 term
       do i=1,4
       do j=1,4
       do k=1,4
       do l=1,4
       if (i/=j.and.i/=k.and.j/=k.and.i/=l.and.j/=l.and.k/=l) then
       ot=angle(as(j),neigh(as(j),k),neigh(neigh(as(j),k),l))-o0
       energy=energy+(r51(as(i),a,as(j),neigh(as(j),k),&
                          neigh(neigh(as(j),k),l))**2)*dble(.5)*&
       bfoo5(tas(i),ta,tas(j),type(neigh(as(j),k)),&
             type(neigh(neigh(as(j),k),l)))*o(i,j)*ot
       end if
       end do
       end do
       end do
       end do
!		foo7 term
      do i=1,4
      do j=1,4
      do k=1,4
      do l=1,4
       if (i/=j.and.i/=k.and.j/=k.and.i/=l.and.j/=l.and.k/=l) then
         ot=angle(a,as(k),neigh(as(k),l))-o0
         energy=energy+(r52(as(i),as(j),as(k),a,neigh(as(k),l))**2)*o(i,j)*ot*&
         bfoo7(tas(i),tas(j),tas(k),ta,type(neigh(as(k),l)))*dble(.5)
       end if
      end do
      end do
      end do
      end do
      return
      end subroutine poti

      real*8 function length(a,b)
      use commvar
      implicit none
      integer a,b,c,d
      real*8 :: dx(3),wrap2
      do d=1,3
        dx(d)=dble(0)
        dx(d)=wrap2(pos(a,d),dpos(a,d),pos(b,d),dpos(b,d))
      end do
      length=sqrt(dx(1)**2+dx(2)**2+dx(3)**2)
      end function length

      real*8 function wrap2(a1,a2,b1,b2)
      use commvar
      implicit none
      real*8 :: a1,a2,b1,b2
      wrap2=a1+a2-b1-b2
      if (a1==dble(0).and.b1==(dble(4*n-1)*r0(1,1)/sqrt(dble(3))))&
                    wrap2=a1+a2-b1-b2+realn
      if (b1==dble(0).and.a1==(dble(4*n-1)*r0(1,1)/sqrt(dble(3))))&
                    wrap2=a1+a2-b1-b2-realn
      end function wrap2

      real*8 function angle(z,b,c)
      use commvar
      implicit none
      integer :: z,b,c,d
      real*8 :: l1,l2,wrap2,v1(3),v2(3),length
      l1=length(z,b)
      l2=length(b,c)
      do d=1,3
        v1(d)=wrap2(pos(z,d),dpos(z,d),pos(b,d),dpos(b,d))
        v2(d)=wrap2(pos(c,d),dpos(c,d),pos(b,d),dpos(b,d))
      end do
      angle=acos((v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3))/(l1*l2))
      if (angle<0) angle=-angle
      if (angle>pi) angle=2*pi-angle
      end function angle

      real*8 function r3(a,b,c)
      use commvar
      implicit none
      integer a,b,c,ta,tb,tc
      ta=type(a)
      tb=type(b)
      tc=type(c)
      r3=(r0(ta,tb)+r0(tb,tc))/dble(2)
      end function r3

      real*8 function r41(a,b,c,d)
      use commvar
      implicit none
      integer a,b,c,d,ta,tb,tc,td
      ta=type(a)
      tb=type(b)
      tc=type(c)
      td=type(d)
      r41=(r0(ta,tb)+r0(ta,tc)+r0(ta,td))/dble(3)
      end function r41

      real*8 function r42(a,b,c,d)
      use commvar
      implicit none
      integer a,b,c,d,ta,tb,tc,td
      ta=type(a)
      tb=type(b)
      tc=type(c)
      td=type(d)
      r42=(r0(ta,tb)+r0(tb,tc)+r0(tc,td))/dble(3)
      end function r42

      real*8 function r51(a,b,c,d,e)
      use commvar
      implicit none
      integer a,b,c,d,e,ta,tb,tc,td,te
      ta=type(a)
      tb=type(b)
      tc=type(c)
      td=type(d)
      te=type(e)
      r51=(r0(ta,tb)+r0(tb,tc)+r0(tc,td)+r0(td,te))/dble(4)
      end function r51

      real*8 function r52(a,b,c,d,e)
      use commvar
      implicit none
      integer a,b,c,d,e,ta,tb,tc,td,te
      ta=type(a)
      tb=type(b)
      tc=type(c)
      td=type(d)
      te=type(e)
      r52=(r0(ta,td)+r0(tb,td)+r0(tc,td)+r0(tc,te))/dble(4)
      end function r52

      subroutine writeout(a)
      use para
      implicit none
      integer :: a,i,j,k,l,m
      if (a==1)open(unit=22,file='start.f90',status='replace')
      if (a==2)open(unit=22,file='param.f90',status='replace')
      if (a==3)open(unit=22,file='minim.f90',status='replace')
      write(22,*)'      module para'
      write(22,*)'      real*8::bfr(0:3,0:3),bfo(0:3,0:3,0:3)'
      write(22,*)'      real*8::bfro(0:3,0:3,0:3),bfrr(0:3,0:3,0:3)'
      write(22,*)'      real*8::bfoo(0:3,0:3,0:3,0:3)'
      write(22,*)'      real*8::bfoo1(0:3,0:3,0:3,0:3)'
      write(22,*)'      real*8::bfoo4(0:3,0:3,0:3,0:3,0:3)'
      write(22,*)'      real*8::bfoo5(0:3,0:3,0:3,0:3,0:3)'
      write(22,*)'      real*8::bfoo7(0:3,0:3,0:3,0:3,0:3)'
      write(22,*)'      contains'
      write(22,*)'      subroutine parameters'
1     format(a10,i1,a1,i1,a2,f13.8)
      do i=0,3
        do j=0,3
          write(22,1)'      bfr(',i,',',j,')=',bfr(i,j)**2  
        end do
      end do
2     format(a10,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
            write(22,2)'      bfo(',i,',',j,',',k,')=',bfo(i,j,k)
           end do
         end do
      end do
3     format(a11,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
             write(22,3)'      bfro(',i,',',j,',',k,')=',bfro(i,j,k)
           end do
         end do
      end do
      do i=0,3
         do j=0,3
           do k=0,3
            write(22,3)'      bfrr(',i,',',j,',',k,')=',bfrr(i,j,k)
           end do
         end do
      end do
4     format(a11,i1,a1,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
              write(22,4)'      bfoo(',i,',',j,',',k,',',l,')=',bfoo(i,j,k,l)
             end do 
           end do
         end do
      end do
5     format(a12,i1,a1,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
              write(22,5)'      bfoo1(',i,',',j,',',k,',',l,')=',bfoo1(i,j,k,l)
             end do
           end do
         end do
      end do
6     format(a12,i1,a1,i1,a1,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
               do m=0,3
      write(22,6)'      bfoo4(',i,',',j,',',k,',',l,',',m,')=',bfoo4(i,j,k,l,m)
               end do
             end do
           end do
         end do
      end do
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
               do m=0,3
      write(22,6)'      bfoo5(',i,',',j,',',k,',',l,',',m,')=',bfoo5(i,j,k,l,m)
               end do
             end do
           end do
         end do
      end do
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
               do m=0,3
      write(22,6)'      bfoo7(',i,',',j,',',k,',',l,',',m,')=',bfoo7(i,j,k,l,m)
               end do
             end do 
           end do
         end do
      end do
      write(22,*)'      end subroutine parameters'
      write(22,*)'      end module para'
      close(22)
      return
      end subroutine writeout

      subroutine ek(lvm)
      use commvar
      implicit none
      real*8 :: lvm
      real*8 :: rr(n1),mass,br(n1,n1),empty(n1,n1),dl(n1),e(n1)
      integer :: ifail,iter(n1),i,j,k,l,m,o,p
      lvm=dble(0)
      do j=1,n1
        rr(j)=dble(0)
        dl(j)=dble(0)
        e(j)=dble(0)
      end do
      do l=1,n3
       do m=1,n3
        do o=1,3
         do p=1,3
          br(3*(l-1)+o,3*(m-1)+p)=dble(0)
          if (l==m.and.o==p) br(3*(l-1)+o,3*(m-1)+p)=mass(type(m))
         end do
        end do
       end do
      end do
      ifail=0
      call f02aef(dv,n1,br,n1,n1,rr,empty,n1,dl,e,ifail)
      if (ifail/=0) write(10,*) 'EK error Ifail=',ifail
      lvm=sqrt(rr(n1))
      return
      end subroutine ek

      subroutine makedv
      use commvar
      implicit none
      real*8 :: swap1,swap2,e1,e2,e3,e4,d,dd
      integer :: i,j,k,l,m
      do i=1,n1
        do j=1,n1
          dv(i,j)=dble(0)
        end do
      end do
      write(10,*) 'Into makedv'
      do i=1,n3
        do j=i,n3
          do k=1,n3
            if (flag(k,i).and.flag(k,j)) then
              if (i==j) then
                do l=1,3
                  do m=l,3
                   if (l==m) then
                    swap1=dpos(i,l)
                    call poti(e1,k)
                    dpos(i,l)=swap1+small+small
                    call poti(e2,k)
                    dpos(i,l)=swap1-small-small
                    call poti(e3,k)
                    dpos(i,l)=swap1
                    dv(3*(i-1)+l,3*(i-1)+l)=dv(3*(i-1)+l,3*(i-1)+l)+&
                      (e2-e1-e1+e3)/(dble(4)*small*small)
                   else
                    d=dd(k,i,l,i,m)
                    dv(3*(i-1)+l,3*(i-1)+m)=dv(3*(i-1)+l,3*(i-1)+m)+d
                    dv(3*(i-1)+m,3*(i-1)+l)=dv(3*(i-1)+m,3*(i-1)+l)+d
                   end if
                  end do
                end do
              else
                do l=1,3
                  do m=1,3
                    d=dd(k,i,l,j,m)
                    dv(3*(i-1)+l,3*(j-1)+m)=dv(3*(i-1)+l,3*(j-1)+m)+d 
                    dv(3*(j-1)+m,3*(i-1)+l)=dv(3*(j-1)+m,3*(i-1)+l)+d
                  end do
                end do
              end if
            end if
          end do
        end do
      end do
      return
      end subroutine makedv
         
      real*8 function dd(a,b,c,d,e)
      use commvar
      implicit none
      integer a,b,c,d,e
      real*8 :: swap1,swap2,e1,e2,e3,e4
      swap1=dpos(b,c)
      swap2=dpos(d,e)
      dpos(b,c)=swap1+small
      dpos(d,e)=swap2+small
      call poti(e1,a)
      dpos(b,c)=swap1-small
      dpos(d,e)=swap2+small
      call poti(e2,a)
      dpos(b,c)=swap1+small
      dpos(d,e)=swap2-small
      call poti(e3,a)
      dpos(b,c)=swap1-small
      dpos(d,e)=swap2-small
      call poti(e4,a)
      dpos(b,c)=swap1
      dpos(d,e)=swap2
      dd=(e1-e2-e3+e4)/(dble(4)*small*small)
      end function dd

      real*8 function mass(i)
      implicit none
      integer :: i
      if (i==1) mass=dble(28.09)
      if (i==2) mass=dble(72.61)
      if (i==3) mass=dble(12.01)
      end function mass