module commvar 
      integer,parameter :: nx=12,ny=12,nz=12,n=nx*ny*nz
      real*8 ::rnx,rny
!  NX & NY MUST BE EVEN 
      integer :: neigh(0:n,4),pos(n,3),type(0:n)
      integer :: testplace,atoms
      real*8 :: rpos(n,3),dpos(n,3),pi,per(3),freq(3,3)
      real*8 :: r0(3,3),o0,cohes(3,3),kt
      logical :: flag(n,n),flag2(n,n),flagto(n),flagfrom(n)
      contains
      subroutine values
      implicit none
      per(1)=dble(1)
      per(2)=dble(0)
      per(3)=dble(0)
      kt=dble(8.617e-5)*dble(600)
      pi=acos(dble(-1))
      freq(1,1)=dble(15.53e12)
      freq(2,2)=dble(9.0e12)
      freq(3,3)=dble(39.95e12)
      freq(1,2)=dble(13.29e12)
      freq(2,1)=freq(1,2)
      freq(1,3)=dble(18.13e12)
      freq(3,1)=freq(1,3)
      freq(2,3)=dble(15.92e12)
      freq(3,2)=freq(2,3)
      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)
      o0=acos(dble(-1)/dble(3))
      pi=acos(dble(-1))
      end subroutine values

      end module commvar

!  MAIN PROGRAM
!   THINGS TO SET & WHERE TO SET THEM
!  
!     per = percentage of each type of atom incoming   line 13-15
!     substrate = number of the type that the substrate is made of  line ****
!     depth of substrate      line
!

      program phd3
      use commvar
      use para
      implicit none
      integer :: goes
      call g05ccf
      call values
      call setup
      call parameters
      goes=0
      do while (atoms>6*nx*ny.and.atoms<8*nx*ny.and.goes<4000)
        call move
        goes=goes+1
      end do
      end program phd3

      subroutine setup
      use commvar
      implicit none
      integer :: i,j,k,l,m,o,place,mn,substrate,nn
      real*8 :: mx,my,mz,sumdis,length
      substrate=1
      rnx=dble(nx)*r0(substrate,substrate)*sqrt(dble(2)/dble(3))
      rny=dble(ny)*r0(substrate,substrate)*sqrt(dble(2))
      type(0)=0
      atoms=0
      do i=0,n
        do j=1,4
           neigh(i,j)=0
        end do
      end do
      do i=1,nx
        do j=1,ny
          do k=1,nz
            place=i+(j-1)*nx+(k-1)*nx*ny
            if (k<8) then
!   don't forget to also put that number in the setting up of flagfrom&flagto
             type(place)=1
             atoms=atoms+1
            end if
            pos(place,1)=i
            pos(place,2)=j
            pos(place,3)=k
            dpos(place,1)=dble(0)
            dpos(place,2)=dble(0)
            dpos(place,3)=dble(0)
             rpos(place,1)=dble(i-1)*sqrt(dble(2)/dble(3))
            if((-1)**(i+j+k)==-1) then
              rpos(place,2)=dble(j)*dsqrt(dble(2))/dble(3)+dble(j-1)*&
                            dsqrt(dble(8))/dble(3)+dble(k-1)*&
                            dsqrt(dble(2))/dble(3)
              rpos(place,3)=dble(k-1)+dble(k)/dble(3)
              mx=i-1
              my=j
              mz=k
              if(mx==0) mx=nx
              mn=mx+(my-1)*nx+(mz-1)*nx*ny
              neigh(place,1)=mn
              mx=i+1
              my=j
              mz=k
              if(mx>nx) mx=1
              mn=mx+(my-1)*nx+(mz-1)*nx*ny
              neigh(place,2)=mn
              mx=i
              my=j+1
              mz=k
              if(my>ny) my=1
              mn=mx+(my-1)*nx+(mz-1)*nx*ny
              neigh(place,3)=mn
              mx=i
              my=j
              mz=k+1
              if(mz>nz) then
                mn=mx+(my-1)*nx
                neigh(place,4)=0
              else
                mn=mx+(my-1)*nx+(mz-1)*nx*ny
                neigh(place,4)=mn
              endif
            else
              rpos(place,2)=sqrt(dble(2))*dble(j-1)+&
                       dble(k-1)*sqrt(dble(2))/dble(3)      
              rpos(place,3)=dble(4)*dble(k-1)/dble(3)
              mx=i+1
              my=j
              mz=k
              if(mx>nx) mx=1
              mn=mx+(my-1)*nx+(mz-1)*nx*ny
              neigh(place,1)=mn
              mx=i-1
              my=j
              mz=k
              if(mx==0) mx=nx
              mn=mx+(my-1)*nx+(mz-1)*nx*ny
              neigh(place,2)=mn
              mx=i
              my=j-1
              mz=k
              if(my==0) my=ny
              mn=mx+(my-1)*nx+(mz-1)*nx*ny
              neigh(place,3)=mn
              mx=i
              my=j
              mz=k-1
              if(mz==0) then
                mn=mx+(my-1)*nx+(nz-1)*nx*ny
                neigh(place,4)=0
              else
                mn=mx+(my-1)*nx+(mz-1)*nx*ny
                neigh(place,4)=mn
              endif
            end if
          end do
        end do
      end do
      do i=1,n
        do j=1,3
          rpos(i,j)=rpos(i,j)*r0(substrate,substrate)
        end do
        write(200,*) neigh(i,1),neigh(i,2),neigh(i,3),neigh(i,4)
        write(200,*) rpos(i,1),rpos(i,2),rpos(i,3)
        write(200,*) pos(i,1),pos(i,2),pos(i,3)
      end do
      do i=1,n
        flagto(i)=.false.
        flagfrom(i)=.false.
        if(pos(i,3)==7.and.pos(neigh(i,4),3)==8) flagfrom(i)=.true.
        if(pos(i,3)==8.and.pos(neigh(i,4),3)==7) flagto(i)=.true.
        do j=1,n
          flag(i,j)=.false.
          flag2(i,j)=.false.
        end do
      end do
      do i=1,n
         do j=1,4
           if(neigh(i,j)/=0) then
             flag(i,neigh(i,j))=.true.
             flag2(i,neigh(i,j))=.true.
             do k=1,4
              if(neigh(neigh(i,j),k)/=0.and.j/=k)then
                flag(i,neigh(neigh(i,j),k))=.true.
                flag2(i,neigh(neigh(i,j),k))=.true.
                do l=1,4
                  if(neigh(neigh(neigh(i,j),k),l)/=0.and.k/=l)then
                    flag(i,neigh(neigh(neigh(i,j),k),l))=.true.
                    flag2(i,neigh(neigh(neigh(i,j),k),l))=.true.
                   do m=1,4
                     if(neigh(neigh(neigh(neigh(i,j),k),l),m)/=0.and.m/=l)then
                       flag2(i,neigh(neigh(neigh(neigh(i,j),k),l),m))=.true.
           do o=1,4
             if(neigh(neigh(neigh(neigh(neigh(i,j),k),l),m),o)/=0.and.m/=o)&
                flag2(i,neigh(neigh(neigh(neigh(neigh(i,j),k),l),m),o))=.true.
           end do 
                     end if
                   end do
                  end if
                end do
              end if 
             end do
           end if
         end do
      end do
      sumdis=dble(0)
      do i=1,n
        if (flag(786,i)) sumdis=sumdis+dble(1)/length(786,i,.true.) 
      end do
      cohes(1,1)=dble(4.64)/sumdis
      cohes(2,2)=(r0(2,2)*dble(3.87))/(r0(1,1)*sumdis)
      cohes(3,3)=(r0(2,2)*dble(7.36))/(r0(1,1)*sumdis)
      do i=1,3
        do j=1,3
          cohes(i,j)=(cohes(i,i)+cohes(j,j))/dble(2)
        end do
      end do
      return
      end subroutine setup


      subroutine funct1(no,x,energy)
      use commvar
      implicit none
      integer :: no
      real*8 :: x(no),energy
      energy=dble(0)
      dpos(testplace,1)=x(1)
      dpos(testplace,2)=x(2)
      dpos(testplace,3)=x(3)
      call pot(energy,testplace)
      return
      end subroutine funct1

      subroutine pot(energy,place)
      use commvar
      implicit none
      integer :: place,i,j
      real*8 :: energy,ebit
      energy=dble(0)
      do i=1,4
        if(type(neigh(place,i))/=0)then
         ebit=dble(0)
         call poti(ebit,neigh(place,i))
         energy=energy+ebit
        end if
        if(neigh(place,i)/=0)then
         do j=1,4
           if(type(neigh(neigh(place,i),j))/=0.and.i/=j)then
            ebit=dble(0)
            call poti(ebit,neigh(neigh(place,i),j))
            energy=energy+ebit
           end if
         end do
        end if
      end do
      ebit=dble(0)
      if (type(place)>0) call poti(ebit,place)
      energy=energy+ebit
      return
      end subroutine pot

      real*8 function angle(a,b,c)
      use commvar
      implicit none
      integer :: a,b,c,d
      real*8 :: l1,l2,wrap,v1(3),v2(3),length
      l1=length(a,b,.false.)
      l2=length(b,c,.false.)
      do d=1,3
        v1(d)=wrap(a,b,d,.false.)
        v2(d)=wrap(c,b,d,.false.)
      end do
      angle=acos((v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3))/(l1*l2))
      if (angle<dble(0)) angle=-angle
      if (angle>acos(dble(-1))) angle=2*acos(dble(-1))-angle
      angle=angle-o0
      end function angle

      real*8 function length(a,b,f)
      use commvar
      implicit none
      integer a,b,d,i,j
      logical :: f
      real*8 :: dx(3),wrap
      do d=1,3
        dx(d)=dble(0)
        dx(d)=wrap(a,b,d,f)
      end do
      length=dsqrt(dx(1)**2+dx(2)**2+dx(3)**2)
      end function length

      real*8 function wrap(a,b,d,f)
      use commvar
      implicit none
      integer :: a,b,d
      logical :: f
      if (f) then
        wrap=abs(rpos(a,d)-rpos(b,d))
        if (d==1) then
          if(abs(rnx-wrap)<wrap) wrap=abs(rnx-wrap)
        end if
        if (d==2) then
          if(abs(rny-wrap)<wrap) wrap=abs(rny-wrap)
        end if
      else
        wrap=rpos(a,d)+dpos(a,d)-rpos(b,d)-dpos(b,d)
        if (d==1) then
          if(pos(a,1)==1.and.pos(b,1)==nx)&
               wrap=wrap+rnx
          if(pos(a,1)==nx.and.pos(b,1)==1)&
               wrap=wrap-rnx
        end if
        if (d==2) then
          if(pos(a,2)==1.and.pos(b,2)==ny)&
               wrap=wrap+rny
          if(pos(a,2)==ny.and.pos(b,2)==1)&
               wrap=wrap-rny
        end if
      end if
      end function wrap

      subroutine calb(eb,p,to)
      use commvar
      implicit none
      real*8 :: eb,length
      integer :: p,i,to
      logical :: empty 
      eb=dble(0)
      empty=.false.
      if (type(p)==0) then
        type(p)=to
        empty=.true.
      end if
      do i=1,n
        if (flag(p,i).and.type(i)>0)then
          eb=eb+cohes(type(p),type(i))/length(p,i,.true.)
        end if
      end do
      if (empty) type(p)=0
      return
      end subroutine calb
      

      subroutine move
      use commvar
      implicit none
      integer,parameter ::tn=3,liw=5,lw=39
      integer :: i,j,k,l,m,ibound,ifail,nn
      integer :: iw(liw),bestplace,bestto,bestplace2
      real*8 :: bpos(3),tpos(3),besttime,testtime
      real*8 :: bl(3),ul(3),w(lw),g05caf,rnd
      real*8 :: ebit1,ebit2,ebitb,length,freqcal
      besttime=dble(10000)
      bestplace=0
      bestplace2=0
      bestto=0
      do i=1,3
        bpos(i)=dble(0)
      end do
      do i=1,n
        testplace=i
        if (flagto(i)) then
          nn=0 
          do j=1,4
            if(type(neigh(i,j))/=0) nn=nn+1
          end do
          ebit1=dble(0)
          call poti(ebit1,i)
           do j=1,3
           if (per(j)/=dble(0)) then
            type(i)=j
            do k=1,3
              bl(k)=dble(-5)
              tpos(k)=dble(0)
              ul(k)=dble(5)
            end do
            ibound=0
            ifail=1
            ebit2=dble(0)
            call e04jaf(tn,ibound,bl,ul,tpos,ebit2,iw,liw,w,lw,ifail)
           call calb(ebitb,i,j)
            dpos(i,1)=dble(0)
            dpos(i,2)=dble(0)
            dpos(i,3)=dble(0)
            testtime=-log(g05caf(rnd))/&
                      (dble(1.0E+014)*per(j)*dexp(-(ebit2-ebit1-ebitb-&
                      dble(3)*kt/dble(2))/kt))
            call best(testtime,i,j,tpos(1),tpos(2),tpos(3),besttime,bestplace,&
                      bestto,bpos,0,bestplace2)
            end if
           end do
          type(i)=0
        end if
        if (flagfrom(i)) then
          nn=0
          do j=1,4
            if(type(neigh(i,j))/=0) nn=nn+1
          end do
          call calb(ebitb,i,0)
          ebit1=dble(0)
          call poti(ebit1,i)
          j=type(i)
          type(i)=0
          ebit2=dble(0)
          call pot(ebit2,i)
          type(i)=j
          tpos(1)=dble(0)
          tpos(2)=dble(0)
          tpos(3)=dble(0)
          freqcal=dble(0)
          do l=1,4
            if (neigh(i,l)/=0) then
             if (type(neigh(i,l))/=0) then
               freqcal=freqcal+freq(j,type(neigh(i,l)))
             end if
            end if
          end do
          freqcal=freqcal/dble(nn)
          testtime=-log(g05caf(rnd))/&
                    (dexp(-((ebit2+ebitb-ebit1)/kt))*freqcal)
           call best(testtime,i,0,dble(0),dble(0),dble(0),besttime,bestplace,&
                    bestto,bpos,0,bestplace2)
        do j=1,n
           if(flag2(i,j).and.flagto(j))&
           call diffuse(i,j,ebit1,ebitb,besttime,bestplace,bestplace2,&
                        bpos,bestto,freqcal)
         end do
        end if
      end do
      if (bestto<4) then
       type(bestplace)=bestto
       if (bestto>0) atoms=atoms+1
       if (bestto==0) atoms=atoms-1
       dpos(bestplace,1)=bpos(1)
       dpos(bestplace,2)=bpos(2)
       dpos(bestplace,3)=bpos(3)
       write(100,*) bestplace,bestto
       write(100,*) bpos(1),bpos(2),bpos(3)
       write(101,*) atoms
       if (bestto==0) then
         flagto(bestplace)=.true.
         flagfrom(bestplace)=.false.
         do i=1,4
           if(type(neigh(bestplace,i))/=0) flagfrom(neigh(bestplace,i))=.true.
         end do
       else
         flagto(bestplace)=.false.
         flagfrom(bestplace)=.true.
         do i=1,4
           if(neigh(bestplace,i)/=0)then
             if(type(neigh(bestplace,i))==0) flagto(neigh(bestplace,i))=.true.
             nn=0
             do j=1,4
             if (type(neigh(neigh(bestplace,i),j))/=0) nn=nn+1
             end do
             if (nn==4) flagto(neigh(bestplace,i))=.false.
           end if
         end do
       end if
      else
       type(bestplace2)=type(bestplace)
       type(bestplace)=0
       flagto(bestplace)=.true.
       flagfrom(bestplace)=.false.
       do i=1,4
         if(type(neigh(bestplace,i))/=0) flagfrom(neigh(bestplace,i))=.true.
       end do
       do k=1,3
         dpos(bestplace,k)=dble(0)
         dpos(bestplace2,k)=bpos(k)
       end do
       write(100,*) bestplace,0
       write(100,*) dble(0),dble(0),dble(0)
       write(100,*) bestplace2,type(bestplace2)
       write(100,*) bpos(1),bpos(2),bpos(3)
       write(101,*) atoms
       flagto(bestplace2)=.false.
       flagfrom(bestplace2)=.true.
       do i=1,4
         if(neigh(bestplace2,i)/=0)then
           if(type(neigh(bestplace2,i))==0) flagto(neigh(bestplace2,i))=.true.
           nn=0
           do j=1,4
             if (type(neigh(neigh(bestplace2,i),j))/=0) nn=nn+1
           end do
           if (nn==4) flagto(neigh(bestplace2,i))=.false.
         end if
       end do
      end if
      return
      end subroutine move

      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(tas(j)==0.or.tas(k)==0) then
         else
            if (j/=k)  o(j,k)=angle(as(j),a,as(k))
         end if
      end do
      end do
      do i=1,4
           if(tas(i)==0) then
             r(i)=dble(0)
           else
             r(i)=length(a,as(i),.false.)-r0(ta,tas(i))
           end if
      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(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.and.neigh(as(j),i)/=0) then
        ot=angle(a,as(j),neigh(as(j),i))
        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
!		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.and.&
                        neigh(as(k),l)/=0) then
         ot=angle(a,as(k),neigh(as(k),l))
         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 r3(a,b,c)
      use commvar
      implicit none
      integer a,b,c,ta,tb,tc
      ta=type(a)
      tb=type(b)
      tc=type(c)
      if(ta==0.or.tb==0.or.tc==0) then
        r3=dble(0)
      else
        r3=(r0(ta,tb)+r0(tb,tc))/dble(2)
      end if
      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)
      if(ta==0.or.tb==0.or.tc==0.or.td==0) then
        r41=dble(0)
      else
        r41=(r0(ta,tb)+r0(ta,tc)+r0(ta,td))/dble(3)
      end if
      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 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 best(tt,p,to,tp1,tp2,tp3,bt,bp,bto,bps,t2,bp2)
      use commvar
      implicit none
      integer p,to,bp,t2,bp2,bto
      real*8 :: tt,bt,tp1,tp2,tp3,bps(3)
      if (tt<bt)then
        bt=tt
        bp=p
        bto=to
        bps(1)=tp1
        bps(2)=tp2
        bps(3)=tp3
        bp2=t2
      end if
      return
      end subroutine best

      subroutine diffuse(from,to,ebit1,ebitb,besttime,bestplace,&
                         bestplace2,bpos,bestto,freqcal)
      use commvar
      implicit none
      integer,parameter ::tn=3,liw=5,lw=39
      integer :: i,j,k,l,m,ibound,ifail,from,to
      integer :: iw(liw),bestplace,bestto,bestplace2
      real*8 :: bpos(3),tpos(3),besttime,testtime
      real*8 :: bl(3),ul(3),w(lw),g05caf,rnd
      real*8 :: ebit1,ebit2,ebitb,ebitb2,length,freqcal
      type(to)=type(from)
      type(from)=0
      do k=1,3
        bl(k)=dble(-5)
        tpos(k)=dble(0)
        ul(k)=dble(5)
      end do
      ibound=0
      ifail=1
      ebit2=dble(0)
      testplace=to
      call e04jaf(tn,ibound,bl,ul,tpos,ebit2,iw,liw,w,lw,ifail)
      dpos(to,1)=tpos(1)
      dpos(to,2)=tpos(2)
      dpos(to,3)=tpos(3)
      call calb(ebitb2,to,0)
      dpos(to,1)=dble(0)
      dpos(to,2)=dble(0)
      dpos(to,3)=dble(0)
      testplace=from
      type(from)=type(to)
      testtime=-log(g05caf(rnd))/(freqcal*dexp(-(ebit2+ebitb-ebit1-ebitb2)/kt))
      type(to)=0
      call best(testtime,from,4,tpos(1),tpos(2),tpos(3),besttime,bestplace,&
                    bestto,bpos,to,bestplace2)
      return
      end subroutine diffuse