      double precision function amcb(alat,p)
      implicit double precision (A-H,O-Z)
c this program simulates total energy of a Cu bulk
c using the embeding atom potential
      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     character*9 set(11)
c
c
      do i=1,n
    	 x(i)=p(i,1)
    	 y(i)=p(i,2)
	     z(i)=p(i,3)
      end do
c ----------------------------------
      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
c
c program for calculation of EAM energy
      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)
c     step=0.005
      amcb=fe1
      return
      end
c
      subroutine daread(pax,pay,pa2,edx,edy,ed2,emx,emy,em2)
      implicit double precision (A-H,O-Z)
      dimension pax(*),pay(*),pa2(*),emx(*),emy(*),em2(*)
      dimension edx(*),edy(*),ed2(*)
      open(1,file='../potential/cu.pair1')
      read(1,*) n
      do i=1,n
	  read(1,*) pax(i),pay(i),pa2(i)
      end do
      close(1)
      open(1,file='../potential/cu.den1')
      read(1,*) n
      do i=1,n
	 read(1,*) emx(i),emy(i),em2(i)
      end do
      close(1)
      open(1,file='../potential/cu.embed1')
      read(1,*) n
      do i=1,n
	 read(1,*) edx(i),edy(i),ed2(i)
      end do
      close(1)
      return
      end
c
c
      double precision function xpot(pax,pay,pa2,x)
      implicit double precision (A-H,O-Z)
      call splint(pax,pay,pa2,3000,x,y)
      xpot=y
      return
      end
c
      double precision function xfat(edx,edy,ed2,x)
      implicit double precision (A-H,O-Z)
      dimension edx(*),edy(*),ed2(*)
      call splint(edx,edy,ed2,3000,x,y)
      xfat=y
      return
      end
c
      double precision function dist(pos,i,j)
      implicit double precision (A-H,O-Z)
      dimension pos(4320,3)
      d=0.0d0
      do l=1,3
	    d=d+(pos(i,l)-pos(j,l))*(pos(i,l)-pos(j,l))
      end do
      dist=dsqrt(d)
      return
      end
c
      subroutine neigh(n,nat,pos,list,rad)
      implicit double precision (A-H,O-Z)
      dimension pos(160,3),list(160,100),llist(100)
      n27=n*27
      do i=1,n
      ne=0
	   do j=1,n27
	    if (i.NE.j) then
	       d=dist(pos,i,j)
	       if (rad.GT.d) then
		  ne=ne+1
		  llist(ne)=j
	       end if
	    end if
	   end do
	   list(i,1)=ne
	   do k=1,ne
	    list(i,k+1)=llist(k)
	   end do
      end do
      return
      end
c
      double precision function fe
     &            (n,pos,pax,pay,pa2,edx,edy,ed2,emx,emy,em2,list)
      implicit double precision (A-H,O-Z)
      dimension pos(4320,3),list(160,100)
      dimension pax(*),pay(*),pa2(*),edx(*),edy(*),ed2(*)
      dimension emx(*),emy(*),em2(*)
      fat=0.0d0
      do i=1,n
       fat=fat+fei(i,pos,pax,pay,pa2,edx,edy,ed2,emx,emy,em2,list)
      end do
      fe=fat
      return
      end
c
      double precision function fei
     &          (i,pos,pax,pay,pa2,edx,edy,ed2,emx,emy,em2,list)
      implicit double precision (A-H,O-Z)
      dimension pos(4320,3),list(160,100)
      dimension pax(*),pay(*),pa2(*),edx(*),edy(*),ed2(*)
      dimension  emx(*),emy(*),em2(*)
      pot=0.0d0
      rho=0.0d0
      nsus=list(i,1)
      do j=1,nsus
	    k=list(i,j+1)
	    d=dist(pos,i,k)
	    pot=pot+xpot(pax,pay,pa2,d)
	    rho=rho+xpot(emx,emy,em2,d)
      end do
      fei=0.5d0*pot+xfat(edx,edy,ed2,rho)+5.d-2
      return
      end
c ============================================
      subroutine multy(n,x,y,z,pos,xp,yp,zp,sh)
      implicit double precision (A-H,O-Z)
      dimension x(*),y(*),z(*),pos(4320,3),xp(*),yp(*),zp(*)
      j=0
      do ix=1,3
	   do iy=1,3
	    do iz=1,3
		do i=1,n
		   j=j+1
		   pos(j,1)=x(i)+xp(ix)
		   pos(j,2)=y(i)+yp(iy)
		   pos(j,3)=z(i)+zp(iz)
		end do
	    end do
	   end do
      end do
      return
      end
c
c
      SUBROUTINE MINPRB(XMP,XA,YA,XB,YB,XC,YC)
      implicit double precision (A-H,O-Z)
      PARAMETER (EPS=1.d-10)
      XBA=XB-XA
      XBC=XB-XC
      F1=XBA*(YB-YC)
      F2=XBC*(YB-YA)
      ANUM=XBA*F1-XBC*F2
      DEN=F1-F2
      DEN=SIGN(MAX(ABS(DEN),EPS),DEN)
      XMP=XB-ANUM/(2.D0*DEN)
      END
c
c
      subroutine symetry(n,x,xd)
      implicit double precision (A-H,O-Z)
      dimension x(*)
      akon=0.5d0*xd+0.1d0
      do i=1,n
	     if (x(i).GT.akon) x(i)=x(i)-xd
      end do
      return
      end
c
      function locate(pax,x)
      implicit double precision (A-H,O-Z)
      dimension pax(*)
      jl=0
      ju=3001
 10   if (ju-jl.GT.1) then
         jm=(ju+jl)/2
         if ((pax(3000).GT.pax(1)).EQV.(x.GT.pax(jm))) then
              jl=jm
         else
              ju=jm
         end if
         goto 10
      end if
      locate=jl
      return
      end     
c
      function locate1(edx,x)
      implicit double precision (A-H,O-Z)
      dimension edx(*)
      jl=0
      ju=2001
 10   if (ju-jl.GT.1) then
         jm=(ju+jl)/2
         if ((edx(2000).GT.edx(1)).EQV.(x.GT.edx(jm))) then
              jl=jm
         else
              ju=jm
         end if
         goto 10
      end if
      locate1=jl
      return
      end            
