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),opos(n3,3)
      integer :: neigh(n3,4)
      real*8 :: r0(3,3),o0,pi
      integer :: type(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='mini',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))
      end subroutine values

      end module commvar

!	MAIN PROGRAM

      program phd4
      use commvar
      use para
      implicit none
      call values
      call parameters
      call setup
      call readin
      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)
      a=1
      b=2
      bfr(a,b)=x(1)
      bfr(b,a)=x(1)
      bfo(a,b,b)=x(2)
      bfo(b,b,a)=x(2)
      bfo(b,a,b)=x(3)
      bfo(a,a,b)=x(4)
      bfo(b,a,a)=x(4)
      bfo(a,b,a)=x(5)
      bfro(a,b,b)=x(6)
      bfro(b,b,a)=x(6)
      bfro(b,a,b)=x(7)
      bfro(a,a,b)=x(8)
      bfro(b,a,a)=x(8)
      bfro(a,b,a)=x(9)
      bfrr(a,b,b)=x(10)
      bfrr(b,b,a)=x(10)
      bfrr(b,a,b)=x(11)
      bfrr(a,a,b)=x(12)
      bfrr(b,a,a)=x(12)
      bfrr(a,b,a)=x(13)
      bfoo(a,b,b,b)=x(14)
      bfoo(a,a,b,b)=x(15)
      bfoo(a,b,b,a)=x(15)
      bfoo(a,a,a,b)=x(16)
      bfoo(a,b,a,a)=x(16)
      bfoo(a,b,a,b)=x(17)
      bfoo(a,a,b,a)=x(18)
      bfoo(b,a,a,a)=x(19)
      bfoo(b,a,b,b)=x(20)
      bfoo(b,b,b,a)=x(20)
      bfoo(b,b,a,b)=x(21)
      bfoo(b,a,b,a)=x(22)
      bfoo(b,a,a,a)=x(23)
      bfoo1(b,a,b,b)=x(24)
      bfoo1(b,b,a,b)=x(24)
      bfoo1(a,a,b,b)=x(25)
      bfoo1(b,b,a,a)=x(25)
      bfoo1(a,b,b,b)=x(26)
      bfoo1(b,b,b,a)=x(26)
      bfoo1(b,a,a,a)=x(27)
      bfoo1(a,a,a,b)=x(27)
      bfoo1(a,a,b,a)=x(28)
      bfoo1(a,b,a,a)=x(28)
      bfoo1(b,a,b,a)=x(29)
      bfoo1(a,b,a,b)=x(29)
      bfoo1(a,b,b,a)=x(30)
      bfoo1(b,a,a,b)=x(31)
      bfoo7(a,a,a,a,b)=x(32)
      bfoo7(a,a,a,b,a)=x(33)
      bfoo7(a,a,a,b,b)=x(34)
      bfoo7(a,a,b,a,a)=x(35)
      bfoo7(a,a,b,a,b)=x(36)
      bfoo7(a,a,b,b,a)=x(37)
      bfoo7(a,a,b,b,b)=x(38)
      bfoo7(a,b,a,a,a)=x(39)
      bfoo7(b,a,a,a,a)=x(39)
      bfoo7(a,b,a,a,b)=x(40)
      bfoo7(b,a,a,a,b)=x(40)
      bfoo7(a,b,a,b,a)=x(41)
      bfoo7(b,a,a,b,a)=x(41)
      bfoo7(a,b,a,b,b)=x(42)
      bfoo7(b,a,a,b,b)=x(42)
      bfoo7(a,b,b,a,a)=x(43)
      bfoo7(b,a,b,a,a)=x(43)
      bfoo7(a,b,b,a,b)=x(44)
      bfoo7(b,a,b,a,b)=x(44)
      bfoo7(a,b,b,b,a)=x(45)
      bfoo7(b,a,b,b,a)=x(45)
      bfoo7(a,b,b,b,b)=x(46)
      bfoo7(b,a,b,b,b)=x(46)
      bfoo7(b,b,a,a,a)=x(47)
      bfoo7(b,b,a,a,b)=x(48)
      bfoo7(b,b,a,b,a)=x(49)
      bfoo7(b,b,a,b,b)=x(50)
      bfoo7(b,b,b,a,a)=x(51)
      bfoo7(b,b,b,a,b)=x(52)
      bfoo7(b,b,b,b,a)=x(53)
      call relax
      do i=2,n3
         do j=1,3
           dpos(i,j)=dpos(i,j)-dpos(1,j)
         end do
      end do
      dpos(1,1)=dble(0)
      dpos(1,2)=dble(0)
      dpos(1,3)=dble(0)
      do i=1,n3
        do j=1,3
          fvecc(3*(i-1)+j)=opos(i,j)-pos(i,j)-dpos(i,j)
          sum=sum+fvecc(3*(i-1)+j)**2
        end do
      end do
      if (sum<mini) then
        mini=sum
        write(101,*) mini
        open(unit=102,file='mini',status='replace')
        write(102,*) mini
        close(102)
        call writeout(3)
      end if
      write(10,*) 'sum=',sum
      write(10,*)
      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(100),lw,ifail
      real*8 :: eta,xtol,stepmx,x(53),fsumsq,fvec(192),fjac(192,53)
      real*8 :: s(53),v(53,53),w(12256)
      real*8 :: fr,fo,frr,fro,foo,foo1,foo7
      integer :: i,j,k,l,m,a,b
      external lsqfun,lsqmon
      ifail=1
      m1=192
      no=53
      liw=100
      lj=m1
      lv=no
      lw=6*no+m1*no+2*m1+no*(no-1)/2
      iprint=1
      maxcal=20000
      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))
          do k=1,3
          bfo(i,j,k)=sqrt(bfo(i,j,k))
          end do
        end do
      end do
      a=1
      b=2
      x(1)=bfr(a,b)
      x(2)=bfo(a,b,b)
      x(3)=bfo(b,a,b)
      x(4)=bfo(a,a,b)
      x(5)=bfo(a,b,a)
      x(6)=bfro(a,b,b)
      x(7)=bfro(b,a,b)
      x(8)=bfro(a,a,b)
      x(9)=bfro(a,b,a)
      x(10)=bfrr(a,b,b)
      x(11)=bfrr(b,a,b)
      x(12)=bfrr(a,a,b)
      x(13)=bfrr(a,b,a)
      x(14)=bfoo(a,b,b,b)
      x(15)=bfoo(a,b,b,a)
      x(16)=bfoo(a,a,a,b)
      x(17)=bfoo(a,b,a,b)
      x(18)=bfoo(a,a,b,a)
      x(19)=bfoo(b,a,a,b)
      x(20)=bfoo(b,a,b,b)
      x(21)=bfoo(b,b,a,b)
      x(22)=bfoo(b,a,b,a)
      x(23)=bfoo(b,a,a,a)
      x(24)=bfoo1(b,a,b,b)
      x(25)=bfoo1(a,a,b,b)
      x(26)=bfoo1(a,b,b,b)
      x(27)=bfoo1(a,a,a,b)
      x(28)=bfoo1(a,a,b,a)
      x(29)=bfoo1(b,a,b,a)
      x(30)=bfoo1(a,b,b,a)
      x(31)=bfoo1(b,a,a,b)
      x(32)=bfoo7(a,a,a,a,b)
      x(33)=bfoo7(a,a,a,b,a)
      x(34)=bfoo7(a,a,a,b,b)
      x(35)=bfoo7(a,a,b,a,a)
      x(36)=bfoo7(a,a,b,a,b)
      x(37)=bfoo7(a,a,b,b,a)
      x(38)=bfoo7(a,a,b,b,b)
      x(39)=bfoo7(a,b,a,a,a)
      x(40)=bfoo7(a,b,a,a,b)
      x(41)=bfoo7(a,b,a,b,a)
      x(42)=bfoo7(a,b,a,b,b)
      x(43)=bfoo7(a,b,b,a,a)
      x(44)=bfoo7(a,b,b,a,b)
      x(45)=bfoo7(a,b,b,b,a)
      x(46)=bfoo7(a,b,b,b,b)
      x(47)=bfoo7(b,b,a,a,a)
      x(48)=bfoo7(b,b,a,a,b)
      x(49)=bfoo7(b,b,a,b,a)
      x(50)=bfoo7(b,b,a,b,b)
      x(51)=bfoo7(b,b,b,a,a)
      x(52)=bfoo7(b,b,b,a,b)
      x(53)=bfoo7(b,b,b,b,a)
      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(a,b)=x(1)
      bfr(b,a)=x(1)
      bfo(a,b,b)=x(2)
      bfo(b,b,a)=x(2)
      bfo(b,a,b)=x(3)
      bfo(a,a,b)=x(4)
      bfo(b,a,a)=x(4)
      bfo(a,b,a)=x(5)
      bfro(a,b,b)=x(6)
      bfro(b,b,a)=x(6)
      bfro(b,a,b)=x(7)
      bfro(a,a,b)=x(8)
      bfro(b,a,a)=x(8)
      bfro(a,b,a)=x(9)
      bfrr(a,b,b)=x(10)
      bfrr(b,b,a)=x(10)
      bfrr(b,a,b)=x(11)
      bfrr(a,a,b)=x(12)
      bfrr(b,a,a)=x(12)
      bfrr(a,b,a)=x(13)
      bfoo(a,b,b,b)=x(14)
      bfoo(a,a,b,b)=x(15)
      bfoo(a,b,b,a)=x(15)
      bfoo(a,a,a,b)=x(16)
      bfoo(a,b,a,a)=x(16)
      bfoo(a,b,a,b)=x(17)
      bfoo(a,a,b,a)=x(18)
      bfoo(b,a,a,a)=x(19)
      bfoo(b,a,b,b)=x(20)
      bfoo(b,b,b,a)=x(20)
      bfoo(b,b,a,b)=x(21)
      bfoo(b,a,b,a)=x(22)
      bfoo(b,a,a,a)=x(23)
      bfoo1(b,a,b,b)=x(24)
      bfoo1(b,b,a,b)=x(24)
      bfoo1(a,a,b,b)=x(25)
      bfoo1(b,b,a,a)=x(25)
      bfoo1(a,b,b,b)=x(26)
      bfoo1(b,b,b,a)=x(26)
      bfoo1(b,a,a,a)=x(27)
      bfoo1(a,a,a,b)=x(27)
      bfoo1(a,a,b,a)=x(28)
      bfoo1(a,b,a,a)=x(28)
      bfoo1(b,a,b,a)=x(29)
      bfoo1(a,b,a,b)=x(29)
      bfoo1(a,b,b,a)=x(30)
      bfoo1(b,a,a,b)=x(31)
      bfoo7(a,a,a,a,b)=x(32)
      bfoo7(a,a,a,b,a)=x(33)
      bfoo7(a,a,a,b,b)=x(34)
      bfoo7(a,a,b,a,a)=x(35)
      bfoo7(a,a,b,a,b)=x(36)
      bfoo7(a,a,b,b,a)=x(37)
      bfoo7(a,a,b,b,b)=x(38)
      bfoo7(a,b,a,a,a)=x(39)
      bfoo7(b,a,a,a,a)=x(39)
      bfoo7(a,b,a,a,b)=x(40)
      bfoo7(b,a,a,a,b)=x(40)
      bfoo7(a,b,a,b,a)=x(41)
      bfoo7(b,a,a,b,a)=x(41)
      bfoo7(a,b,a,b,b)=x(42)
      bfoo7(b,a,a,b,b)=x(42)
      bfoo7(a,b,b,a,a)=x(43)
      bfoo7(b,a,b,a,a)=x(43)
      bfoo7(a,b,b,a,b)=x(44)
      bfoo7(b,a,b,a,b)=x(44)
      bfoo7(a,b,b,b,a)=x(45)
      bfoo7(b,a,b,b,a)=x(45)
      bfoo7(a,b,b,b,b)=x(46)
      bfoo7(b,a,b,b,b)=x(46)
      bfoo7(b,b,a,a,a)=x(47)
      bfoo7(b,b,a,a,b)=x(48)
      bfoo7(b,b,a,b,a)=x(49)
      bfoo7(b,b,a,b,b)=x(50)
      bfoo7(b,b,b,a,a)=x(51)
      bfoo7(b,b,b,a,b)=x(52)
      bfoo7(b,b,b,b,a)=x(53)
      return
      end subroutine leastsquares

      subroutine setup
      use commvar
      implicit none
      integer :: count,set=n*n*n,alter
      integer :: i,j,k,l,m
      real*8 :: match1,match2,match3,rnd,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
      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
        if (c0>n3) then
          print *,'cock up',c1,c2,c3
          stop 
        end if
        goto 10
      end if
      neigh(b,a)=c0
      neigh(c0,a)=b
      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))**2)*((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 readin
       use commvar
       implicit none
       integer ::a,i,j,k,e,f,g,h,place
       real*8 :: b,c,d,opos2(n3,3)
       open(unit=45,file='sige.dat',status='unknown')
100    format(i3,i3,i4,i4,i4,i4,f15.8,f15.8,f14.8)
       do i=1,n3
         read(45,100)a,type(i),j,k,e,f,opos(i,1),opos(i,2),opos(i,3)
         dpos(i,1)=opos(i,1)-pos(i,1)
         dpos(i,2)=opos(i,2)-pos(i,2)
         dpos(i,3)=opos(i,3)-pos(i,3)
!        if (type(i)==2) type(i)=3
       end do
       open(unit=45,file='sige.txt',status='unknown')
       do i=1,n3
         read(45,*)a,opos2(i,1),opos2(i,2),opos2(i,3),b,c,d
       end do
       place=1
       do i=1,43
10       if (type(place)==1) then
           opos(place,1)=opos2(i,1)
           opos(place,2)=opos2(i,2)
           opos(place,3)=opos2(i,3)
           place=place+1
         else
           place=place+1
           goto 10
         end if 
       end do
       place=1
       do i=44,64
20       if (type(place)==2) then
           opos(place,1)=opos2(i,1)
           opos(place,2)=opos2(i,2)
           opos(place,3)=opos2(i,3)
           place=place+1
         else
           place=place+1
           goto 20
         end if
       end do
       do i=1,n3
         do j=1,3
           opos(i,j)=dble(.529177249)*opos(i,j)
         end do
       end do
       open(unit=101,file='newloc.dat',status='replace')
       do i=1,n3
         write(101,*) i,opos(i,1),opos(i,2),opos(i,3)
       end do
       close(101)    
       do i=2,n3
         do j=1,3
           opos(i,j)=opos(i,j)-opos(1,j)
         end do
       end do
       opos(1,1)=dble(0)
       opos(1,2)=dble(0)
       opos(1,3)=dble(0)
       return
       end subroutine readin

       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)**2
           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