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

!	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
!	neigh(i,j) = j th neighbour of ith atom
!	flag(i,j) = relationship between atom i and atom j
!	r0(i,j)= perfect lattice length between atoms of type=i and type=j
!	o0= 109'
!	pi 
!	sf= scale factor up to angstroms
!	wrap = boundary conditions function
!	dv= double differntial of potential matrix
!	type = ith atom is type i 1=Si 2=Ge 3=C
!	go5caf,rnd random number generator

      contains

      subroutine values
      implicit none
      integer :: i,j
      base=3
      twos=0
      small=dble(.0001)
      r0(1,1)=dble(2.35156)
      r0(2,2)=dble(2.447)
      r0(3,3)=dble(1.545)
      r0(1,3)=dble(1.888)
      r0(3,1)=r0(1,2)
      r0(2,3)=(r0(2,2)+r0(3,3))/dble(2)
      do i=1,3
        do j=1,3
           r0(i,j)=r0(i,j)*(sqrt(dble(3))/dble(4))/r0(base,base)
        end do 
      end do
      o0=acos(dble(-1)/dble(3))
      pi=acos(dble(-1))
      return
      end subroutine values

      end module commvar

!	MAIN PROGRAM

      program phd1
      use commvar
      use para
      call values
      call parameters
      call g05ccf
      call setup
      call relax
      call makedv
      call ek(1,1,1)
!     call ek(1,1,0)
!     call ek(1,0,0)
      end program phd1

      subroutine setup
      use commvar
      implicit none
      integer :: count,set=n*n*n,alter
      integer :: i,j,k,l,m
      real*8 :: match1,match2,match3
      real*8:: g05caf,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
!	Variable flag
!	f = no connection
!	t = connection 
      do i=1,n3
            do j=1,n3
                 if(i==j) flag(i,j)=.true.
                 if(i/=j) flag(i,j)=.false.
             end do
      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,4
      do k=1,4
       if (j/=k) then
       if (n>1) flag(i,neigh(neigh(i,j),k))=.true.
      do l=1,4
      if(l/=j.and.l/=k) then
!   CHANGE
            if(n>1)flag(i,neigh(neigh(neigh(i,k),j),l))=.true.
!   CHANGE
            if(n>1)flag(i,neigh(neigh(i,k),l))=.true.
                 do m=1,4
                 if(m/=j.and.m/=k.and.m/=l) then
!   CHANGE
              if(n>1) flag(i,neigh(neigh(neigh(i,k),l),m))=.true.
              if(n>1)  flag(i,neigh(neigh(i,l),m))=.true.
                 end if
                 end do
      end if
      end do
      end if
      end do
      end do
      end do
      do i=1,n3
      type(i)=base
      end do
      do i=1,twos
20		alter=int(g05caf(rnd)*dble(n3))+1
                   if (type(alter)==2) goto 20
                   type(alter)=2
      end do
!      do i=(n3/2)+1,n3
!      type(i)=2
!      end do
      open(unit=11,file='nn',status='replace')
      write(11,*) neigh
      close(11)
      open(unit=12,file='perfectpos',status='replace')
      write(12,*) pos
      close(12)
      open(unit=14,file='type',status='replace')
      write(14,*) type
      close(14)
      open(unit=15,file='flag',status='replace')
      write(15,*) flag
      close(15)
      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
      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,n2
      bl(i)=dble(-0.2)
      x(i)=dble(0)
      ul(i)=dble(0.2)
      end do
      bl(n2)=realn-dble(0.2)
      x(n2)=realn
      ul(n2)=realn+dble(0.2)
      ifail=1
      ibound=0
      output=.true.
      call e04jaf(n2,ibound,bl,ul,x,f,iw,liw,w,lw,ifail)
      output=.false.
      write(21,*) 'Final value =',f
      do i=1,n3
            do j=1,3
                 dpos(i,j)=x(3*(i-1)+j)
             end do
       end do
       realn=x(n2)
      write(21,*) 'final realn =',realn
      open(unit=15,file='dpos',status='replace')
      write(15,*) dpos
      close(15)
      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
       if (output) write(21,*) energy,realn
       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(as(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))*(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(tas(i),ta,tas(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(tas(i),ta,tas(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(ta,tas(i),tas(j),tas(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
        if (n==1) goto 1000
!		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(tas(i),ta,tas(j),type(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
!   CHANGE
!       if (n==2) goto 2000
        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(tas(i),ta,tas(j),type(neigh(as(j),i)),&
!   CHANGE
!       type(neigh(neigh(as(j),i),k)))**2)*o(i,j)*dble(.5)*&
        type(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(tas(i),ta,tas(j),type(neigh(as(j),k)),&
                           type(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
2000     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(tas(i),tas(j),tas(k),ta,type(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
1000        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),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,d)
          use commvar
          implicit none
          real*8 :: a1,a2,b1,b2,swap
          integer ::d
          wrap2=a1+a2-b1-b2
          if (a1==dble(0).and.b1==dble(n)-dble(0.25))&
                    wrap2=a1+a2-b1-b2+realn
          if (b1==dble(0).and.a1==dble(n)-dble(0.25))&
                    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),d)
           v2(d)=wrap2(pos(c,d),dpos(c,d),pos(b,d),dpos(b,d),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
            r3=(r0(a,b)+r0(b,c))/dble(2)
            end function r3

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

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

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

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

          subroutine makedv
          use commvar
          implicit none
          real*8 :: swap1,swap2,e1,e2,e3,e4,d,dd
          real*8 :: length
          integer ::phaseik(3),phasekj(3),tph(3)
          integer :: i,j,k,l,m,o
          do i=1,n1
          do j=1,n1
          do k=-1,1
          do l=-1,1
          do m=-1,1
               dv(i,j,k,l,m)=dble(0)
          end do
          end do
          end do
          end do
          end do
          write(21,*) ''
          write(21,*) 'Into makedv'
          write(21,*) ''
          open(unit=44,file='dv',status='replace')
          do i=1,n3
          write(21,*) i
          do j=i,n3
          do k=1,n3
              if (flag(k,i).and.flag(k,j)) then
              call phasedif(i,k,phaseik)
              call phasedif(k,j,phasekj)
              do l=1,3
              tph(l)=phaseik(l)+phasekj(l)
              end do
                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,tph(1),tph(2),tph(3))=&
                      dv(3*(i-1)+l,3*(i-1)+l,tph(1),tph(2),tph(3))+&
                      (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,tph(1),tph(2),tph(3))=&
                      dv(3*(i-1)+l,3*(i-1)+m,tph(1),tph(2),tph(3))+d
                    dv(3*(i-1)+m,3*(i-1)+l,tph(1),tph(2),tph(3))=&
                      dv(3*(i-1)+m,3*(i-1)+l,tph(1),tph(2),tph(3))+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,tph(1),tph(2),tph(3))=&
                        dv(3*(i-1)+l,3*(j-1)+m,tph(1),tph(2),tph(3))+d
                     dv(3*(j-1)+m,3*(i-1)+l,-tph(1),-tph(2),-tph(3))=&
                        dv(3*(j-1)+m,3*(i-1)+l,-tph(1),-tph(2),-tph(3))+d
                   end do
                   end do
                end if
              end if
         end do
         end do
         end do
         do i=1,n1
         do j=1,n1
         do k=-1,1
         do l=-1,1
         do m=-1,1
         if (dv(i,j,k,l,m)/=dble(0)) write(44,*) i,j,k,l,m,dv(i,j,k,l,m)
         end do
         end do
         end do
         end do
         end do
         write(44,*) 0,0,0,0,0,dble(0)
         close(44)
         return
         end subroutine makedv
         
         subroutine phasedif(z,y,x)
         use commvar
         implicit none
         real*8 :: dif
         integer :: z,y,x(3),w
         do w=1,3
         x(w)=0
         dif=abs(pos(z,w)-pos(y,w))
         if (n==1.and.dif>dble(.25)) then
         if (pos(z,w)<pos(y,w)) x(w)=1
         if (pos(y,w)<pos(z,w)) x(w)=-1
         end if
!   CHANGE
!        if (n==2.and.dif>dble(.5)) then
!        if (pos(z,w)<pos(y,w)) x(w)=1
!        if (pos(y,w)<pos(z,w)) x(w)=-1
!        end if 
         if (n>1.and.dif>dble(.75)) then
         if (pos(z,w)<pos(y,w)) x(w)=1
         if (pos(y,w)<pos(z,w)) x(w)=-1
         end if
         end do
         return
         end subroutine phasedif


         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

         subroutine ek(k1,k2,k3)
          use commvar
          implicit none
          integer :: k1,k2,k3,o,p,q,r,s,file
          real*8 :: ai(n1,n1),ar(n1,n1),rr(n1),ri(n1),mass,q1,q2,q3
          real*8 :: bi(n1,n1),br(n1,n1),empty(n1,n1),beta(n1)
          integer :: iter(n1),ifail
          integer :: i,j,k,l,m
          write(21,*) ''
          write(21,*) 'Into ek',k1,k2,k3
          write(21,*) ''
          file=n*(k1*100+k2*10+k3)
          do i=0,100
          write(21,*) i
          q1=pi*dble(i)*dble(4*k1)/dble(100)
          q2=pi*dble(i)*dble(4*k2)/dble(100)
          q3=pi*dble(i)*dble(4*k3)/dble(100)
          do j=1,n1
          rr(j)=dble(0)
          ri(j)=dble(0)
          do k=1,n1
          ai(j,k)=dble(0)
          ar(j,k)=dble(0)
          end do
          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)
          bi(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(l))
          do q=-1,1
          do r=-1,1
          do s=-1,1
          ar(3*(l-1)+o,3*(m-1)+p)=ar(3*(l-1)+o,3*(m-1)+p)+&
            dv(3*(l-1)+o,3*(m-1)+p,q,r,s)*cos(q1*dble(q)+q2*dble(r)+&
             q3*dble(s))
          ai(3*(l-1)+o,3*(m-1)+p)=ai(3*(l-1)+o,3*(m-1)+p)+&
             dv(3*(l-1)+o,3*(m-1)+p,q,r,s)*sin(q1*dble(q)+q2*dble(r)+&
             q3*dble(s))
          end do
           end do
           end do
           end do
           end do
           end do
           end do
           ifail=0
           call f02gjf(n1,ar,n1,ai,n1,br,n1,bi,n1,dble(-1),rr,ri,beta,&
                        .false.,empty,n1,empty,n1,iter,ifail)
           if (ifail/=0) write(21,*) 'EK error Ifail=',ifail,q1,q2,q3
          do j=1,n1
                if (abs(ri(j)/beta(j))>1e-6) then
                    write(12,*) q1,q2,q3
                    write(12,*) rr(j)/beta(j),ri(j)/beta(j)
                else
                if (rr(j)/beta(j)<-1e-6) then
                     write(12,*) q1,q2,q3
                     write(12,*) rr(j)/beta(j),ri(j)/beta(j)
                else
                if (rr(j)/beta(j)>dble(0)) write(file,*) i,sqrt(rr(j)/beta(j))
                if (rr(j)/beta(j)<=dble(0)) write(file,*) i,dble(0)
                end if
                end if
             end do
             end do
              return
             end subroutine ek


             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