c n.b H(x) is the Heaviside step function c Pair potentials V(r) V_Fe-Fe = (9.7342365892908E+03/r) # *( 0.18180*exp(-2.8616724320005E+01*r) # +0.50990*exp(-8.4267310396064E+00*r) # +0.28020*exp(-3.6030244464156E+00*r) # +0.02817*exp(-1.8028536321603E+00*r) # )*H(1.0000-r) # +exp( 7.4122709384068E+00 # -6.4180690713367E-01*r # -2.6043547961722E+00*r**2 # +6.2625393931230E-01*r**3 # )*H(r-1.0000)*H(2.0500-r) # +(-2.7444805994228E+01*(H(2.2-r)*(2.2-r)**3) # +1.5738054058489E+01*(H(2.3-r)*(2.3-r)**3) # +2.2077118733936E+00*(H(2.4-r)*(2.4-r)**3) # -2.4989799053251E+00*(H(2.5-r)*(2.5-r)**3) # +4.2099676494795E+00*(H(2.6-r)*(2.6-r)**3) # -7.7361294129713E-01*(H(2.7-r)*(2.7-r)**3) # +8.0656414937789E-01*(H(2.8-r)*(2.8-r)**3) # -2.3194358924605E+00*(H(3.0-r)*(3.0-r)**3) # +2.6577406128280E+00*(H(3.3-r)*(3.3-r)**3) # -1.0260416933564E+00*(H(3.7-r)*(3.7-r)**3) # +3.5018615891957E-01*(H(4.2-r)*(4.2-r)**3) # -5.8531821042271E-02*(H(4.7-r)*(4.7-r)**3) # -3.0458824556234E-03*(H(5.3-r)*(5.3-r)**3) # )*H(r-2.0500) V_Fe-P = (5.6159057245908E+03/r) # *( 0.1818*exp(-2.6329090970098E+01*r) # +0.5099*exp(-7.7530945066009E+00*r) # +0.2802*exp(-3.3149971099538E+00*r) # +0.02817*exp(-1.6587327311161E+00*r) # )*H(1.0000-r) # +exp( 1.0761854424880E+01 # -1.0004045788895E+01*r # +4.9854254472397E+00*r**2 # -1.2599788569372E+00*r**3 # )*H(r-1.0000)*H(2.0000-r) # +(-3.3136605743629E+00*(H(5.3-r)*(5.3-r)**4) # +1.2625238193602E+01*(H(5.3-r)*(5.3-r)**5) # -2.0361693308072E+01*(H(5.3-r)*(5.3-r)**6) # +1.7629292543942E+01*(H(5.3-r)*(5.3-r)**7) # -8.8120728047659E+00*(H(5.3-r)*(5.3-r)**8) # +2.5494288609989E+00*(H(5.3-r)*(5.3-r)**9) # -3.9698390783403E-01*(H(5.3-r)*(5.3-r)**10) # +2.5779015833433E-02*(H(5.3-r)*(5.3-r)**11) # )*H(r-2.0000) V_P-P = (3.2399456103409E+03/r) # *( 0.18180*exp(-2.3822786399080E+01*r) # +0.50990*exp(-7.0150661324540E+00*r) # +0.28020*exp(-2.9994377000591E+00*r) # +0.02817*exp(-1.5008355431420E+00*r) # )*H(0.9000-r) # +exp( 9.9382842499617E+00 # -8.5637164272526E+00*r # +3.4519627285990E+00*r**2 # -6.1453831350215E-01*r**3 # )*H(r-0.9000)*H(2.5000-r) # +(-7.8293794709143E-02*(H(5.3-r)*(5.3-r)**4) # +3.7557214911646E-02*(H(5.3-r)*(5.3-r)**5) # )*H(r-2.5000) c c Density functions: c phi_FeFe = 1.1686859407970E+01*(H(2.4-r)*(2.4-r)**3) # -1.4710740098830E-02*(H(3.2-r)*(3.2-r)**3) # +4.7193527075943E-01*(H(4.2-r)*(4.2-r)**3) phi_FeP = phi_FeFe*(21./24.)**2 phi_PP = phi_FeFe*(21./24.)**4 rho is defined as the sum of phi_XY for a given atom. c c Embedding energy (F(rho) c F_Fe= -rho**0.5 # -6.7314115586063E-04*(rho**2) # +7.6514905604792E-08*(rho**4) F_P = -rho**0.5+1.1950274540243E-03*(rho**2) Total energy on an atom is then given by 1/2 Sum_V(r) + F(rho) For checking purposes, here are the results obtained with this potential U = -4.01299E+00 eV/atom U1 = -4.01299E+00 eV/atom H = -4.01299E+00 eV/atom Current values of ro: romin = 2.631E+01; romax = 2.631E+01 The nearest neighbor distance = 2.47 A Average number of neighbors: 60; Energy of the system without vacancy: -4.01298231 eV/atom Unrelaxed formation energy of the vacancy of the type 1: 1.8369 eV/atom Elastic constants: 2nd derivatives of the pairwise functions: -(1-1)'': -2.7444805994228E+01*(-H(2.2-r)*6*(2.2-r))+1.5738054058489E+01*(-H(2.3-r)*6*(2.3-r))+2.2077118733936E+00*(-H(2.4-r)*6*(2.4-r))-2.4989799053251E+00*(-H(2.5-r)*6*(2.5-r))+4.2099676494795E+00*(-H(2.6-r)*6*(2.6-r))-7.7361294129713E-01*(-H(2.7-r)*6*(2.7-r))+8.0656414937789E-01*(-H(2.8-r)*6*(2.8-r))-2.3194358924605E+00*(-H(3.0-r)*6*(3.0-r))+2.6577406128280E+00*(-H(3.3-r)*6*(3.3-r))-1.0260416933564E+00*(-H(3.7-r)*6*(3.7-r))+3.5018615891957E-01*(-H(4.2-r)*6*(4.2-r))-5.8531821042271E-02*(-H(4.7-r)*6*(4.7-r))-3.0458824556234E-03*(-H(5.3-r)*6*(5.3-r)) 2nd derivatives of the density functions: -(1-1)'': 1.1686859407970E+01*(-H(2.4-r)*6*(2.4-r))-1.4710740098830E-02*(-H(3.2-r)*6*(3.2-r))+4.7193527075943E-01*(-H(4.2-r)*6*(4.2-r)) 2nd derivatives of the embedding energy: (1)'': 1.0000000000000E+00*(0.5*0.5/ro^1.5)-6.7314115586063E-04*(2)+7.6514905604792E-08*(12*ro^2) C11=243.33 GPa=1.519 eV/A^3 C12=145.05 GPa=0.905 eV/A^3 C44=116.15 GPa=0.725 eV/A^3 C33=243.33 GPa=1.519 eV/A^3 C13=145.05 GPa=0.905 eV/A^3 C66=116.15 GPa=0.725 eV/A^3 B=177.81 GPa=1.110 eV/A^3