#!/usr/bin/env python from math import sqrt from fractions import gcd from sys import argv import math import os import sys import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt #------------------------------------------------------------------------------- # # define function to read parameters from TE directory # #------------------------------------------------------------------------------- def read_params(x): x.modelname = raw_input('Model name has been read\n') x.latticetype = raw_input('Lattice type has been read\n') x.species = raw_input('Species has been read\n') x.a_0 = raw_input('Initial guess of lattice parameter has been read\n') x.ingamma = raw_input('Shear period has been read\n') x.svector = raw_input('Shear direction has been read\n') x.nvector = raw_input('Shear plane normal has been read\n') x.flag = raw_input('Relaxation flag has been read\n') #------------------------------------------------------------------------------- # # Define a function that return the G.C.D for more than two numbers # Need import argv from sys # #----------------------------------------------------------------------------- def gcdm(*args): return reduce(gcd,args) #------------------------------------------------------------------------------ # # define the function to get input file # template: file name of template input file # cbkiminput: cbkim input file name # #----------------------------------------------------------------------------------------- def get_input_file(x): replacements = {'sed_latticematrix':x.latticematrix, 'sed_nbasis':str(x.nbasis), 'sed_basispos':x.basispos, 'sed_mass':str(x.massarray), 'sed_model_string':str(x.modelname), 'sed_F':x.DefG} TDdir = os.path.dirname(os.path.realpath(__file__)) # Get the name of the dir that contains the TD file_input = open(TDdir+'/cbkim.in.template','r') file_write = open('output/cbkim.in','w') for line in file_input: if 'sed_species_string' in line: s=line s11=s.replace('sed_species_string' , x.species) s21=s11.replace('sed_type_string',x.latticetype) file_write.write(s21) elif 'sed_key_string' in line: s=line s11=s.replace('sed_key_string',x.flag) s21=s11.replace('sed_step_number',str(x.stepnumber)) file_write.write(s21) else: for src, target in replacements.iteritems(): line = line.replace(src, target) file_write.write(line) file_input.close() file_write.close() #-------------------------------------------------------------------------------------- # # Define the function to write result file # #-------------------------------------------------------------------------------------- def get_output_file(x): TDdir = os.path.dirname(os.path.realpath(__file__)) file_output = open('output/cbkim.out','r') f = open(TDdir+'/result.edn.template','r') f_result = open('output/results.edn','w') x.sd=[] x.sisd=[] x.gammalis=[] for line in file_output: if 'energy density' in line: s=line s_list=s.split() energy=float(s_list[-1])*1.602176474*math.pow(10,-19) x.sd.append(float(s_list[-1])) x.sisd.append(energy) if 'Relaxed lattice vectors' in line: s=line s=next(file_output) s_list=s.split() a1=s_list[-1] a2=s_list[-3] if x.latticetype in ('diamond','fcc'): lat=float(a1)*2 if x.latticetype in ('sc','b2','bcc'): lat=float(a2) latm=lat*math.pow(10,10) latm=str(latm) lat=str(lat) p=0 y=1 for y in range(1,x.stepnumber+2): p=float(x.ingamma)/float(x.stepnumber)*(y-1) x.gammalis.append(p) y=y+1 for line in f: if '_latticetype_' in line: s=line f_result.write(s.replace('_latticetype_', x.latticetype)) elif '_source-value_' in line: s=line f_result.write(s.replace('_source-value_',lat)) elif '_species_' in line: s=line f_result.write(s.replace('_species_',x.species)) elif '_letter_' in line: s=line f_result.write(s.replace('_letter_',x.wyckoffmultletter)) elif '_spacegroup_' in line: s=line f_result.write(s.replace('_spacegroup_',x.space_group)) elif '_shear-plane-normal_' in line: s=line f_result.write(s.replace('_shear-plane-normal_',x.nvector)) elif '_shear-direction_' in line: s=line f_result.write(s.replace('_shear-direction_',x.svector)) elif '_shear-parameter_' in line: s=line f_result.write(s.translate(None,'_shear-parameter_')) y=0 for y in range(0,x.stepnumber+1): s=' {} \n'.format(x.gammalis[y]) y=y+1 f_result.write(s) elif '_sourceenergy_' in line: s=line f_result.write(s.translate(None,'_sourceenergy')) y=0 for y in range(0,x.stepnumber+1): s=' {} \n'.format(x.sd[y]) y=y+1 f_result.write(s) elif '_basisatomcoor_' in line: s=line f_result.write(s.replace('_basisatomcoor_',x.basisatomcoords)) elif '_wyckoffcoor_' in line: s=line f_result.write(s.replace('_wyckoffcoor_',x.wyckoffcoords)) elif '_specieslist_' in line: s=line f_result.write(s.replace('_specieslist_',x.specieslist)) else: s=line f_result.write(s) #-------------------------------------------------------------------------------------- # # define a function to plot the energy figure, and save the image to output directory # #--------------------------------------------------------------------------------------- def plot_energy(x): plt.figure(1) plt.plot(x.gammalis,x.sd, 'b.-') plt.xlabel('Shear parameter') plt.ylabel('Potentional energy density') plt.savefig('output/energy.png') #-------------------------------------------------------------------------------------- # # Define a function to set some imformation about crystical according to the input # #-------------------------------------------------------------------------------------- def set_lattice(x): # get every component of shear direction and shear plane normal s_list=x.svector.split() s1=float(s_list[-3]) s2=float(s_list[-2]) s3=float(s_list[-1]) s_list=x.nvector.split() n1=float(s_list[-3]) n2=float(s_list[-2]) n3=float(s_list[-1]) dot=n1*s1+n2*s2+n3*s3 if dot != 0: sys.exit('Error: shear direction and shear plane normal is not perpendicular') # Divide s1 s2 s3 and n1 n2 n3 by their G.C.D (if s=[0 2 4], change it to [0 1 2]) args=[s1,s2,s3] dividerS=gcdm(*args) s1=s1/dividerS s2=s2/dividerS s3=s3/dividerS args=[n1, n2, n3] dividerN=gcdm(*args) n1=n1/dividerN n2=n2/dividerN n3=n3/dividerN #--------------------------------------------------------------------------------------- # Calculate shear parameter to have one period crystal invariance shear for given s and v. # Take diamond(nonrelax) as an example to explain how to identify the gamma to get one period lattice invariance shear. # # Assume that s=[a b c], v=[d e f], shear parameter is y # # So, F=[ a*d*y + 1, a*e*y, a*f*y] # [ b*d*y, b*e*y + 1, b*f*y] # [ c*d*y, c*e*y, c*f*y + 1] # Express F in the fractional c.s. we get # Ftilte = [ (a*d*y)/2 + (a*e*y)/2 + (b*d*y)/2 + (b*e*y)/2 - (c*d*y)/2 - (c*e*y)/2 + 1, (a*d*y)/2 + (b*d*y)/2 + (a*f*y)/2 - (c*d*y)/2 + (b*f*y)/2 - (c*f*y)/2, (a*e*y)/2 + (a*f*y)/2 + (b*e*y)/2 + (b*f*y)/2 - (c*e*y)/2 - (c*f*y)/2] #[ (a*d*y)/2 + (a*e*y)/2 - (b*d*y)/2 - (b*e*y)/2 + (c*d*y)/2 + (c*e*y)/2, (a*d*y)/2 - (b*d*y)/2 + (a*f*y)/2 + (c*d*y)/2 - (b*f*y)/2 + (c*f*y)/2 + 1, (a*e*y)/2 + (a*f*y)/2 - (b*e*y)/2 - (b*f*y)/2 + (c*e*y)/2 + (c*f*y)/2] #[ (b*d*y)/2 - (a*e*y)/2 - (a*d*y)/2 + (b*e*y)/2 + (c*d*y)/2 + (c*e*y)/2, (b*d*y)/2 - (a*d*y)/2 - (a*f*y)/2 + (c*d*y)/2 + (b*f*y)/2 + (c*f*y)/2, (b*e*y)/2 - (a*f*y)/2 - (a*e*y)/2 + (b*f*y)/2 + (c*e*y)/2 + (c*f*y)/2 + 1] # If there's shuffle relaxation, the components of Ftilte is integer is good enough to guarentee the lattice reproduce. Because the shuffle will find its energy minimization position, which is the lattice reproduce position, by itself. # However, if there's no shuffle relaxation, the shuffle should also go back to it's position. # The shuffle vector of diamond is P1=[0.25 0.25 0.25] # The shuffle vector after deformation is P=F*P1 should go to [0.25+Z, 0.25+Z, 0.25+Z], where Z is integer. # P= # (a*d*y)/4 + (a*e*y)/4 + (a*f*y)/4 + 1/4 # (b*d*y)/4 + (b*e*y)/4 + (b*f*y)/4 + 1/4 # (c*d*y)/4 + (c*e*y)/4 + (c*f*y)/4 + 1/4 # # In summary, to get the y that make one period lattice invariance shear for diamond(nonrelax) lattice, we need to get the smallest y that satisfy the following requirments: # 1. a*y*(d+e+f), b*y*(d+e+f), c*y*(d+e+f) should be the multiple of 4 # 2. (d+e)*(a+b-c)*y, (d+f)*(a+b-c)*y, (e+f)*(a+b-c)*y, (d+e)*(a+c-b)*y, (d+f)*(a+c-b)*y, (e+f)*(a+c-b)*y, (d+e)*(c+b-a)*y, (d+f)*(c+b-a)*y, (e+f)*(c+b-a)*y should be the multiple of 2 # # As the components of P are some polynomials over 4, the components of Ftilter are some polynomials over 2, which means the y that guarantee the components of Ftilter are integers might not guarantee the components of P are integers, but the y that guarantee the components of P are integers can also guarantee the components of Ftilter are integers. # So, for diamond nonrelax, we first need to make sure the components of P are integers. How to do that? # 1. We get the greatest common divisor (gcd) of [s1*(n1+n2+n3)/4, s2*(n1+n2+n3)/4, s3*(n1+n2+n3)/4], which we call it A. A times y must be the smallest integer, namely, 1. # 2. To get y, we just need to use 1 divide A. # 3. But, sometimes [s1*(n1+n2+n3)/4, s2*(n1+n2+n3)/4, s3*(n1+n2+n3)/4] can be all zeros. In this case, we need to consider the components of Ftilter. # 4. We get the (gcd) of every component of Ftilter, which is [(n1+n2)*(s1+s3-s2)/2,(n1+n3)*(s1+s3-s2)/2,(n2+n3)*(s1+s3-s2)/2,(n1+n2)*(s1+s2-s3)/2,(n1+n3)*(s1+s2-s3)/2,(n2+n3)*(s1+s2-s3)/2,(n1+n2)*(s2+s3-s1)/2,(n1+n3)*(s2+s3-s1)/2,(n2+n3)*(s2+s3-s1)/2]. We call it B. Again, B times y must be the smallest integer, which is 1. # 5. To get y, we just need to use 1 divide B. if x.latticetype in ('bcc'): args=[(s1*n1+s1*n3/2), (s1*n2+s1*n3/2), s1*n3/2, (s2*n1+s2*n3/2), (s2*n2+s2*n3/2), s2*n3/2, (2*s3*n1-s2*n1-s1*n3/2-s1*n1-s2*n3/2+s3*n3), (2*s3*n2-s1*n3/2-s2*n2-s2*n3/2-s1*n2+s3*n3), (s3*n3-s2*n3/2-s1*n3/2)] tgamma=1/gcdm(*args) if x.latticetype in ('diamond'): if x.flag in ('nrlx'): args=[s1*(n1+n2+n3)/4,s2*(n1+n2+n3)/4,s3*(n1+n2+n3)/4] args2=[(n1+n2)*(s1+s3-s2)/2,(n1+n3)*(s1+s3-s2)/2,(n2+n3)*(s1+s3-s2)/2,(n1+n2)*(s1+s2-s3)/2,(n1+n3)*(s1+s2-s3)/2,(n2+n3)*(s1+s2-s3)/2,(n1+n2)*(s2+s3-s1)/2,(n1+n3)*(s2+s3-s1)/2,(n2+n3)*(s2+s3-s1)/2] if abs(gcdm(*args))!=0: tgamma=1/abs(gcdm(*args)) elif abs(gcdm(*args))==0: tgamma=1/abs(gcdm(*args2)) elif x.flag in ('relx'): args = [(n1+n2)*(s1+s3-s2)/2,(n1+n3)*(s1+s3-s2)/2,(n2+n3)*(s1+s3-s2)/2,(n1+n2)*(s1+s2-s3)/2,(n1+n3)*(s1+s2-s3)/2,(n2+n3)*(s1+s2-s3)/2,(n1+n2)*(s2+s3-s1)/2,(n1+n3)*(s2+s3-s1)/2,(n2+n3)*(s2+s3-s1)/2] tgamma=1/abs(gcdm(*args)) if x.latticetype in ('fcc'): args = [(n1+n2)*(s1+s3-s2)/2,(n1+n3)*(s1+s3-s2)/2,(n2+n3)*(s1+s3-s2)/2,(n1+n2)*(s1+s2-s3)/2,(n1+n3)*(s1+s2-s3)/2,(n2+n3)*(s1+s2-s3)/2,(n1+n2)*(s2+s3-s1)/2,(n1+n3)*(s2+s3-s1)/2,(n2+n3)*(s2+s3-s1)/2] tgamma=1/abs(gcdm(*args)) if x.latticetype in ('b2'): if x.flag in ('relx'): tgamma=1 if x.flag in ('nrlx'): args=[s1*(n1+n2+n3)/2,s2*(n1+n2+n3)/2,s3*(n1+n2+n3)/2] tgamma=1/abs(gcdm(*args)) #-------------------------------------------------------------------------------------------------- #Calculate step number based on ingamma, which is the shear period in the input file x.stepnumber=int(x.ingamma)*100 #Calculate deformation gradient DefG from s and v gamma=tgamma*int(x.ingamma) F11=1+gamma*s1*n1 F12=gamma*s1*n2 F13=gamma*s1*n3 F21=gamma*s2*n1 F22=gamma*s2*n2+1 F23=gamma*s2*n3 F31=gamma*s3*n1 F32=gamma*s3*n2 F33=gamma*s3*n3+1 x.DefG= '{} {} {}\n{} {} {}\n{} {} {}\n'.format(F11,F12,F13,F21,F22,F23,F31,F32,F33) #-------------------------------------------------------------------------------------------------- # Get the element number of species if x.species in ['Al']: x.elementnum=13 elif x.species in ['Cu']: x.elementnum=29 elif x.species in ['Si']: x.elementnum=14 elif x.species in ['Ar']: x.elementnum=18 elif x.species in ['Ni']: x.elementnum=28 else: print('not-supported species') #-------------------------------------------------------------------------------------------------- # Set lattice information base on different lattice type. if x.latticetype in ['bcc']: x.latticetype = 'bcc' x.space_group = 'Im-3m' x.wyckoffmultletter= '2a' x.wyckoffcoords = '[0 0 0]' x.wyckoffspecies = x.species x.basisatomcoords = '[0 0 0] \n [0.5 0.5 0.5]' x.specieslist = '"{}" \n "{}"'.format(x.species,x.species) x.nbasis = '1' x.latticematrix = '{} 0 0\n0 {} 0\n{} {} {}\n'.format(x.a_0, x.a_0, float(x.a_0)/2, float(x.a_0)/2, float(x.a_0)/2) x.massarray = '1' x.basispos = '{} 0 0 0'.format(x.elementnum) elif x.latticetype in ['fcc']: x.latticetype = 'fcc' x.space_group = 'Fm-3m' x.wyckoffcoords = '[0 0 0]' x.wyckoffmultletter = '4a' x.wyckoffspecies = x.species x.basisatomcoords = '[0 0 0] \n [0 0.5 0.5] \n [0.5 0 0.5] \n [0.5 0.5 0] \n' x.specieslist = '"{}" \n "{}" \n "{}" \n "{}"'.format(x.species,x.species,x.species,x.species) x.nbasis = '1' x.latticematrix = '0 {} {}\n {} 0 {}\n {} {} 0\n'.format(x.a_0,x.a_0,x.a_0,x.a_0,x.a_0,x.a_0) x.massarray = '1' x.basispos = '{} 0 0 0 '.format(x.elementnum) elif x.latticetype in ['diamond']: x.latticetype = 'diamond' x.space_group = 'Fd-3m' x.wyckoffmultletter = '8a' x.wyckoffcoords = '[0 0 0]' x.wyckoffspecies = x.species x.basisatomcoords = '[0 0 0] \n [0 0.5 0.5] \n [0.5 0.5 0] \n [0.5 0 0.5] \n [0.75 0.25 0.75] \n [0.25 0.25 0.25] \n [0.25 0.75 0.75] \n [0.75 0.75 0.25]' x.specieslist = '"{}" \n "{}" \n "{}" \n "{}" \n "{}" \n "{}" \n "{}"'.format(x.species,x.species,x.species,x.species,x.species,x.species,x.species) x.nbasis = '2' x.latticematrix = '0 {} {} \n {} 0 {} \n {} {} 0'.format(float(x.a_0)/2, float(x.a_0)/2, float(x.a_0)/2, float(x.a_0)/2, float(x.a_0)/2, float(x.a_0)/2) x.massarray = '1\n 1' x.basispos = '{} 0 0 0\n {} {} {} {}'.format(x.elementnum, x.elementnum, float(x.a_0)/4, float(x.a_0)/4, float(x.a_0)/4) elif x.latticetype in ['sc']: x.latticetype = 'sc' x.spece_group = 'Pm-3m' x.wyckoffmultleter = '1a' x.wyckoffcoords = '[0 0 0]' x.wyckoffspecies= x.species x.basisatomcoords='[0 0 0]' x.specieslist = '{}\n'.format(x.species) x.nbasis = '1' x.latticematrix = '{} 0 0\n0 {} 0\n0 0 {}\n'.format(x.a_0,x.a_0,x.a_0) x.massarray='1' x.basispos = '{} 0 0 0'.format(x.elementnum) #-------------------------------------------------------------------------------------- # # Define a function to get opt.in from input.in (Don't inculde the deformation part) # #-------------------------------------------------------------------------------------- def get_opt_file(x): src=open('output/cbkim.in','r') opt=open('output/opt.in','w') for line in src: s=line opt.write(s) if 'opti' in line: li=next(src) opt.write(li) opt.write('end\nstop') break src.close() opt.close() #-------------------------------------------------------------------------------------- # # Main function # #-------------------------------------------------------------------------------------- class Input: pass x = Input() # Read input read_params(x) set_lattice(x) # Replace the string in cbkim.in.template, and generate input file cbkim.in get_input_file(x) get_opt_file(x) # Write the opt.in input to test the lattice after optimization. # If the lattice changes after optimization, return an error. cmd1='cb-kim.x < output/opt.in > output/opt.out' os.system(cmd1) comp=open('output/opt.out','r') for line in comp: if 'Relaxed lattice vectors' in line: l1=next(comp) l1_list=l1.split() l2=next(comp) l2_list=l2.split() l3=next(comp) l3_list=l3.split() if x.latticetype in ('fcc','diamond'): error=(math.fabs(float(l1_list[-1])-float(l1_list[-2]))+math.fabs(float(l2_list[-1])-float(l2_list[-3]))+math.fabs(float(l3_list[-2])-float(l3_list[-3])))/3 if x.latticetype in ('b2','sc'): error=(math.fabs(float(l1_list[-3])-float(l2_list[-2]))+math.fabs(float(l2_list[-2])-float(l3_list[-1]))+math.fabs(float(l3_list[-1])-float(l1_list[-3])))/3 if x.latticetype in ('bcc'): error=(math.fabs(float(l1_list[-1])-float(l2_list[-1]))+math.fabs(float(l1_list[-1])-float(l3_list[-1]))+math.fabs(float(l2_list[-1])-float(l3_list[-1]))+math.fabs(float(l1_list[-3])-float(l2_list[-2])))/4 if error > 0.0000001: sys.exit('Error message: the lattice changes after optimazation') # If the lattice doesn't change after optimization, run the cbkim.in cmd='cb-kim.x < output/cbkim.in > output/cbkim.out' os.system(cmd) # Change the result.edn.template to results.edn get_output_file(x) # Plot the strain energy density vs. shear parameter plot_energy(x)