M.E.A.M Python Script

From EVOCD
Jump to: navigation, search

< Back to Interatomic Potential Generation

Purpose

This python script is designed to accept FCC Aluminum properties obtained in DFT calculations or from experimentation and use them to run the Modified Embedded Atom Method.

Implementation

  • The default values for FCC Aluminum have been included in the "Default Settings" section of the code.
  • Replace the cohesive energy Ecoh, equilibrium lattice constant a, and bulk modulus Ebulk.
#!/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()

Personal tools
Namespaces

Variants
Actions
home
Materials
Material Models
Design
Resources
Projects
Education
Toolbox