M.E.A.M Python Script
From EVOCD
Revision as of 13:19, 14 March 2016 by Huddleston (Talk | contribs)
< 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()
