c ========================== PROGRAM BULK ==============================
c    
c    Program for validation and comparison of EAM potentials. 
c    It works with FCC, BCC and HCP structure
c
c    Output data: lattice parameter, bulk modulus, cohesive energy,
c    non-relaxed free surface energy {001} and vacancy formation energy,
c    energy differences between structures.
c    Only for FCC and BCC: elastic coefficients c11,c12,c44, tetragonal
c    shear constant C'.
c    Transformation of FCC structure to BCC structure.
c
c ======================================================================
      PROGRAM BULK_ATAVAHA
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*80 PAIR,DEN,EMBED,INPUT
      CHARACTER*80 EMBED_IN,PAIR_IN,DEN_IN
      CHARACTER*10 POT_ID
      CHARACTER*2 ELEMENT_ID
      PARAMETER (NAT=256,NAT27=NAT*27,NP=40)
      DIMENSION X(NP),Y(NP),SIGMA(NP),A(3),Y1(NP),Y2(NP)
      DIMENSION P(0:NAT,3),P27(NAT27,3),P227(NAT27,3)
c
      PAIR='out.pair'
      DEN='out.den'
      EMBED='out.embed'
      PI=DATAN(1.0d0)*4.0d0 !pi
      E=1.60217733d-19 !elementary charge
c--------------------- Reading input file ------------------------------
      OPEN (1,FILE='bulk.dat')
      READ(1,*) ELEMENT_ID
      READ(1,*)SAL_BCC,SAL_FCC,SAL_HCP
      READ(1,*) POT_ID
      IF (POT_ID.EQ.'dynamo') THEN 
         READ (1,*) INPUT
         CALL T1 (INPUT,RAD,NPD)
      ELSE IF (POT_ID.EQ.'tabulated') THEN
         READ (1,*) EMBED_IN
         READ (1,*) PAIR_IN
         READ (1,*) DEN_IN
         READ (1,*) NPD
         CALL T2 (EMBED_IN,PAIR_IN,DEN_IN,RAD,NPD)
      ELSE IF (POT_ID.EQ.'bulk') THEN
         READ (1,*) EMBED
         READ (1,*) PAIR
         READ (1,*) DEN
         READ (1,*) NPD,RAD
      ELSE
        WRITE (*,*) 'Please check potential format in input file.'
        STOP
      END IF
      CLOSE (1)
      CALL daread (NPD,PAIR,DEN,EMBED)!reads EAM potential
c ---------------------- testing boundary condition---------------------
      IF (4.0d0*SAL_BCC.LT.RAD) THEN
         WRITE (*,*) 'Periodical boundary condition not met.'
         STOP
      END IF
      IF (4.0d0*SAL_FCC.LT.RAD) THEN
         WRITE (*,*) 'Periodical boundary condition not met.'
         STOP
      END IF      
      IF (4.0d0*SAL_HCP.LT.RAD) THEN
         WRITE (*,*) 'Periodical boundary condition not met.'
         STOP
      END IF
c******************************** bcc **********************************
c
      ALAT=SAL_BCC ! start lattice parameter
      OPEN (UNIT=1,FILE='bcc-i.dat')!a_bcc(A),E_bcc(eV)
      DO I=1,NP
	 CALL LATT_BCC(ALAT,NAT,NAT27,P27,P) !lattice
         X(I)=ALAT
	 Y(I)=AMCB(NAT,NAT27,RAD,P27,NPD) !energy
         WRITE (1,*) X(I),Y(I)/REAL(NAT)
	 SIGMA(I)=0.0d0
         ALAT=ALAT+1.0d-3
      END DO
      CLOSE(1)
      CALL POLFIT(X,Y,SIGMA,NP,3,0,A,CHI) !quadratic approximation
      XMIN=-A(2)/(2.0d0*A(3))
c
c----------- lattice parameter, bulk modulus, cohesive energy ----------
c
      NI=0 !refined calculation of lattice parameter
      OPEN (UNIT=1,FILE='bcc-r.dat')!a_bcc(A),E_bcc(eV)
      DO AA=0.998d0,1.002d0,0.0002d0
         NI=NI+1
         CALL LATT_BCC(XMIN*AA,NAT,NAT27,P27,P) !lattice
         X(NI)=XMIN*AA
         Y(NI)=AMCB(NAT,NAT27,RAD,P27,NPD) !energy
         WRITE (1,*) X(NI),Y(NI)/REAL(NAT)
	 SIGMA(NI)=0.0d0
      END DO
      CLOSE (1)
      CALL POLFIT(X,Y,SIGMA,NI,3,0,A,CHI) !quadratic approximation
      A_BCC=-A(2)/(2.0d0*A(3))!lattice parameter (A)
      VOL0=REAL(NAT)/2.0d0*A_BCC*A_BCC*A_BCC*1.0d-30
      DO I=1,NI
    	 X(I)=1.d-10*X(I) !calculation A->m
    	 X(I)=REAL(NAT)/2.0d0*X(I)*X(I)*X(I) !volume (m^3)
         Y(I)=E*Y(I) !calculation eV->J
      END DO
      CALL POLFIT(X,Y,SIGMA,NI,3,0,A,CHI) !quadratic approximation
      B_BCC=VOL0*2.d0*A(3)/1.0d9 !bulk modulus (GPa)
      CALL LATT_BCC(A_BCC,NAT,NAT27,P27,P)
      E_BCC=AMCB(NAT,NAT27,RAD,P27,NPD)/REAL(NAT)!minimal energy (eV)
c
c --------------------------- vac. form. energy ------------------------
c
      P(200,3)=P(200,3)*1.0d3
      CALL multy(NAT,NAT27,P,P27)
      EP=AMCB(NAT,NAT27,RAD,P27,NPD)
      VFE_BCC=EP-E_BCC*REAL(NAT-1)
c
c ******************************* fcc **********************************
c
      ALAT=SAL_FCC ! start lattice parameter
      OPEN (UNIT=1,FILE='fcc-i.dat')!a_fcc(A),E_fcc(eV)
      DO I=1,NP
	     CALL LATT_FCC(ALAT,NAT,NAT27,P27,P) !lattice
         X(I)=ALAT
	     Y(I)=AMCB(NAT,NAT27,RAD,P27,NPD) !energy
         WRITE (1,*) X(I),Y(I)/REAL(NAT)
	     SIGMA(I)=0.0d0
         ALAT=ALAT+1.0d-3
      END DO
      CLOSE(1)
      CALL POLFIT(X,Y,SIGMA,NP,3,0,A,CHI) !quadratic approximation
      XMIN=-A(2)/(2.0d0*A(3))
c
c----------- lattice parameter, bulk modulus, cohesive energy ----------
c
      NI=0 !refined calculation of lattice parameter
      OPEN (UNIT=1,FILE='fcc-r.dat')!a_fcc(A),E_fcc(eV)
      DO AA=0.998d0,1.002d0,0.0002d0
         NI=NI+1
         CALL LATT_FCC(XMIN*AA,NAT,NAT27,P27,P) !lattice
         X(NI)=XMIN*AA
         Y(NI)=AMCB(NAT,NAT27,RAD,P27,NPD) !energy
         WRITE (1,*) X(NI),Y(NI)/REAL(NAT)
	     SIGMA(NI)=0.0d0
      END DO
      CLOSE (1)
      CALL POLFIT(X,Y,SIGMA,NI,3,0,A,CHI) !quadratic approximation
      A_FCC=-A(2)/(2.0d0*A(3))!lattice parameter (A)
      VOL0=REAL(NAT)/4.0d0*A_FCC*A_FCC*A_FCC*1.0d-30
      DO I=1,NI
    	 X(I)=1.d-10*X(I) !calculation A->m
    	 X(I)=REAL(NAT)/4.0d0*X(I)*X(I)*X(I) !volume (m^3)
         Y(I)=E*Y(I) !calculation eV->J
      END DO
      CALL POLFIT(X,Y,SIGMA,NI,3,0,A,CHI) !quadratic approximation
      B_FCC=VOL0*2.d0*A(3)/1.0d9 !bulk modulus (GPa)
      CALL LATT_FCC(A_FCC,NAT,NAT27,P27,P)
      E_FCC=AMCB(NAT,NAT27,RAD,P27,NPD)/REAL(NAT)!minimal energy (eV)
c
c --------------------------- vac. form. energy ------------------------
c
      P(200,3)=P(200,3)*1.0d3
      CALL multy(NAT,NAT27,P,P27)
      EP=AMCB(NAT,NAT27,RAD,P27,NPD)
      VFE_FCC=EP-E_FCC*REAL(NAT-1)
c
c ******************************* hcp **********************************
c
      ALAT=SAL_HCP ! start lattice parameter
      OPEN (UNIT=1,FILE='hcp-i.dat')!a_hcp(A),E_hcp(eV)
      DO I=1,NP
	     CALL LATT_HCP(ALAT,NAT,NAT27,P27,P) !lattice
         X(I)=ALAT
	     Y(I)=AMCB(NAT,NAT27,RAD,P27,NPD) !energy
         WRITE (1,*) X(I),Y(I)/REAL(NAT)
	     SIGMA(I)=0.0d0
         ALAT=ALAT+1.0d-3
      END DO
      CLOSE(1)
      CALL POLFIT(X,Y,SIGMA,NP,3,0,A,CHI) !quadratic approximation
      XMIN=-A(2)/(2.0d0*A(3))
c
c----------- lattice parameter, bulk modulus, cohesive energy ----------
c
      NI=0 !refined calculation of lattice parameter
      OPEN (UNIT=1,FILE='hcp-r.dat')!a_hcp(A),E_hcp(eV)
      DO AA=0.998d0,1.002d0,0.0002d0
         NI=NI+1
         CALL LATT_HCP(XMIN*AA,NAT,NAT27,P27,P) !lattice
         X(NI)=XMIN*AA
         Y(NI)=AMCB(NAT,NAT27,RAD,P27,NPD) !energy
         WRITE (1,*) X(NI),Y(NI)/REAL(NAT)
	     SIGMA(NI)=0.0d0
      END DO
      CLOSE (1)
      CALL POLFIT(X,Y,SIGMA,NI,3,0,A,CHI) !quadratic approximation
      A_HCP=-A(2)/(2.0d0*A(3))!lattice parameter (A)
      VOLC=A_HCP*A_HCP*A_HCP*1.633d0*DSIN(PI/3.0d0) ! volume of elementary cell in A^3
      VOL0=REAL(NAT)/2.0d0*VOLC*1.0d-30 !sedi
      DO I=1,NI
    	 X(I)=1.d-10*X(I) !calculation A->m
    	 X(I)=REAL(NAT)/2.0d0*X(I)*X(I)*X(I)*1.633d0*DSIN(PI/3.0d0)!volume (m^3)
         Y(I)=E*Y(I) !calculation eV->J
      END DO
      CALL POLFIT(X,Y,SIGMA,NI,3,0,A,CHI) !quadratic approximation
      B_HCP=VOL0*2.d0*A(3)/1.0d9 !bulk modulus (GPa)
      CALL LATT_HCP(A_HCP,NAT,NAT27,P27,P)
      E_HCP=AMCB(NAT,NAT27,RAD,P27,NPD)/REAL(NAT)!minimal energy (eV)
c
c --------------------------- vac. form. energy ------------------------
c
      P(200,3)=P(200,3)*1.0d3
      CALL multy(NAT,NAT27,P,P27)
      EP=AMCB(NAT,NAT27,RAD,P27,NPD)
      VFE_HCP=EP-E_HCP*REAL(NAT-1)
c
c ******************* transformation FCC to BCC ************************
c
      OPEN(1,FILE='fcctobcc.dat')
      XD=((A_BCC*2.0d0)/(A_FCC*DSQRT(2.0d0))-1.0d0)/20.0d0
      ZD=(A_BCC/A_FCC-1.0d0)/20.0d0
      XX=1.0d0-15.0d0*XD
      ZZ=1.0d0-15.0d0*ZD
c
      DO I=1,40
         CALL LATT_FCC(A_FCC,NAT,NAT27,P27,P)
         CALL FCCtoBCC(NAT27,P27,XX,XX,ZZ) !yy=xx
         EN=AMCB(NAT,NAT27,RAD,P27,NPD)/REAL(NAT)
         WRITE(1,*)I,ZZ*A_FCC,EN
         XX=XX+XD
         ZZ=ZZ+ZD
      END DO
      CLOSE(1)
c
c ************************* Elastic constants BCC **********************
c
c ------------------------------ BCC  ----------------------------------
      NI=0
      OPEN(1,FILE='c11bcc.dat')
      OPEN(2,FILE='c12bcc.dat')
      OPEN(3,FILE='c44bcc.dat')
      DO EPS=-0.002d0,0.002d0,0.0002d0 !deformation
         NI=NI+1
         X(NI)=EPS
         CALL LATT_BCC(A_BCC,NAT,NAT27,P27,P)
c
         CALL C11(NAT27,P27,1.0d0+EPS,P227)! deformation for c11
         Y(NI)=AMCB(NAT,NAT27,RAD,P227,NPD)
         CALL C12(NAT27,P27,1.0d0+EPS,P227)! deformation for c12
         Y1(NI)=AMCB(NAT,NAT27,RAD,P227,NPD)     
         CALL C44(NAT27,P27,EPS,P227) ! deformation for c44
         Y2(NI)=AMCB(NAT,NAT27,RAD,P227,NPD)
c   
         SIGMA(NI)=0.0d0
         WRITE (1,*) X(NI),Y(NI)
         WRITE (2,*) X(NI),Y1(NI)
         WRITE (3,*) X(NI),Y2(NI)
         Y(NI)=Y(NI)*E !energy J->eV
         Y1(NI)=Y1(NI)*E !energy J->eV
         Y2(NI)=Y2(NI)*E !energy J->eV
      END DO
      CLOSE(1)
      CLOSE(2)
      CLOSE(3)
      VOL0=REAL(NAT)/2.0d0*A_BCC*A_BCC*A_BCC*1.0d-30 !volume m^3
      CALL POLFIT(X,Y,SIGMA,NI,3,0,A,CHI)!quadratic approximation
      C11_BCC=2.d0*A(3)/VOL0/1.0d9 !c11 (GPa)
      CALL POLFIT(X,Y1,SIGMA,NI,3,0,A,CHI)!quadratic approximation
      C12_BCC=(A(3)/VOL0/1.0d9)-C11_BCC !c12 (GPa)
      C_BCC=(C11_BCC-C12_BCC)/2.0d0 !C' (GPa)
      CALL POLFIT(X,Y2,SIGMA,NI,3,0,A,CHI)!quadratic approximation
      C44_BCC=2.d0*A(3)/VOL0/1.0d9 !c44 (GPa)
c ------------------------------ FCC  ----------------------------------
      NI=0
      OPEN(1,FILE='c11fcc.dat')
      OPEN(2,FILE='c12fcc.dat')
      OPEN(3,FILE='c44fcc.dat')
      DO EPS=-0.002d0,0.002d0,0.0002d0 !deformation
         NI=NI+1
         X(NI)=EPS
         CALL LATT_FCC(A_FCC,NAT,NAT27,P27,P)
c
         CALL C11(NAT27,P27,1.0d0+EPS,P227)! deformation for c11
         Y(NI)=AMCB(NAT,NAT27,RAD,P227,NPD)
         CALL C12(NAT27,P27,1.0d0+EPS,P227)! deformation for c12 
         Y1(NI)=AMCB(NAT,NAT27,RAD,P227,NPD)     
         CALL C44(NAT27,P27,EPS,P227) ! deformation for c44
         Y2(NI)=AMCB(NAT,NAT27,RAD,P227,NPD)
c   
         SIGMA(NI)=0.0d0
         WRITE (1,*) X(NI),Y(NI)
         WRITE (2,*) X(NI),Y1(NI)
         WRITE (3,*) X(NI),Y2(NI)
         Y(NI)=Y(NI)*E !energy J->eV
         Y1(NI)=Y1(NI)*E !energy J->eV
         Y2(NI)=Y2(NI)*E !energy J->eV
      END DO
      CLOSE(1)
      CLOSE(2)
      CLOSE(3)
      VOL0=REAL(NAT)/4.0d0*A_FCC*A_FCC*A_FCC*1.0d-30 !volume m^3
      CALL POLFIT(X,Y,SIGMA,NI,3,0,A,CHI)!quadratic approximation
      C11_FCC=2.d0*A(3)/VOL0/1.0d9 !c11 (GPa)
      CALL POLFIT(X,Y1,SIGMA,NI,3,0,A,CHI)!quadratic approximation
      C12_FCC=(A(3)/VOL0/1.0d9)-C11_FCC !c12 (GPa)
      C_FCC=(C11_FCC-C12_FCC)/2.0d0 !C' (GPa)
      CALL POLFIT(X,Y2,SIGMA,NI,3,0,A,CHI)!quadratic approximation
      C44_FCC=2.d0*A(3)/VOL0/1.0d9 !c44 (GPa)
c
c ********************* free surface energy ****************************
c
c     -------------- BCC -----------------
      CALL LATT_BCC(A_BCC,NAT,NAT27,P27,P)
      ES1=AMCB(NAT,NAT27,RAD,P27,NPD)*E ! energy of ideal crystal
      CALL SURFACE (NAT,NAT27,P,P27,RAD) ! makes free surface in (001) plane
      ES2=AMCB(NAT,NAT27,RAD,P27,NPD)*E
      ES_BCC=(ES2-ES1)/(2.0d0*16.0d0*A_BCC*A_BCC*1.0d-20)!(J/m^2)
c     -------------- FCC -----------------
      CALL LATT_FCC(A_FCC,NAT,NAT27,P27,P)
      ES1=AMCB(NAT,NAT27,RAD,P27,NPD)*E ! energy of ideal crystal
      CALL SURFACE (NAT,NAT27,P,P27,RAD) ! makes free surface in (001) plane
      ES2=AMCB(NAT,NAT27,RAD,P27,NPD)*E
      ES_FCC=(ES2-ES1)/(2.0d0*16.0d0*A_FCC*A_FCC*1.0d-20) !(J/m^2)
c     -------------- HCP -----------------
      CALL LATT_HCP(A_HCP,NAT,NAT27,P27,P)
      ES1=AMCB(NAT,NAT27,RAD,P27,NPD)*E ! energy of ideal crystal
      CALL SURFACE (NAT,NAT27,P,P27,RAD) ! makes free surface in (001) plane
      ES2=AMCB(NAT,NAT27,RAD,P27,NPD)*E
      ES_HCP=(ES2-ES1)/(2*32.0d0*A_HCP*A_HCP*
     &                 DSIN(PI/3.0d0)*1.0d-20) !(J/m^2) 
c
c ---------------------- results -> desktop ----------------------------
      WRITE(*,*) ELEMENT_ID
      WRITE(*,*)''
      WRITE(*,10)'COMPUTED PARAMETER','BCC','FCC','HCP','UNIT'
      WRITE(*,20)'lattice parameter',A_BCC,A_FCC,A_HCP,'A'
      WRITE(*,20)'cohesive energy',E_BCC,E_FCC,E_HCP,'eV'
      WRITE(*,20)'bulk modulus',B_BCC,B_FCC,B_HCP,'GPa'
      WRITE(*,20)'non.rel.vac.form.en.',VFE_BCC,VFE_FCC,VFE_HCP,'eV'
      WRITE(*,20)'free surf. en. (100)',ES_BCC,ES_FCC,ES_HCP,'J/m^2'
      WRITE(*,*)''
      WRITE(*,30)'delta Ec [eV]'
      WRITE(*,35)'BCC-FCC',E_BCC-E_FCC
      WRITE(*,35)'BCC-HCP',E_BCC-E_HCP
      WRITE(*,35)'FCC-HCP',E_FCC-E_HCP
      WRITE(*,*)''
      WRITE(*,15)'ELASTIC PARAMETER','BCC','FCC','UNIT'
      WRITE(*,25)'c11',C11_BCC,C11_FCC,'GPa'
      WRITE(*,25)'c12',C12_BCC,C12_FCC,'GPa'
      WRITE(*,25)'c44',C44_BCC,C44_FCC,'GPa'
      WRITE(*,25)'C''',C_BCC,C_FCC,'GPa'
      WRITE(*,*)''
c ------------------------- results -> results.dat ---------------------
      OPEN (1,FILE='results.dat')
      WRITE (1,*) ELEMENT_ID
      WRITE (1,*) '' 
      WRITE(1,10)'COMPUTED PARAMETER','BCC','FCC','HCP','UNIT'
      WRITE(1,20)'lattice parameter',A_BCC,A_FCC,A_HCP,'A'
      WRITE(1,20)'cohesive energy',E_BCC,E_FCC,E_HCP,'eV'
      WRITE(1,20)'bulk modulus',B_BCC,B_FCC,B_HCP,'GPa'
      WRITE(1,20)'non.rel.vac.form.en.',VFE_BCC,VFE_FCC,VFE_HCP,'eV'
      WRITE(1,20)'free surf. en. (100)',ES_BCC,ES_FCC,ES_HCP,'J/m^2'
      WRITE(1,*)''
      WRITE(1,30)'delta Ec [eV]'
      WRITE(1,35)'BCC-FCC',E_BCC-E_FCC
      WRITE(1,35)'BCC-HCP',E_BCC-E_HCP
      WRITE(1,35)'FCC-HCP',E_FCC-E_HCP
      WRITE(1,*)''
      WRITE(1,15)'ELASTIC PARAMETER','BCC','FCC','UNIT'
      WRITE(1,25)'c11',C11_BCC,C11_FCC,'GPa'
      WRITE(1,25)'c12',C12_BCC,C12_FCC,'GPa'
      WRITE(1,25)'c44',C44_BCC,C44_FCC,'GPa'
      WRITE(1,25)'C''',C_BCC,C_FCC,'GPa'
      WRITE(1,*)''
      CLOSE(1)
C ----------------------------------------------------------------------
   10 FORMAT (A20,A9,A10,A10,A10)
   20 FORMAT (A20,F10.3,F10.3,F10.3,A8)
   30 FORMAT (A20)
   35 FORMAT (A17,F13.3)
   15 FORMAT (A17,A8,A10,A10)
   25 FORMAT (A10,F17.3,F10.3,A8)
      END
c 
c ========================== subroutines ===============================
c
c
      subroutine daread (npd,pair,den,embed)
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine reads embedding atom potential
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (A-H,O-Z)
      CHARACTER *80 pair, den, embed
      common /potenc/
     1 pax(10000),pay(10000),pa2(10000),
     2 dex(10000),dey(10000),de2(10000),
     3 emx(10000),emy(10000),em2(10000)
      open(1,file=pair)
      do i=1,npd
         read(1,*) pax(i),pay(i),pa2(i)
      end do
      close(1)
      open(1,file=den)
      do i=1,npd
         read(1,*) dex(i),dey(i),de2(i)
      end do
      close(1)
      open(1,file=embed)
      do i=1,npd
         read(1,*) emx(i),emy(i),em2(i)
      end do
      close(1)
      end
c
c
      SUBROUTINE LATT_BCC(ALAT,NAT,NAT27,POS27,POS)
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine genarates coordinates of atoms
c     for bcc-lattice.
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION POS(0:NAT,3),POS1(2,3),POS27(NAT27,3)
      ALAT2=0.5d0*ALAT
c ------1st atom-----
      POS1(1,1)=0.0d0
      POS1(1,2)=0.0d0 
      POS1(1,3)=0.0d0
c ------2nd atom-----
      POS1(2,1)=ALAT2
      POS1(2,2)=ALAT2
      POS1(2,3)=ALAT2
c -------------------------------------------------
      N=0
      DO I=1,4 !X direction
	     DO J=1,4 !Y direction
	       DO K=1,8 !Z direction
             DO L=1,2
	       N=N+1
	       POS(N,1)=POS1(L,1)+DBLE((I-1))*ALAT
	       POS(N,2)=POS1(L,2)+DBLE((J-1))*ALAT
	       POS(N,3)=POS1(L,3)+DBLE((K-1))*ALAT
             END DO
	       END DO
         END DO
      END DO
      POS(0,1)=4.0d0*ALAT
      POS(0,2)=4.0d0*ALAT
      POS(0,3)=8.0d0*ALAT
c--------------------- multiplying --------------
      CALL multy(NAT,NAT27,POS,POS27)
      END
c
      SUBROUTINE LATT_FCC(ALAT,NAT,NAT27,POS27,POS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine genarates co-ordinates of atoms
c     for fcc-lattice.
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION POS(0:NAT,3),POS1(4,3),POS27(NAT27,3)
      ALAT2=0.5d0*ALAT
c ------1st atom-----
      POS1(1,1)=0.0d0
      POS1(1,2)=0.0d0 
      POS1(1,3)=0.0d0
c ------2nd atom-----
      POS1(2,1)=ALAT2
      POS1(2,2)=0.0d0
      POS1(2,3)=ALAT2
c ------3rd atom-----
      POS1(3,1)=ALAT2
      POS1(3,2)=ALAT2
      POS1(3,3)=0.0d0
c ------4th atom-----
      POS1(4,1)=0.0d0
      POS1(4,2)=ALAT2
      POS1(4,3)=ALAT2
c -------------------------------------------------
      N=0
      DO I=1,4 !X direction
	     DO J=1,4 !Y direction
	       DO K=1,4 !Z direction
             DO L=1,4
	       N=N+1
	       POS(N,1)=POS1(L,1)+DBLE((I-1))*ALAT
	       POS(N,2)=POS1(L,2)+DBLE((J-1))*ALAT
	       POS(N,3)=POS1(L,3)+DBLE((K-1))*ALAT
             END DO
	       END DO
         END DO
      END DO
      POS(0,1)=4.0d0*ALAT
      POS(0,2)=4.0d0*ALAT
      POS(0,3)=4.0d0*ALAT
c--------------------- multiplying --------------
      CALL multy(NAT,NAT27,POS,POS27)
      END
c
      SUBROUTINE LATT_HCP(ALAT,NAT,NAT27,POS27,POS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine genarates co-ordinates of atoms
c     for hcp-lattice.
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION POS(0:NAT,3),POS1(4,3),POS27(NAT27)
      ALAT2=0.5d0*ALAT
      PI=DATAN(1.0d0)*4.0d0
      C=1.63299d0 ! optimal HCP c/a ratio
c --------1st atom---------------------
      POS1(1,1)=0.0d0
      POS1(1,2)=0.0d0
      POS1(1,3)=0.0d0
c --------2nd atom---------------------
      POS1(2,1)=ALAT2
      POS1(2,2)=ALAT*DSIN(pi/3.0d0)
      POS1(2,3)=0.0d0
c --------3rd atom---------------------
      POS1(3,1)=ALAT2
      POS1(3,2)=ALAT*DSIN(PI/3.0d0)/3.0d0
      POS1(3,3)=C*ALAT2
c --------4th atom---------------------
      POS1(4,1)=ALAT
      POS1(4,2)=4.0d0/3.0d0*ALAT*DSIN(PI/3.0d0)
      POS1(4,3)=C*ALAT2
c -------------------------------------------------
      N=0
      DO I=1,4 !X direction
	     DO J=1,4 !Y direction
	       DO K=1,4 !Z direction
             DO L=1,4
	       N=N+1
           POS(N,1)=POS1(L,1)+DBLE(I-1)*ALAT
           POS(N,2)=POS1(L,2)+DBLE(J-1)*2.0*ALAT*DSIN(PI/3.0d0)
           POS(N,3)=POS1(L,3)+DBLE(K-1)*C*ALAT
             END DO
	       END DO
         END DO
      END DO
      POS(0,1)=4.0d0*ALAT
      POS(0,2)=4.0d0*2.d0*ALAT*DSIN(PI/3.0d0)
      POS(0,3)=4.0d0*C*ALAT
c--------------------- multiplying --------------
      CALL multy(NAT,NAT27,POS,POS27)
      END
c
      SUBROUTINE FCCtoBCC (NAT27,POS27,XX,YY,ZZ)
ccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine makes deformation from FCC
c     to BCC structure
ccccccccccccccccccccccccccccccccccccccccccccccccc
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION POS27(NAT27,3)
      DO I=1,NAT27
         POS27(I,1)=POS27(I,1)*XX
         POS27(I,2)=POS27(I,2)*YY
         POS27(I,3)=POS27(I,3)*ZZ
      END DO
      END
c
c
      double precision function dist(pos,i,j,nat27)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This function calculates distance between i and j atom.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (A-H,O-Z)
      dimension pos(nat27,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)
      end
c
      subroutine neigh(nat,nat27,pos,list,rad)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine finds neighbours for each atom
c     in a given radius (radius is specified in variable rad).
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (A-H,O-Z)
      dimension pos(nat27,3),list(nat,200),llist(200)
      do i=1,nat
        ne=0
	do j=1,nat27
	   d=dist(pos,i,j,nat27) 
           if (rad.gt.d) then
              if (i.ne.j) then
	         ne=ne+1
                 llist(ne)=j !position of neighbour
              end if
           end if
	end do
	list(i,1)=ne !number of neighbours for i-atom
	do k=1,ne !positions of neighbours of i-atom
	   list(i,k+1)=llist(k) 
	end do
      end do
      end
c
c
      double precision function fei (i,pos,list,nat,nat27,NPD)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      This function simulates energy of atom i.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (A-H,O-Z)
      common /potenc/
     1 pax(10000),pay(10000),pa2(10000),
     2 dex(10000),dey(10000),de2(10000),
     3 emx(10000),emy(10000),em2(10000)
      dimension pos(nat27,3),list(nat,200)
c
      pot=0.0d0
      rho=0.0d0
      nsus=list(i,1)
      do j=1,nsus
         k=list(i,j+1)
         dis=dist(pos,i,k,nat27)
         call splint(pax,pay,pa2,NPD,dis,y1)
         call splint(dex,dey,de2,NPD,dis,y2)
         pot=pot+y1
         rho=rho+y2
      end do
      call splint(emx,emy,em2,NPD,rho,emb1)
      fei=0.5d0*pot+emb1
      end
c 
      subroutine multy(nat,nat27,pos,pos27)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine makes supercell of NAT27 atoms.
c     It multiplies bulk of NAT atoms 3 times in x-,y-,z-direction.  
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (A-H,O-Z)
      dimension pos(0:nat,3),pos27(nat27,3)
      dimension xp(3),yp(3),zp(3)
c--------------------------------
      xp(1)=0.0d0
      yp(1)=0.0d0
      zp(1)=0.0d0
      xp(2)=-pos(0,1)
      yp(2)=-pos(0,2)
      zp(2)=-pos(0,3)
      xp(3)=pos(0,1)
      yp(3)=pos(0,2)
      zp(3)=pos(0,3) 
c--------------------------------
      j=0
      do ix=1,3
	     do iy=1,3
	       do iz=1,3
		     do i=1,nat
		        j=j+1
                        pos27(j,1)=pos(i,1)+xp(ix)
                        pos27(j,2)=pos(i,2)+yp(iy)
                        pos27(j,3)=pos(i,3)+zp(iz)
		     end do
	       end do
	     end do
      end do
      end     
c
      double precision function amcb(nat,nat27,rad,pos,NPD)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This function simulates total energy of a bulk
c     using the embedding atom potential.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (A-H,O-Z)
      dimension list(nat,200),pos(nat27,3)
      call neigh(nat,nat27,pos,list,rad)
      fat=0.0d0
      do i=1,nat
         fat=fat+fei(i,pos,list,nat,nat27,NPD)
      end do
      amcb=fat
      end
c
c
      SUBROUTINE C11(NAT27,POS27,DELTA,POS227)
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine makes deformation of lattice
c     to compute parameter c11
ccccccccccccccccccccccccccccccccccccccccccccccccccc
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION POS27(NAT27,3),POS227(NAT27,3)
       DO I=1,NAT27
          POS227(I,1)=POS27(I,1)*DELTA
          POS227(I,2)=POS27(I,2)
          POS227(I,3)=POS27(I,3)
       END DO
      END
c
c
      SUBROUTINE C12(NAT27,POS27,DELTA,POS227)
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine makes deformation of lattice
c     to compute parameter c12
ccccccccccccccccccccccccccccccccccccccccccccccccccc
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION POS27(NAT27,3),POS227(NAT27,3)
       DO I=1,NAT27
          POS227(I,1)=POS27(I,1)*DELTA
          POS227(I,2)=POS27(I,2)*DELTA
          POS227(I,3)=POS27(I,3)
       END DO
      END
c
c
      SUBROUTINE C44(NAT27,POS27,DELTA,POS227)
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine makes deformation of lattice
c     to compute parameter c44
ccccccccccccccccccccccccccccccccccccccccccccccccccc
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION POS27(NAT27,3),POS227(NAT27,3)
       PI=DATAN(1.0d0)*4.0d0
C       BETA=PI/4.0d0-DATAN(DELTA/2.0d0)
       BETA=PI/4.0d0-DELTA/2.0d0
       AA=DSQRT(2.0d0)*DCOS(BETA)
       BB=DSQRT(2.0d0)*DSIN(BETA)
       CALL ROTACIA(POS27,NAT27,0.0,0.0,45.0,POS227)
       DO I=1,NAT27
          POS227(I,1)=POS227(I,1)
          POS227(I,2)=POS227(I,2)*AA
          POS227(I,3)=POS227(I,3)*BB
       END DO
      END
c
c
      SUBROUTINE SURFACE (NAT,NAT27,POS,POS27,R)
cccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This subroutine makes free surface in crystal
c     It splits crystal in (001) direction
cccccccccccccccccccccccccccccccccccccccccccccccccccc
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION POS27(NAT27,3),POS(0:NAT,3)
       HR=0.5d0*POS(0,3)
       P=R+0.2d0 !vacuum layer
       DO I=1,NAT
          IF (POS(I,3).GT.HR) POS(I,3)=POS(I,3)+P
       END DO
       POS(0,3)=POS(0,3)+P
       CALL multy(NAT,NAT27,POS,POS27) 
      END
c
c
      SUBROUTINE ROTACIA (P,N,APSI,AFI,ATETA,PR)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Subroutine for rotational transformation.
c     P - array of coordinates of atoms
c     N - number of atoms in array P
c     APSI,AFI,ATETA - rotational angles in degrees
c     PR - array of rotated coordinates of atoms
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         PARAMETER (PI=3.141592654)
         REAL*8 PR(N,3),P(N,3)
         REAL*8 E(3,3)
         REAL*8 PSI,FI,TETA
c *********************** degrees to radians ***************************
         PSI=APSI/180*PI
         FI=AFI/180*PI
         TETA=ATETA/180*PI
c *************************** euler matrix *****************************
         E(1,1)=DCOS(PSI)*DCOS(FI)-DCOS(TETA)*DSIN(PSI)*DSIN(FI)
         E(1,2)=DCOS(FI)*DSIN(PSI)+DCOS(TETA)*DCOS(PSI)*DSIN(FI)
         E(1,3)=DSIN(FI)*DSIN(TETA)
         E(2,1)=-DSIN(FI)*DCOS(PSI)-DCOS(TETA)*DSIN(PSI)*DCOS(FI)
         E(2,2)=-DSIN(FI)*DSIN(PSI)+DCOS(TETA)*DCOS(PSI)*DCOS(FI)
         E(2,3)=DCOS(FI)*DSIN(TETA)
         E(3,1)=DSIN(TETA)*DSIN(PSI)
         E(3,2)=-DSIN(TETA)*DCOS(PSI)
         E(3,3)=DCOS(TETA)
c ********************* multiplication of matrixes *********************
         DO K=1,N
            DO I=1,3
               SUMA=0.
               DO J=1,3
                  SUMA=SUMA+P(K,J)*E(I,J)   
               END DO
               PR(K,I)=SUMA
            END DO
         END DO
      END
C
C
c ==================== downloaded subroutines =========================
c http://www.irb.hr/users/maks/v5/splint.for
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XA(N),Y2A(N),YA(N)
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) pause 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.
      return
      END
c                               
C SUBROUTINE POLFIT.F 
C
C SOURCE
C   BEVINGTON, PAGES 140-142.
C
C PURPOSE
C   MAKE A LEAST-SQUARES FIT TO DATA WITH A POLYNOMIAL CURVE
C      Y = A(1) + A(2)*X + A(3)*X**2 + A(3)*X**3 + . . .
C
C USAGE 
C   CALL POLFIT (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, CHISQR)
C
C DESCRIPTION OF PARAMETERS
C   X	   - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C   Y	   - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C   SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS
C   NPTS   - NUMBER OF PAIRS OF DATA POINTS
C   NTERMS - NUMBER OF COEFFICIENTS (DEGREE OF POLYNOMIAL + 1)
C   MODE   - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
C	     +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
C	      0 (NO WEIGHTING) WEIGHT(I) = 1.
C	     -1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
C   A	   - ARRAY OF COEFFICIENTS OF POLYNOMIAL
C   CHISQR - REDUCED CHI SQUARE FOR FIT 
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
C   DETERM (ARRAY, NORDER)
C      EVALUATES THE DETERMINANT OF A SYMETRICAL
C      TWO-DIMENSIONAL MATRIX OF ORDER NORDER
C
C COMMENTS
C   DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
C
      SUBROUTINE POLFIT (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR)
	  implicit double precision (A-H,O-Z)
c        DOUBLE PRECISION SUMX,SUMY,XTERM,YTERM,ARRAY,CHISQ
	  DIMENSION X(1),Y(1),SIGMAY(1),A(1)
	  DIMENSION SUMX(19),SUMY(10),ARRAY(10,10)
C
C ACCUMULATE WEIGHTED SUMS
C
  11	NMAX=2*NTERMS-1 
    	DO 13 N=1,NMAX
  13	SUMX(N)=0.d0
    	DO 15 J=1,NTERMS
  15	SUMY(J)=0.d0
    	CHISQ=0.d0
  21	DO 50 I=1,NPTS
    	XI=X(I) 
	    YI=Y(I) 
  31	IF (MODE) 32,37,39
  32	IF (YI) 35,37,33
  33	WEIGHT=1./YI
    	GOTO 41 
  35	WEIGHT=1.d0/(-YI) 
    	GOTO 41 
  37	WEIGHT=1.d0
    	GOTO 41 
  39	WEIGHT=1.d0/SIGMAY(I)**2
  41	XTERM=WEIGHT
    	DO 44 N=1,NMAX
	       SUMX(N)=SUMX(N)+XTERM
  44	XTERM=XTERM*XI
  45	YTERM=WEIGHT*YI 
    	DO 48 N=1,NTERMS
	       SUMY(N)=SUMY(N)+YTERM
  48	YTERM=YTERM*XI
  49	CHISQ=CHISQ+WEIGHT*YI**2
  50	CONTINUE
C
C CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS 
C
  51	DO 54 J=1,NTERMS
    	DO 54 K=1,NTERMS
	    N=J+K-1 
  54	ARRAY(J,K)=SUMX(N)
    	DELTA=DETERM (ARRAY,NTERMS)
	    IF (DELTA) 61,57,61
  57	CHISQR=0.d0
    	DO 59 J=1,NTERMS
  59	A(J)=0.d0 
    	GOTO 80 
  61	DO 70 L=1,NTERMS
  62	DO 66 J=1,NTERMS
    	DO 65 K=1,NTERMS
	    N=J+K-1 
  65	ARRAY(J,K)=SUMX(N)
  66	ARRAY(J,L)=SUMY(J)
  70	A(L)=DETERM (ARRAY,NTERMS) /DELTA
C
C CALCULATE CHI SQUARE
C
  71	DO 75 J=1,NTERMS
    	CHISQ=CHISQ-2.d0*A(J)*SUMY(J)
	    DO 75 K=1,NTERMS
	    N=J+K-1 
  75	CHISQ=CHISQ+A(J)*A(K)*SUMX(N)
  76	FREE=NPTS-NTERMS
  77	CHISQR=CHISQ/FREE
  80	RETURN
    	END
c
C FUNCTION DETERM.F
C
C SOURCE
C   BEVINGTON, PAGE 294.
C
C PURPOSE
C   CALCULATE THE DETERMINANT OF A SQUARE MATRIX
C
C USAGE 
C   DET = DETERM (ARRAY, NORDER)
C
C DESCRIPTION OF PARAMETERS
C   ARRAY  - MATRIX
C   NORDER - ORDER OF DETERMINANT (DEGREE OF MATRIX)
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
C   NONE
C
C COMMENTS
C   THIS SUBPROGRAM DESTROYS THE INPUT MATRIX ARRAY
C   DIMENSION STATEMENT VALID FOR NORDER UP TO 10
C
    	double precision FUNCTION DETERM (ARRAY,NORDER)
	    implicit double precision (A-H,O-Z)
c        DOUBLE PRECISION ARRAY,SAVE
    	DIMENSION ARRAY(10,10)
C
  10	DETERM=1.d0
  11	DO 50 K=1,NORDER
C
C INTERCHANGE COLUMNS IF DIAGONAL ELEMENT IS ZERO
C
    	IF (ARRAY(K,K)) 41,21,41
  21	DO 23 J=K,NORDER
    	IF (ARRAY(K,J)) 31,23,31
  23	CONTINUE
    	DETERM=0.d0
	    GOTO 60 
  31	DO 34 I=K,NORDER
    	SAVE=ARRAY(I,J) 
	    ARRAY(I,J)=ARRAY(I,K)
  34	ARRAY(I,K)=SAVE 
    	DETERM=-DETERM
C
C SUBTRACT ROW K FROM LOWER ROWS TO GET DIAGONAL MATRIX 
C
  41	DETERM=DETERM*ARRAY(K,K)
    	IF (K-NORDER) 43,50,50
  43	K1=K+1
    	DO 46 I=K1,NORDER
	    DO 46 J=K1,NORDER
  46	ARRAY(I,J)=ARRAY(I,J)-ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
  50	CONTINUE
  60	RETURN
    	END
c ======================================================================
c
c                     Last modified: 18.03.2013              
c                                 
c                 E. Vitkovska (atavaha@hotmail.com)
c                                &
c                  P. Ballo (peter.ballo@stuba.sk)
c
c=======================================================================
