      program bulk
      implicit double precision (A-H,O-Z)
      parameter (np=40)
      dimension x(80),y(80),sigma(80),a(3),p(160,3),p1(160,3)
      dimension p3(4320,3)
      alat=3.6 ! lattice constant for copper
      xmin=3.61
      go to 100
      open(unit=4,file='copper.dat')
      do i=1,np
    	 call latt(alat,p)
         x(i)=alat
	     yp=amcb(alat,p)
	     y(i)=yp
	     write(*,20) x(i),y(i)
	     write(4,20) x(i),y(i)
	     sigma(i)=0.0d0
	     alat=alat+1.0d-3
      end do
      close(4)
      call polfit(x,y,sigma,np,3,0,a,chi)
      xmin=-a(2)/(2.0d0*a(3))
      vol0=40.0d0*xmin*xmin*xmin*1.0d-30
       do i=1,np
    	x(i)=1.d-10*x(i)
    	x(i)=40.0d0*x(i)*x(i)*x(i)
        y(i)=1.60219d-19*y(i)
       end do
      call polfit(x,y,sigma,np,3,0,a,chi)
      write(*,*) (a(j),j=1,3)
      const=vol0*2.d0*a(3)
      en0=amcb(xmin,p)
      go to 100
c coeficient C'
      i=0
      alat=xmin
      open(unit=4,file='copper-c.dat')
      do aa=0.99d0,1.01d0,0.001d0
         i=i+1
    	 call lattc(alat,aa,p)
	     yp=amcc(alat,aa,p)
         x(i)=alat*aa
	     y(i)=yp
	     write(*,20) x(i),y(i)
	     write(4,20) x(i),y(i)
	     sigma(i)=0.0d0
      end do
      close(4)
      call polfit(x,y,sigma,i,3,0,a,chi)
      xmin=-a(2)/(2.0d0*a(3))
      write(*,*) xmin,alat
      vol0=40.0d0*xmin*xmin*xmin*1.0d-30
       do i=1,np
    	x(i)=1.d-10*x(i)
    	x(i)=40.0d0*x(i)*x(i)*x(i)
        y(i)=1.60219d-19*y(i)
       end do
      call polfit(x,y,sigma,np,3,0,a,chi)
      write(*,*) (a(j),j=1,3)
      const2=vol0*2.d0*a(3)/4.d0
      c12=(3.d0*const-2.d0*const2)/3.d0
      c11=2.d0*const2+c12
c  c44
      ni=0
      open(unit=4,file='copper-c44.dat')
      do aa=0.99d0,1.01d0,0.001d0
         ni=ni+1
    	 call lattc44(alat,aa,p)
	     yp=amcc(alat,aa,p)
         x(ni)=alat*aa
	     y(ni)=yp
	     write(*,20) x(ni),y(ni)
	     write(4,20) x(ni),y(ni)
	     sigma(ni)=0.0d0
      end do
      close(4)
      call polfit(x,y,sigma,ni,3,0,a,chi)
      xmin=-a(2)/(2.0d0*a(3))
      write(*,*) xmin,alat
      vol0=40.0d0*xmin*xmin*xmin*1.0d-30
       do i=1,ni
    	x(i)=1.d-10*x(i)
    	x(i)=40.0d0*x(i)*x(i)*x(i)
        y(i)=1.60219d-19*y(i)
       end do
      call polfit(x,y,sigma,ni,3,0,a,chi)
      write(*,*) (a(j),j=1,3)
      const3=vol0*2.d0*a(3)/4.d0
c
c      en0=amcb(xmin,p)
      coef=159.d0/160.d0
      p(100,1)=p(100,1)+1.d3
      en1=amcb(xmin,p)
      en1s=0.d0
      en1v=en1-coef*en0
      en2v=en1s-coef*en0
c bcc energy
      alat=4.05d0
      open(unit=4,file='copper-bcc.dat')
      do i=1,np
    	 call latt_bcc(alat,p)
         x(i)=alat
	     yp=amcb_bcc(alat,p)
	     y(i)=yp
	     write(*,20) x(i),y(i)
	     write(4,20) x(i),y(i)
	     sigma(i)=0.0d0
	     alat=alat+1.0d-3
      end do
      close(4)
      call polfit(x,y,sigma,np,3,0,a,chi)
      xmin_bcc=-a(2)/(2.0d0*a(3))
      call latt_bcc(xmin_bcc,p)
      e_bcc=amcb_bcc(xmin_bcc,p)
      vol0=40.0d0*xmin_bcc*xmin_bcc*xmin_bcc*1.0d-30
       do i=1,np
    	x(i)=1.d-10*x(i)
    	x(i)=40.0d0*x(i)*x(i)*x(i)
        y(i)=1.60219d-19*y(i)
       end do
      call polfit(x,y,sigma,np,3,0,a,chi)
      write(*,*) (a(j),j=1,3)
      const_bcc=vol0*2.d0*a(3)
c SC energy
      alat=1.99d0
      open(unit=4,file='copper-sc.dat')
      do i=1,np
    	 call latt_bcc(alat,p)
         call rotation(p,0.d0,45.d0,0.d0)
         call latt_sc(alat,p)
        call rotation(p,0.d0,-45.d0,0.d0)
         x(i)=alat
	     yp=amcb_sc(alat,p)
	     y(i)=yp
	     write(*,20) x(i),y(i)
	     write(4,20) x(i),y(i)
	     sigma(i)=0.0d0
	     alat=alat+1.0d-3
      end do
      close(4)
      call polfit(x,y,sigma,np,3,0,a,chi)
      xmin_sc=-a(2)/(2.0d0*a(3))
      call latt_bcc(xmin_sc,p)
      call rotation(p,0.d0,45.d0,0.d0)
      call latt_sc(xmin_sc,p)
      call rotation(p,0.d0,-45.d0,0.d0)
      e_sc=amcb_sc(xmin_sc,p)
      vol0=40.0d0*xmin_sc*xmin_sc*xmin_sc*1.0d-30
       do i=1,np
    	x(i)=1.d-10*x(i)
    	x(i)=40.0d0*x(i)*x(i)*x(i)
        y(i)=1.60219d-19*y(i)
       end do
      call polfit(x,y,sigma,np,3,0,a,chi)
      write(*,*) (a(j),j=1,3)
      const_sc=vol0*2.d0*a(3)
 100  continue
      xx=0.d0
      yy=0.d0
      zz=0.d0
      call latt_small(xmin,p3,xx,yy,zz,m)
      yb=amcb_small(p3)
      write(*,*) yb/108.d0,en0/160.d0
      open(10,file='smyk.dat')
      do i=1,70
         call latt_small(xmin,p3,xx,yy,zz,m)
         yb=amcb_small(p3)
         xx=xx+0.01d0
         yy=yy+0.01d0
         zz=zz-0.02d0
         write(*,*)  xx,yb/108.d0
         write(10,*) xx,yb/108.d0
      end do
      close(10)
      sfa=4.5d0*dsqrt(3.d0)*xmin*xmin*1.d-20
      sfe=en0-
      stop  
c results    
      write(*,*) 
      write(*,*) 'minimum energy at ',xmin
      write(*,*) 'energy = ',en0,' eV'
      write(*,*) 'energy per atom = ',en0/160.0d0,'eV/at'
      write(*,*) 'optimal volume = ',vol0,' m^3'
      write(*,*) 'computed bulk modulus'
      write(*,5) const
      write(*,*) 'computed C'
      write(*,5) const2
      write(*,*) 'computed C44'
      write(*,5) const3
      write(*,*) 'computed c11,c12 coeficients'
      write(*,10) c11,c12
      write(*,*) 'non rel. vacancy formation energy = ',en1v,' eV'
c      write(*,*) 'relaxed  vacancy formation energy = ',en2v,' eV'
      write(*,*) 'bcc-lattice constant'
      write(*,5)  dsqrt(2.0d0)*xmin_bcc/2.0d0 
      write(*,*) 'E_bcc - E_fcc'
      write(*,5) (e_bcc-en0)/160.d0 
      write(*,*) 'sc-lattice constant'
      write(*,5)  xmin_sc 
      write(*,*) 'E_sc - E_fcc'
      write(*,5) (e_sc-en0)/160.d0  
   5  format (E12.6E2)
  10  format (2(E12.6E2,3X))
  20  format (2X,F12.8,3X,F14.8)
      end
c  


