      double precision function amces(p,alatt)
      implicit double precision (A-H,O-Z)
c this program simulates total energy of a Al interface
c using the embeding potential
c
      parameter (n=160)
      parameter (rad=5.5d0,nat=4320)
      dimension x(n),y(n),z(n),xn(n),yn(n),zn(n)
      dimension pos(nat,3),upos(500)
      dimension list(n,100)
      dimension xp(3),yp(3),zp(3),np(8),p(n,3)
      dimension pax(3000),edx(3000),emx(3000),pay(3000)
      dimension emy(3000),edy(3000) 
      dimension pa2(3000),em2(3000),ed2(3000)
c
      do i=1,n
         x(i)=p(i,1)
         y(i)=p(i,2)
         z(i)=p(i,3)
      end do
      xd=4.0d0*alat
      xdn=-xd
      yd=2.d0*alat
      ydn=-yd
      zd=5.0d0*alat
      zdn=-zd
      xp(1)=0.d0
      yp(1)=0.d0
      zp(1)=0.d0
      xp(2)=xdn
      yp(2)=ydn
      zp(2)=zdn
      xp(3)=xd
      yp(3)=yd
      zp(3)=zd   
      call daread(pax,pay,pa2,edx,edy,ed2,emx,emy,em2)
      call multy(n,x,y,z,pos,xp,yp,zp,0.0d0)
      call neigh(n,nat,pos,list,rad)
      fe1=fe(n,pos,pax,pay,pa2,edx,edy,ed2,emx,emy,em2,list)  
      write(*,*) fe1
      ili=0
      ac=0.0d0
      beta=20.0d0
      do imc=1,25000
      ili=ili+1
      if (ili.EQ.1000) then
         ratio=ac/1000.0d0
	 ili=0
         write(*,*) imc,fe2,ratio,step
c         write(10,*) imc,temp,fe2
         step=0.85d0*step
c         call neigh(n,nat,pos,list,rad)
	 ac=0.0d0
      end if
      call RNECUY(upos,3*n)
       do i=1,n
	 xn(i)=x(i)+step*(1.0d0-2.0d0*upos(i))
	 yn(i)=y(i)+step*(1.0d0-2.0d0*upos(n+i))
	 zn(i)=z(i)+step*(1.0d0-2.0d0*upos(2*n+i))
       end do
c      if (imc.LE.100000) then
c         do in=1,8
c            xn(np(in))=x(np(in))
c         end do
c      end if
      call multy(n,xn,yn,zn,pos,xp,yp,zp,sh)
      fe2=fe(n,pos,pax,pay,pa2,edx,edy,ed2,emx,emy,em2,list) 
      if (fe2.LT.fe1) then
         fe1=fe2
         ac=ac+1.0d0
         do i=1,n
          x(i)=xn(i)
          y(i)=yn(i)
          z(i)=zn(i)
         end do
      end if
      end do  ! imc
 300  continue
      write(*,*) shr,fe1,fem
  11  format(3(F12.7,3X))
  12  format(I3,3X,3(F12.7,3X))
      amces=fe1
      return
      end
c
      subroutine splint(XA,YA,Y2A,N,X,Y) 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XA(N),YA(N),Y2A(N)
        KLO=1
        KHI=n
  1   IF (KHI-KLO.GT.1) THEN
        K=0.5*(KHI+KLO)
        IF(XA(K).GT.X)THEN
          KHI=K
        ELSE
          KLO=K
        ENDIF
      GOTO 1
      ENDIF
      H=XA(KHI)-XA(KLO)
      HI=1.d0/H
      A=(XA(KHI)-X)*HI
      B=(X-XA(KLO))*HI
      Y=A*YA(KLO)+B*YA(KHI)+
     *      ((A*A*A-A)*Y2A(KLO)+
     *      (B*B*B-B)*Y2A(KHI))*(H*H)*0.16666667D0
      RETURN
      END
c ============= random generator ==================
      SUBROUTINE RNECUY(U,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(N)
      PARAMETER (M1=2147483563,M1MIN1=2147483562,IA1=40014,
     +           IQ1=53668,IR1=12211)
      PARAMETER (M2=2147483399,IA2=40692,IQ2=52774,IR2=2791)
      PARAMETER (XM1INV=4.656613D-10)
      DATA IX1,IX2 /123456, 654321/
      DO 10 I=1,N
C Produce integer random number X1 from first MLCG
        K=IX1/IQ1
        IX1=IA1*(IX1-K*IQ1)-K*IR1
        IF (IX1.LT.0) IX1=IX1+M1
C Produce integer random number X2 from second MLCG
        K=IX2/IQ2
        IX2=IA2*(IX2-K*IQ2)-K*IR2
        IF (IX2.LT.0) IX2=IX2+M2
C Combine
        IZ=IX1-IX2
        IF(IZ.LT.1) IZ=IZ+M1MIN1
C Normalize and transform to floating point
        U(I)=DBLE(IZ)*XM1INV
   10 CONTINUE
      RETURN
C Input of seed
      ENTRY RNE2IN(ISEED1,ISEED2)
      IX1=ISEED1
      IX2=ISEED2
      RETURN
C Output of seed
      ENTRY RNE2OT(ISEED1,ISEED2)
      ISEED1=IX1
      ISEED2=IX2
      RETURN
      END
c ============================================                                
