#!/usr/bin/env python from math import sqrt import math import os import sys #-------------------------------------------------------------------------- # # define a function to read the input # #-------------------------------------------------------------------------- def read_input(x): x.modelname=raw_input('Enter the extended KIM ID of a Model:\n') x.latticetype=raw_input('Enter the Strukturbericht code of the crystal structure:\n') x.species1=raw_input('Enter the chemical symbol of species 1:\n') x.species2=raw_input('Enter the chemical symbol of species 2:\n') x.a_0=float(raw_input('Enter an initial guess of the lattice parameter:\n')) #-------------------------------------------------------------------------- # # define function to get the element number for given element # #-------------------------------------------------------------------------- def get_element_number(x): if x.species1 in ['Al']: x.elementnum1=13 elif x.species1 in ['Cu']: x.elementnum1=29 elif x.species1 in ['Si']: x.elementnum1=14 elif x.species1 in ['Ar']: x.elementnum1=18 elif x.species1 in ['Ni']: x.elementnum1=28 elif x.species1 in ['Au']: x.elementnum1=79 elif x.species1 in ['Ti']: x.elementnum1=22 elif x.species1 in ['Cd']: x.elementnum1=48 else: print('not-supported species1') if x.species2 in ['Al']: x.elementnum2=13 elif x.species2 in ['Cu']: x.elementnum2=29 elif x.species2 in ['Si']: x.elementnum2=14 elif x.species2 in ['Ar']: x.elementnum2=18 elif x.species2 in ['Ni']: x.elementnum2=28 elif x.species2 in ['Au']: x.elementnum2=79 elif x.species2 in ['Ti']: x.elementnum2=22 elif x.species2 in ['Cd']: x.elementnum2=48 else: print('not-supported species2') #------------------------------------------------------------------------------------------------- # # define function to get lattice information base on the given lattice type. # #------------------------------------------------------------------------------------------------- def get_lattice_info(x): if x.latticetype in ['B1']: x.species = '{}-{}'.format(x.species1,x.species2) x.space_group = 'Fm-3barm' x.specieslist = '"{}" \n "{}" \n "{}" \n "{}" \n "{}" \n "{}" \n "{}" \n "{}" '.format(x.species1, x.species1, x.species1, x.species1, x.species2, x.species2, x.species2, x.species2) x.basisatomcoords = '[0 0 0] \n [0 0.5 0.5] \n [0.5 0 0.5] \n [0.5 0.5 0] \n [0.5 0 0] \n [0 0.5 0] \n [0.5 0.5 0] \n [0.5 0.5 0.5]' x.nbasis = '2' x.latticematrix = '{} 0 0\n 0 {} 0\n 0 0 {}\n'.format(x.a_0, x.a_0, x.a_0) x.massarray = '1\n 1' x.basispos = '{} 0 0 0\n {} {} {} {}'.format(x.elementnum1,x.elementnum2, x.a_0/2 ,x.a_0/2 ,x.a_0/2) x.wyckoffmultletter = '4a 4b' x.wyckoffcoords = '[0 0 0] [0.5 0.5 0.5]' x.wyckoffspecies = x.species1 + ' ' + x.species2 elif x.latticetype in ['B2']: x.species = '{}-{}'.format(x.species1,x.species2) x.space_group = 'Pm-3barm' x.specieslist = '"{}" \n "{}"'.format(x.species1,x.species2) x.basisatomcoords = '[0 0 0] \n [0 0.5 0.5] ' x.nbasis = '2' x.latticematrix = '{} 0 0\n 0 {} 0\n 0 0 {}\n'.format(x.a_0, x.a_0, x.a_0) x.massarray = '1\n 1' x.basispos = '{} 0 0 0\n {} {} {} {}'.format(x.elementnum1,x.elementnum2, x.a_0/2, x.a_0/2, x.a_0/2) x.wyckoffmultleter = '1a 1b' x.wyckoffcoords = '[0 0 0] [0.5 0.5 0.5]' x.wyckoffspecies = x.species1 + ' ' + x.species2 elif x.latticetype in ['B3']: x.species = '{}-{}'.format(x.species1, x.species2) x.space_group = 'F4bar-3m' x.specieslist = '"{}" \n "{}" \n "{}" \n "{}" \n "{}" \n "{}" \n "{}" \n "{}"'.format(x.species1, x.species1, x.species1, x.species1, x.species2, x.species2, x.species2, x.species2) x.basisatomcoords = '[0 0 0] \n [0 0.5 0.5] \n [0.5 0 0.5] \n [0.5 0.5 0] \n [0.25 0.25 0.25] \n [0.25 0.25 0.75] \n [0.25 0.75 0.25] \n [0.75 0.25 0.25]' x.nbasis = '2' x.latticematrix = '0 {} {}\n {} 0 {}\n {} {} 0\n'.format(x.a_0/2, x.a_0/2, x.a_0/2, x.a_0/2, x.a_0/2, x.a_0/2) x.massarray = '1\n 1' x.basispos = '{} 0 0 0 \n {} {} {} {}'.format(x.elementnum1, x.elementnum2, x.a_0/2, x.a_0/2, x.a_0/2) elif x.latticetype in ['L12']: x.species = '{}-{}'.format(x.species1,x.species2) x.space_group = 'Pm-3barm' x.specieslist = '"{}" \n "{}" \n "{}" \n "{}"'.format(x.species1, x.species2, x.species2, x.species2) x.basisatomcoords = '[0 0 0] \n [0 0.5 0.5] \n [0.5 0 0.5] \n [0.5 0.5 0]' x.nbasis = '4' x.latticematrix = '{} 0 0 \n 0 {} 0 \n 0 0 {}'.format(x.a_0, x.a_0, x.a_0) x.massarray = '1\n 1\n 1\n 1' x.basispos = '{} 0 0 0\n {} {} {} {}\n {} {} {} {}\n {} {} {} {}'.format(x.elementnum1, x.elementnum2, 0, x.a_0/2, x.a_0/2, x.elementnum2, x.a_0/2, 0, x.a_0/2, x.elementnum2, x.a_0/2, x.a_0/2, 0) x.wyckoffmultleter = '1a 3c' x.wyckoffcoords = '[0 0 0] [0 0.5 0.5]' x.wyckoffspecies = x.species1 + ' ' + x.species2 #------------------------------------------------------------------------------------------------------ # # define a function to replace the string in cbkim.in.template, and generate input file cbkim.in # #------------------------------------------------------------------------------------------------------ def get_input(x): 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_latticematrix' in line: s=line file_write.write(s.replace('sed_latticematrix',x.latticematrix)) elif 'sed_nbasis' in line: s=line file_write.write(s.replace('sed_nbasis',x.nbasis)) elif 'sed_basispos' in line: s=line file_write.write(s.replace('sed_basispos',x.basispos)) elif 'sed_mass' in line: s=line file_write.write(s.replace('sed_mass',x.massarray)) elif 'sed_model_string' in line: s=line file_write.write(s.replace('sed_model_string',x.modelname)) else: s=line file_write.write(s) file_input.close() file_write.close() #--------------------------------------------------------------------------------------------------------------------- # # define a function to generate input file for check equilibrium # #--------------------------------------------------------------------------------------------------------------------- def get_input_equi(): src = open('output/cbkim.in','r') equi = open('output/equi.in','w') for line in src: s = line equi.write(s) if 'equi,,' in line: li=next(src) equi.write(li) equi.write('end\nstop') break src.close() equi.close() #--------------------------------------------------------------------------------------------------------------------- # # define a function to check equilibrium of given crystial structure # #--------------------------------------------------------------------------------------------------------------------- def check_equi(x): src = open('output/equi.out','r') for line in src: if 'dwde for CheckEquilibrium' in line: for i in range(0,int(x.nbasis)): s = src.next() s_list = s.split() for j in range(1,4): if abs(float(s_list[-j]))>1e-6: sys.exit('Error : dwde is too large. The initial lattice structure is not in equilibrium config.') #--------------------------------------------------------------------------------------------------------------------- # # define a function to read the cbkim.out and generate result.edn # #--------------------------------------------------------------------------------------------------------------------- def get_output(x): TDdir = os.path.dirname(os.path.realpath(__file__)) # Get the name of the dir that contains the TD file_output = open('output/cbkim.out','r') f = open(TDdir+'/result.edn.template','r') f_result = open('output/results.edn','w') for line in file_output: if 'Relaxed lattice vectors' in line: s=file_output.next() s_list=s.split() if x.latticetype in ['B1','B3']: lat=float(s_list[-1])*math.sqrt(2) if x.latticetype in ['B2','L12']: lat=float(s_list[-3]) if 'elastic stiffness matrix' in line: s=file_output.next() s_list = s.split() c11=float(s_list[-6]) c12=float(s_list[-5]) c13=float(s_list[-4]) c14=float(s_list[-3]) c15=float(s_list[-2]) c16=float(s_list[-1]) s=file_output.next() s_list = s.split() c21=float(s_list[-6]) c22=float(s_list[-5]) c23=float(s_list[-4]) c24=float(s_list[-3]) c25=float(s_list[-2]) c26=float(s_list[-1]) s=file_output.next() s_list = s.split() c31=float(s_list[-6]) c32=float(s_list[-5]) c33=float(s_list[-4]) c34=float(s_list[-3]) c35=float(s_list[-2]) c36=float(s_list[-1]) s=file_output.next() s_list=s.split() c41=float(s_list[-6]) c42=float(s_list[-5]) c43=float(s_list[-4]) c44=float(s_list[-3]) c45=float(s_list[-2]) c46=float(s_list[-1]) s=file_output.next() s_list=s.split() c51=float(s_list[-6]) c52=float(s_list[-5]) c53=float(s_list[-4]) c54=float(s_list[-3]) c55=float(s_list[-2]) c56=float(s_list[-1]) s=file_output.next() s_list=s.split() c61=float(s_list[-6]) c62=float(s_list[-5]) c63=float(s_list[-4]) c64=float(s_list[-3]) c65=float(s_list[-2]) c66=float(s_list[-1]) if 'Kirchhoff stress' in line: s=file_output.next() s_list=s.split() s11=float(s_list[-3])*160.22 s12=float(s_list[-2])*160.22 s13=float(s_list[-1])*160.22 s=file_output.next() s_list=s.split() s22=float(s_list[-2])*160.22 s23=float(s_list[-1])*160.22 s=file_output.next() s_list=s.split() s33=float(s_list[-3])*160.22 # change unit Gpa to kg/m s^2 source_stress='{}\n {}\n {}\n {}\n {}\n {}'.format(s11,s12,s13,s22,s23,s33) ibm = (c11+2*c12)/3*160.22 excess = math.sqrt((c12-c21)**2 + (c13-c31)**2 + (c14-c41)**2 + (c51-c15)**2 + (c61-c16)**2 + (c23-c32)**2 + (c24-c42)**2 + (c25-c52)**2 + (c26-c62)**2 + (c34-c43)**2 + (c35-c53)**2 + (c36-c63)**2 + (c45-c54)**2 + (c46-c64)**2 + (c56-c65)**2) # write output file for line in f: if '_latticetype_' in line: s=line f_result.write(s.replace('_latticetype_', x.latticetype)) elif '_source-ibm_' in line: s=line f_result.write(s.replace('_source-ibm_',str(ibm))) elif '_source-excess_' in line: s=line f_result.write(s.replace('_source-excess_',str(excess))) elif '_source-value_' in line: s=line f_result.write(s.replace('_source-value_',str(lat))) elif '_spacegroup_' in line: s=line f_result.write(s.replace('_spacegroup_',x.space_group)) elif '_basisatomcoor_' in line: s=line f_result.write(s.replace('_basisatomcoor_',x.basisatomcoords)) elif '_specieslist_' in line: s=line f_result.write(s.replace('_specieslist_',x.specieslist)) elif '_source-c44_' in line: s=line f_result.write(s.replace('_source-c44_',str(c44*160.22))) elif '_source-c11_' in line: s=line f_result.write(s.replace('_source-c11_',str(c11*160.22))) elif '_source-c12_' in line: s=line f_result.write(s.replace('_source-c12_',str(c12*160.22))) elif '_source-stress_' in line: s=line f_result.write(s.replace('_source-stress_',source_stress)) else: s=line f_result.write(s) #---------------------------------------------------------------------------------- # # Main function # #---------------------------------------------------------------------------------- class Input: pass x = Input() read_input(x) get_element_number(x) get_lattice_info(x) get_input(x) get_input_equi() # the path of cb-kim.x file might be different for different machine. It depends on where you install the cb-kim. cmd = 'cb-kim.x < output/equi.in > output/equi.out' os.system(cmd) check_equi(x) cmd = 'cb-kim.x < output/cbkim.in > output/cbkim.out' os.system(cmd) get_output(x)