View source for ICME 2012 POT-GEN-1

[[ICME_2012_HW1| < Homework 1]]

Author(s): [https://sites.google.com/site/lsiliyanage/home Laalitha Liyanage]

==License==
By using the codes provided here you accept the the [[CodeRepository:License Agreement|Mississippi State University's license agreement]]. Please read the agreement carefully before usage.

==Overview==
The script calculates cohesive energy, lattice constant, bulk modulus, equilibrium volume, elastic constants (c11,c12,c44), surface formation energies (for (100),(110), (111)) and defect energies (vacancy formation and self-interstitial energies) for fcc structure.
* This script is a Python script. Therefore save it as <pre><name.py></pre> and do <pre>chmod u+x <name>.py</pre> to make it executable. To run the script type <pre>./<name>.py</pre>

==Inputs==
To simplify execution of the script all inputs are contained in the '''Default Settings''' section. It is at the top of the script and allows you to set options for all these calculations.These include in order

* Name of file to write geometrical data of simulation box
* Beginning lattice parameter
* Parameters of interatomic potential
** pair_style will specify interatomic potential formulation type such as EAM or MEAM.
** pair_coeff will specify the file names which contain coefficients for the interatomic potential formulation. Default file names are library.meam and Al.meam. For workability of all the scripts it is suggested not to change these names.
* Unit conversion constants 
**To convert units given out by LAMMPS to experimentally comparable units.
* Command to execute LAMMPS
**By default it is assumed you have placed the LAMMPS exe in ~/bin. If it's not the case please modify this command. Also for LAMMPS to work you should have <pre>swsetup openmpi-intel-64</pre> in your ~/.bashrc file. Remember the number of processors allocated in this command, since you will have to use the same number when preparing a pbs command script.
*Percentage strain for elastic constant calculation and number of intervals at which the energy is calculated.
*Length of the vacuum for surface energy calculation

==Output==

The outputs of the script include log files and dump files of all 15 calculation and a final response file (''response.txt'') which contains the target properties with comparison to experimental or DFT data.

An example output file looks like

<pre>
E_coh           -3.410657       (-3.35) eV/atom
a_0             4.045260        (4.05)  Angstrom(1e-10m)
V_0             16.549286       (16.61) Angstrom^3/atom
C11             109.928085      (107)   GPa
C12             61.245716       (61)    GPa
C44             32.590087       (28)    GPa
B               77.473172       (82)    GPa
E_100           496.526733      (890)   mJ/m^2
E_110           581.321360      (960)   mJ/m^2
E_111           427.811200      (780)   mJ/m^2
GSFE_I          127.662730      (133)   mJ/m^2
GSFE_I          128.321278      (133)   mJ/m^2
E_vac           0.657647        (0.5)   eV
E_Oct           2.399749        (2.8)   eV
E_Tetra         2.728030        (3.3)   eV
</pre>

==The script==

<pre>
#!/usr/local/bin/python
#Authors: Laalitha Liyanage, Sungho Kim
#Purpose: Calculate material properties of fcc system

import os
import re
import sys
import math
import numpy as np
#========================================================================================================================================
#--------------------------------------------------Default setting-----------------------------------------------------------------------
#========================================================================================================================================

datafile="datafile"

# Default lattice parameter
a=4.05 

#Paramters of interatomic potential
pair_style = "meam"
pair_coeff = "* * library.meam Al Al.meam Al"

#Constants to convert units
unit_conversion1 = 160.217646 # 1 eV/A^3 = 1.60217646E-19 / 1E-30 Pa = 160.217646 GPa
unit_conversion2 = 16021.765 #1 eV/A^3 to mJ/m^2

#Command to execute LAMMPS
executable = "mpirun -np 8 ~/bin/lmp_talon.static"

#Percentage strain and no. of points for elastic constant calculation
perc = .1
nt = 6

#Length of vacuum for surface energy calculation
vacuum =  15.0

#=============================================================================================================================#
def get_fcc_coords(a,nx=1,ny=1,nz=1,vac=0):#Generate coordinates for of FCC
 
  xa = []; ya = []; za = []; ty = []
  ax = a
  ay = a
  az = a
  for i in range(nx):
    for j in range(ny):
      for k in range(nz):        
        if vac == 0:
                xa.append( (0.0+i)*ax ); ya.append( (0.0+j)*ay ); za.append( (0.0+k)*az ); ty.append(1)
                xa.append( (0.0+i)*ax ); ya.append( (0.5+j)*ay ); za.append( (0.5+k)*az ); ty.append(1)
                xa.append( (0.5+i)*ax ); ya.append( (0.0+j)*ay ); za.append( (0.5+k)*az ); ty.append(1)
                xa.append( (0.5+i)*ax ); ya.append( (0.5+j)*ay ); za.append( (0.0+k)*az ); ty.append(1)
        if vac == 1:
                if i==nx/2 and j==ny/2 and k==nz/2:
                        xa.append( (0.0+i)*ax ); ya.append( (0.5+j)*ay ); za.append( (0.5+k)*az ); ty.append(1)
                        xa.append( (0.5+i)*ax ); ya.append( (0.0+j)*ay ); za.append( (0.5+k)*az ); ty.append(1)
                        xa.append( (0.5+i)*ax ); ya.append( (0.5+j)*ay ); za.append( (0.0+k)*az ); ty.append(1)

                else:
                        xa.append( (0.0+i)*ax ); ya.append( (0.0+j)*ay ); za.append( (0.0+k)*az ); ty.append(1)
                        xa.append( (0.0+i)*ax ); ya.append( (0.5+j)*ay ); za.append( (0.5+k)*az ); ty.append(1)
                        xa.append( (0.5+i)*ax ); ya.append( (0.0+j)*ay ); za.append( (0.5+k)*az ); ty.append(1)
                        xa.append( (0.5+i)*ax ); ya.append( (0.5+j)*ay ); za.append( (0.0+k)*az ); ty.append(1)

  return ty,xa,ya,za,ax*nx,ay*ny,az*nz

def distort_simbox(xa,ya,za,box,strain_tensor=[0,0,0,0,0,0]):#Apply strain and calculate the 9 triclinic box parameters for LAMMPS

  #Initial orthorgonal box edges taken as lattice translation vectors
  ax,ay,az = box[0],box[1],box[2]

  #Strains
  e1,e2,e3,e4,e5,e6 = strain_tensor

  #Application of strain to the lattice vectors
  ax_p = [(1.0+e1)*ax,0.5*e6*ax,0.5*e5*ax]
  ay_p = [0.5*e6*ay,(1.0+e2)*ay,0.5*e4*ay]
  az_p = [0.5*e5*az,0.5*e4*az,(1.0+e3)*az]
   
  #The new lattice vectors (cell vectors)
  cell =  [ax_p,ay_p,az_p]
  
  #Get length of cell vectors
  a = np.linalg.norm(cell[0]) 
  b = np.linalg.norm(cell[1]) 
  c = np.linalg.norm(cell[2]) 
  
  #Calculate angles between cell vectors
  alpha = np.arccos( np.vdot(cell[1], cell[2]) / (b * c) )
  beta = np.arccos( np.vdot(cell[0], cell[2]) / (a * c) )
  gamma = np.arccos( np.vdot(cell[0], cell[1]) / (a * b) )
 
  # a_LAMMPS = (xhi-xlo,0,0); b_LAMMPS = (xy,yhi-ylo,0); c_LAMMPS = (xz,yz,zhi-zlo)
  xlo = ylo = zlo = 0.0
  # this choice of origin simplifies things:
  # a_LAMMPS = (xhi,0,0); b_LAMMPS = (xy,yhi,0); c_LAMMPS = (xz,yz,zhi)
  
  xhi = a # a_LAMMPS
  xy = np.cos(gamma) * b# b_LAMMPS
  yhi = np.sin(gamma) * b
  xz = np.cos(beta) * c# c_LAMMPS
  yz = ( b * c * np.cos(alpha) - xy * xz ) / yhi
  zhi = np.sqrt( c**2 - xz**2 - yz**2 )
 
  #LAMMPS vectors
  
  cell_lammps = np.array([[xhi-xlo,0,0],[xy,yhi-ylo,0],[xz,yz,zhi-zlo]])

  bx = [xhi-xlo,0,0]
  by = [xy,yhi-ylo,0]
  bz = [xz,yz,zhi-zlo]
 
  #print cell, cell_lammps
  
  # IMPORTANT: need vector-(rotation-)matrix product (instead of matrix-vector product) here,
  #            cell vectors are ROW VECTORS (also see above)
  
  rotation = np.dot(np.linalg.inv(cell), cell_lammps)

  #print rotation

  for i in range(len(xa)):
    r = [xa[i],ya[i],za[i]]
    [x,y,z] = np.dot(r,rotation)
    #print r, [x,y,z]
    xa[i]=x;ya[i]=y;za[i]=z

  return xa,ya,za,bx,by,bz

def gen_datafile(ty,xa,ya,za,bx,by,bz):
  fout = open(datafile,"w")
  fout.write("cementite a = %f, b = %f, c = %f\n\n" % (bx[0],by[1],bz[2]))
  fout.write("%d atoms\n"%len(xa))
  fout.write("1 atom types\n")
  fout.write(" 0.0  %22.16f   xlo xhi\n"%bx[0])
  fout.write(" 0.0  %22.16f   ylo yhi\n"%by[1])
  fout.write(" 0.0  %22.16f   zlo zhi\n"%bz[2])
  fout.write(" %22.16f  %22.16f %22.16f xy xz yz\n"%(by[0],bz[0],bz[1]))
  fout.write("\nAtoms\n\n")
  for i in range(len(xa)):
    fout.write("%4d %3d %22.16f %22.16f %22.16f\n"%(i+1,ty[i],xa[i],ya[i],za[i]))
  fout.close()
  return len(xa)

def gen_infile(relax=0,box_relax=0):# Create input command script for lammps
  if relax == 1:
      fname = "infile_relax"
  else:
      fname = "infile_static"

  fout = open(fname,'w')
  fout.write("units           metal\n")
  fout.write("boundary        p p p\n")
  fout.write("atom_style      atomic\n")
  fout.write("read_data       datafile\n")
  fout.write("pair_style      %s\n"%pair_style)
  fout.write("pair_coeff      %s\n"%pair_coeff)
  fout.write("neighbor        2.0 bin\n")
  fout.write("neigh_modify    delay 10\n")
  fout.write("dump            1 all custom 1 dump id type x y z\n")
  fout.write("log             log\n")
  fout.write("thermo_style    custom step atoms pe ke etotal press vol lx ly lz\n")
  fout.write("thermo_modify   line one format float %.16g\n")
  
  if relax ==0 and box_relax ==0:
      fout.write("run             0\n")
  elif relax == 1 and box_relax == 0:
      fout.write("min_style       cg\n")
      fout.write("min_modify      line quadratic\n")
      fout.write("minimize        1.0e-30 1.0e-20 100000 1000000\n")
  elif relax == 1 and box_relax ==1:
      fout.write("fix             1 all box/relax iso 0.0 vmax 0.001\n")
      fout.write("min_style       cg\n")
      fout.write("min_modify      line quadratic\n")
      fout.write("minimize        1.0e-30 1.0e-20 100000 1000000\n")
  fout.close()
  return fname

def get_field_from_log_lammps(field):# Get data from Lammps log file
  os.system("grep Loop log -B1|grep -v Loop|awk '{print $%s}' > eout"%(field))
  fin = file("eout",'r')
  line = fin.readline().split()
  fin.close()
  os.system("rm eout")
  return (float(line[0]))

def get_E0():#Relax structure and get equilibrium energy and volume.
  ty,xa,ya,za,ax,ay,az = get_fcc_coords(a)
  xa,ya,za,bx,by,bz = distort_simbox(xa,ya,za,[ax,ay,az],strain_tensor=[0,0,0,0,0,0])
  gen_datafile(ty,xa,ya,za,bx,by,bz)
  os.system("%s -in %s > report.lammps"%(executable,gen_infile(relax=1,box_relax=1)))
  N =  get_field_from_log_lammps("2")
  E0 = get_field_from_log_lammps("3")
  V0 = get_field_from_log_lammps("7")
  a0 = get_field_from_log_lammps("8")
  os.system("mv log log.E0")
  os.system("mv dump dump.E0")
  return(N,E0, V0, a0)

##---------------Function definitions relating to Elastic Constant calculation-----------------------------------##

def get_cxx(cname,a0,a,b,c):
  energy = []
  os.system("echo %s > report.lammps"%(cname))
  t0,t1 = -perc/100,perc/100
  dt = (t1 - t0)/nt

  ty,xa,ya,za,ax,ay,az = get_fcc_coords(a0)

  for i in range(nt+1):

    t = t0+i*dt

    #Generate lattice according to strain tensor
    if cname == 'c11': 
      xa,ya,za,bx,by,bz = distort_simbox(xa,ya,za,[ax,ay,az],strain_tensor=[t,0,0,0,0,0])
    elif cname == 'c12':
      xa,ya,za,bx,by,bz = distort_simbox(xa,ya,za,[ax,ay,az],strain_tensor=[t,-t,0,0,0,0])
    elif cname == 'c44':
      xa,ya,za,bx,by,bz = distort_simbox(xa,ya,za,[ax,ay,az],strain_tensor=[0,0,0,2*t,0,0])

    gen_datafile(ty,xa,ya,za,bx,by,bz)
    #Run LAMMPS and get energy
    os.system("echo 't=%f' >> report.lammps "%(t))
    os.system("%s -in %s >> report.lammps"%(executable,gen_infile(relax=1,box_relax=0)))
    energy.append(get_field_from_log_lammps("5"))

  fout = open("summary","w")
  for i in range(nt+1):
    fout.write("%22.16f %22.16f\n"%(t0+i*dt,energy[i]))
  fout.close()
  gpfile = gen_gpfile_for_parabola_fit(a*a0**3,b,c,cname)
  os.system("gnuplot %s > gnuplot.report 2>&1"%gpfile)
  a,b,c = get_fitted_param(a*a0**3,b,c,"fit.log")
  os.system("mv fit.log %s_fit.log"%cname)
  os.system("mv report.lammps %s_report.lammps"%cname)
  return a/a0**3

def gen_gpfile_for_parabola_fit(a,b,c,prefix):
#Generate input script for GNUplot for fitting
  os.system("rm -f fit.log")
  fn = "%s_fit.gp"%prefix
  fout = open(fn,"w")
  fout.write('set term gif\n')
  fout.write('set output "%s_fit.gif"\n'%prefix)
  fout.write("f(x) = a*x**2 + b*x+ c\n")
  b = 0.0
  fout.write("a = %f; b = %f; c = %f\n"%(a,b,c))
  fout.write("FIT_LIMIT = 1e-9\n")
  if b == 0.0:
    fout.write("fit f(x) 'summary' via a,c\n")
  else:
    fout.write("fit f(x) 'summary' via a,b,c\n")
  fout.write('fx_title = sprintf("(%f)*x**2+(%f)*x+(%f)",a,b,c)\n')
  fout.write("plot 'summary' t 'MEAM' with linespoints 1 5,f(x) t fx_title")
  #fout.write("plot 'summary' with linespoints 1 5,f(x)")
  fout.close()
  return fn

def get_fitted_param(a,b,c,fn):
#Extract fitted parameters from GNUplot output
  fin = file(fn,"r")
  while True:
    line = fin.readline()
    if line == "": sys.exit(123)
    if re.search("^Final set of parameters",line): break
  for i in range(2): line = fin.readline()
  for i in range(2):
    line = fin.readline().split()
    if   len(line)>0 and line[0] == "a": a = float(line[2])
    elif len(line)>0 and line[0] == "b": b = float(line[2])
    elif len(line)>0 and line[0] == "c": c = float(line[2])
  fin.close()
  return a,b,c

##---------------Function definitions relating to Surface Energy calculation-------------------------------##

#--------------------------Surface (100)---------------------------------------------------------

def gen_coords_fcc_100(a,nx=1,ny=1,nz=1):#Generate coordinates for fcc (100) surface structure
  xa=[]; ya=[]; za=[];ty = []
  bx = [0]*3;by = [0]*3; bz = [0]*3

  bx[0],by[1],bz[2] = a*nx,a*ny,a*nz+vacuum
  x,y,z = bx,by,bz

  for i in range(nx):
    for j in range(ny):
      for k in range(nz):
        xa.append(0 + i*a); ya.append(  0 + j*a); za.append(  0 + k*a); ty.append(1)
        xa.append(  0 + i*a); ya.append(a/2 + j*a); za.append(a/2 + k*a);ty.append(1)
        xa.append(a/2 + i*a); ya.append(  0 + j*a); za.append(a/2 + k*a);ty.append(1)
        xa.append(a/2 + i*a); ya.append(a/2 + j*a); za.append(  0 + k*a);ty.append(1)
  return ty,xa,ya,za,bx,by,bz

def get_surf_energy_100(a0,E0):#Calculate fcc (100) surface formation energy
  ty,xa,ya,za,bx,by,bz = gen_coords_fcc_100(a0,2,2,4)
  gen_datafile(ty,xa,ya,za,bx,by,bz)
  os.system("%s -in %s > report.lammps"%(executable,gen_infile(relax=1,box_relax=0)))
  N =  get_field_from_log_lammps("2")
  E = get_field_from_log_lammps("3")
  e100 = (E-N*E0)/2/(get_field_from_log_lammps("8")*get_field_from_log_lammps("9"))
  os.system("mv log log.surf100")
  os.system("mv dump dump.surf100")
  return(e100)

#--------------------------Surface (110)--------------------------------------------------------------

def gen_coords_fcc_110(a,nx=1,ny=1,nz=1):#Generate coordinates for fcc (110) surface structure
  xa=[]; ya=[]; za=[]; ty=[]
  bx = [0]*3;by = [0]*3; bz = [0]*3

  ax = a*math.sqrt(2)/2
  ay = a*math.sqrt(6)/2
  az = a*math.sqrt(3)
  x0 = 0.0
  x2 = math.sqrt(2)/4. * a
  y2 = math.sqrt(6)/4. * a
  y3 = math.sqrt(6)/6. * a
  y4 = math.sqrt(6)*5./12. * a
  y5 = math.sqrt(6)*2./6. * a
  y6 = math.sqrt(6)/12 * a
  z3 = math.sqrt(3)/3. * a
  z5 = math.sqrt(3)*2./3. * a
  bx[0],by[1],bz[2] = ax*nx + vacuum, ay*ny, az*nz
  for i in range(nx):
    for j in range(ny):
      for k in range(nz):
        xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(x0+k*az);ty.append(1)
        xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(x0+k*az);ty.append(1)
        xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(z3+k*az);ty.append(1)
        xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(z3+k*az);ty.append(1)
        xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(z5+k*az);ty.append(1)
        xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(z5+k*az);ty.append(1)
  
  return ty,xa,ya,za,bx,by,bz

def get_surf_energy_110(a0,E0):#Calculate fcc (110) surface formation energy
  ty,xa,ya,za,bx,by,bz = gen_coords_fcc_110(a0,4,2,2)
  gen_datafile(ty,xa,ya,za,bx,by,bz)
  os.system("%s -in %s > report.lammps"%(executable,gen_infile(relax=1,box_relax=0)))
  N =  get_field_from_log_lammps("2")
  E = get_field_from_log_lammps("3")
  e110 = (E-N*E0)/2/(get_field_from_log_lammps("9")*get_field_from_log_lammps("10"))
  os.system("mv log log.surf110")
  os.system("mv dump dump.surf110")
  return e110

#--------------------------Surface (111)-----------------------------------------------------------------------

def gen_coords_fcc_111(a,nx=1,ny=1,nz=1):
  """ Generate datafile of FCC surface: 110:x, 112:y, 111:z """
  xa=[]; ya=[]; za=[];ty =[]
  bx = [0]*3;by = [0]*3; bz = [0]*3

  ax = a*math.sqrt(2)/2
  ay = a*math.sqrt(6)/2
  az = a*math.sqrt(3)
  x0 = 0.0
  x2 = math.sqrt(2)/4 * a
  y2 = math.sqrt(6)/4 * a
  y3 = math.sqrt(6)/6 * a
  y4 = math.sqrt(6)*5/12 * a
  y5 = math.sqrt(6)*2/6 * a
  y6 = math.sqrt(6)/12 * a
  bx[0],by[1],bz[2] = ax*nx, ay*ny, az*nz+vacuum
  for i in range(nx):
    for j in range(ny):
      layer = 0
      for k in range(nz):
        xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
        xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
        xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)

  return ty,xa,ya,za,bx,by,bz

def get_surf_energy_111(a0,E0):#Calculate fcc (111) surface formation energy
  ty,xa,ya,za,bx,by,bz = gen_coords_fcc_111(a0,2,2,4)
  gen_datafile(ty,xa,ya,za,bx,by,bz)
  os.system("%s -in %s > report.lammps"%(executable,gen_infile(relax=1,box_relax=0)))
  N =  get_field_from_log_lammps("2")
  E = get_field_from_log_lammps("3")
  e111 = (E-N*E0)/2/(get_field_from_log_lammps("8")*get_field_from_log_lammps("9"))
  os.system("mv log log.surf111")
  os.system("mv dump dump.surf111")
  return e111

##---------------Function definitions relating to Generalized Stacking Fault Energy------------------------------------##

def gen_data_for_stacking_fault(a,nx=1,ny=1,type="I"):#Generate coordinates for FCC surface: 110:x, 112:y, 111:z
  
  xa=[]; ya=[]; za=[];ty = []
  bx = [0]*3;by = [0]*3; bz = [0]*3

  ax = a*math.sqrt(2)/2
  ay = a*math.sqrt(6)/2
  az = a*math.sqrt(3)
  x0 = 0.0
  x2 = math.sqrt(2)/4 * a
  y2 = math.sqrt(6)/4 * a
  y3 = math.sqrt(6)/6 * a
  y4 = math.sqrt(6)*5/12 * a
  y5 = math.sqrt(6)*2/6 * a
  y6 = math.sqrt(6)/12 * a
  if type == "I":
    nlayer = 11
    nfaults = 1
  elif type == "E":
    nlayer = 13
    nfaults = 1

  bx[0],by[1],bz[2] = ax*nx, ay*ny, az/3.0*nlayer

  for i in range(nx):
    for j in range(ny):
      layer = 0
      for k in range(2):
        xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
        xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
        xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
      if type == "I":
        xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
        xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
        xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
        xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)            ;ty.append(1) 
        xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
        xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
      elif type == "E":                                                           
        xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
        xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
        for k in range(2):
          xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
          xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
          xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
          xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
          xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)            ;ty.append(1)
          xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1;ty.append(1)
  return ty,xa,ya,za,bx,by,bz,nfaults

def get_stacking_fault(a0,E0,type='I'):
  ty,xa,ya,za,bx,by,bz,nfaults = gen_data_for_stacking_fault(a0,2,2,type)
  gen_datafile(ty,xa,ya,za,bx,by,bz)
  os.system("%s -in %s > report.lammps"%(executable,gen_infile(relax=1,box_relax=0)))
  N =  get_field_from_log_lammps("2")
  E = get_field_from_log_lammps("3")
  gsfe = (E-N*E0)/(get_field_from_log_lammps("8")*get_field_from_log_lammps("9"))
  os.system("mv log log.gsfe.%s"%type)
  os.system("mv dump dump.gsfe.%s"%type)
  return gsfe

##-------------------Calculation of vacancy formation energy-----------------------------------------------------##

def get_vac_form_energy(a0,E0):
  ty,xa,ya,za,ax,ay,az = get_fcc_coords(a,4,4,4,vac=1)
  xa,ya,za,bx,by,bz = distort_simbox(xa,ya,za,[ax,ay,az],strain_tensor=[0,0,0,0,0,0])
  gen_datafile(ty,xa,ya,za,bx,by,bz)
  os.system("%s -in %s > report.lammps"%(executable,gen_infile(relax=1,box_relax=1)))
  N =  get_field_from_log_lammps("2")
  E = get_field_from_log_lammps("3")
  e_vac = E-E0*N
  os.system("mv log log.vac")
  os.system("mv dump dump.vac")
  return(e_vac)
  
##-------------------------Calculation of interstitial energy-----------------------------------------------------##

def get_fcc_coords_with_int(a,nx=1,ny=1,nz=1,int='oct'):#Generate coordinates for of FCC with interstitial atom
 
  xa = []; ya = []; za = []; ty = []
  ax = a
  ay = a
  az = a
  for i in range(nx):
    for j in range(ny):
      for k in range(nz):
        if int == 'oct' and i==nx/2 and j==ny/2 and k==nz/2:
           xa.append( (0.5+i)*ax ); ya.append( (0.0+j)*ay ); za.append( (0.0+k)*az ); ty.append(1)  
        if int == 'tet' and i==nx/2 and j==ny/2 and k==nz/2:
           xa.append( (0.25+i)*ax ); ya.append( (0.25+j)*ay ); za.append( (0.25+k)*az ); ty.append(1)  
        xa.append( (0.0+i)*ax ); ya.append( (0.0+j)*ay ); za.append( (0.0+k)*az ); ty.append(1)
        xa.append( (0.0+i)*ax ); ya.append( (0.5+j)*ay ); za.append( (0.5+k)*az ); ty.append(1)
        xa.append( (0.5+i)*ax ); ya.append( (0.0+j)*ay ); za.append( (0.5+k)*az ); ty.append(1)
        xa.append( (0.5+i)*ax ); ya.append( (0.5+j)*ay ); za.append( (0.0+k)*az ); ty.append(1)
  return ty,xa,ya,za,ax*nx,ay*ny,az*nz

def get_int_energy(a0,E0):
#Octahedral interstitial energy calculation
  ty,xa,ya,za,ax,ay,az = get_fcc_coords_with_int(a,4,4,4,int='oct')
  xa,ya,za,bx,by,bz = distort_simbox(xa,ya,za,[ax,ay,az],strain_tensor=[0,0,0,0,0,0])
  gen_datafile(ty,xa,ya,za,bx,by,bz)
  os.system("%s -in %s > report.lammps"%(executable,gen_infile(relax=1,box_relax=1)))
  N =  get_field_from_log_lammps("2")
  E = get_field_from_log_lammps("3")
  e_oct = E-E0*N
  os.system("mv log log.oct")
  os.system("mv dump dump.oct")
#Tetrahedral interstitial energy calculation
  ty,xa,ya,za,ax,ay,az = get_fcc_coords_with_int(a,4,4,4,int='tet')
  xa,ya,za,bx,by,bz = distort_simbox(xa,ya,za,[ax,ay,az],strain_tensor=[0,0,0,0,0,0])
  gen_datafile(ty,xa,ya,za,bx,by,bz)
  os.system("%s -in %s > report.lammps"%(executable,gen_infile(relax=1,box_relax=1)))
  N =  get_field_from_log_lammps("2")
  E = get_field_from_log_lammps("3")
  e_tet = E-E0*N
  os.system("mv log log.tet")
  os.system("mv dump dump.tet")
  return(e_oct, e_tet)

#Main program

#Get equilibrium crystal structure data
N,E0,V,a0 = get_E0()
#Get elastic constants
c11 = get_cxx("c11",a0,0.5,1.0,E0)*2*unit_conversion1
k12 = get_cxx("c12",a0,0.5,1.0,E0)*2*unit_conversion1
c22=c11
c12 = (c11+c22-k12)/2
c44 = get_cxx("c44",a0,2.0,1.0,E0)/2*unit_conversion1
B=(c11+2*c12)/3
#Get surface energies
E100 = get_surf_energy_100(a0,E0/N)*unit_conversion2
E110 = get_surf_energy_110(a0,E0/N)*unit_conversion2
E111 = get_surf_energy_111(a0,E0/N)*unit_conversion2
#Get stacking fault energies
gsfe_i = get_stacking_fault(a0,E0/N,'I')*unit_conversion2
gsfe_e = get_stacking_fault(a0,E0/N,'E')*unit_conversion2
#Get defect energies
E_vac =  get_vac_form_energy(a0,E0/N)
E_oct,E_tet = get_int_energy(a0,E0/N)

#Write results to file

fout = open("responses.txt",'w')

S = "E_coh\t\t%f\t(-3.35)\teV/atom\n" % (E0/N);fout.write(S)
S = "a_0\t\t%f\t(4.05)\tAngstrom(1e-10m)\n" % (a0);fout.write(S)
S = "V_0\t\t%f\t(16.61)\tAngstrom^3/atom\n" % (V/N);fout.write(S)
S = "C11\t\t%f\t(107)\tGPa\n" % (c11);fout.write(S)
S = "C12\t\t%f\t(61)\tGPa\n" % (c12);fout.write(S)
S = "C44\t\t%f\t(28)\tGPa\n" % (c44);fout.write(S)
S = "B\t\t%f\t(82)\tGPa\n" % (B);fout.write(S)
S = "E_100\t\t%f\t(890)\tmJ/m^2\n" % (E100);fout.write(S)
S = "E_110\t\t%f\t(960)\tmJ/m^2\n" % (E110);fout.write(S)
S = "E_111\t\t%f\t(780)\tmJ/m^2\n" % (E111);fout.write(S)
S = "GSFE_I\t\t%f\t(133)\tmJ/m^2\n" % (gsfe_i);fout.write(S)
S = "GSFE_I\t\t%f\t(133)\tmJ/m^2\n" % (gsfe_e);fout.write(S)
S = "E_vac\t\t%f\t(0.5)\teV\n" % (E_vac);fout.write(S)
S = "E_Oct\t\t%f\t(2.8)\teV\n" % (E_oct);fout.write(S)
S = "E_Tetra\t\t%f\t(3.3)\teV\n" % (E_tet);fout.write(S)

fout.close()

</pre>