#!/usr/bin/env python3 # """ # This Test Driver computes the dislocation core energy of a cubic crystal at zero temperature and a given stress state for a specified dislocation core cut-# off radius. # First, it generates dipole dislocation cores with distances of 10*a 20*a 30*a 40*a 50*a (where a is the lattice constant), by using MD++. # Second, by specifying a dislocation core radius, it will calculate the dislocation core energy. # Author: Yichen Qian, Villanova University 2020.Mar # """ import os import math import numpy as np import matplotlib.pyplot as plt import subprocess import gzip import shutil from kim_query import get_lattice_constant_cubic from kim_query import get_cohesive_energy_cubic from kim_query import get_elastic_constants_isothermal_cubic from kim_property import kim_property_create, kim_property_modify, kim_property_dump def stripquotes(matchobj): return matchobj.group(1) # Define the lammps running function def run_lammps(infile, outfile): """Run LAMMPS with given input file and write the output to outfile""" with open(outfile, "w") as outfl: try: subprocess.check_call(["lammps", "-in", infile], shell=False, stdout=outfl) except subprocess.CalledProcessError: extrainfo = "" try: with open("log.lammps") as f: extrainfo = f.read() except IOError: extrainfo = "no log file" raise Exception("LAMMPS did not exit properly:\n" + extrainfo) # Define the lattice type by full name(will use in generating dipole structure) def get_lattice_type_full_name(lattice_type): full_names = {"bcc": "body-centered-cubic", "fcc": "face-centered-cubic", "diamond": "diamond-cubic", "sc": "simple_cubic"} return full_names[lattice_type] def get_basis_atom_coordinates(lattice_type): basis_atom_coordinates = { "bcc": [[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]], "diamond": [ [0.0, 0.0, 0.0], [0.25, 0.25, 0.25], [0.5, 0.5, 0.0], [0.75, 0.75, 0.25], [0.5, 0.0, 0.5], [0.75, 0.25, 0.75], [0.0, 0.5, 0.5], [0.25, 0.75, 0.75], ], "fcc": [[0.0, 0.0, 0.0], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]], "sc": [[0.0, 0.0, 0.0]], } return basis_atom_coordinates[lattice_type] def format_basis_atom_coordinates(basis_atom_coordinates): """ Takes a nested list of basis atom coordinates and creates a list that can be splatted into a kim_property modify command """ output = [] for ind, atom in enumerate(basis_atom_coordinates): output += ["source-value", str(ind + 1), "1:3", *list(map(str, atom))] return output def format_species(element, lattice_type): """ Produces a list of appropriate length containing the (single) element provided as input that can be splatted into a kim_property modify command """ num_basis_atoms = {"bcc": 2, "diamond": 8, "fcc": 4, "sc": 1} Natoms = num_basis_atoms[lattice_type] return ["source-value", f"1:{Natoms}", *[element] * Natoms] def get_space_group(lattice_type): space_groups = {"fcc": "Fm-3m", "bcc": "Im-3m", "sc": "Pm-3m", "diamond": "Fd-3m"} return space_groups[lattice_type] # # File: madsum_iso.m # # # # Last Modified : Sun Aug 8 16:09:39 2004 # # Wei Cai, caiwei@stanford.edu # # # # References: # # # # Wei Cai, Vasily V. Bulatov, Jinpeng Chang, Ju Li and Sidney Yip, # # Periodic image effects in dislocation modelling, # # Philosophical Magazine A, 83, 539 (2003). # # ( for regularization of conditional convergence ) # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # This function calculates the elastic energy of a dislocation # dipole in a supercell under periodic boundary conditions # # Geometry: # # (x1,y1) A------------------------------------- # \ \ # \ P \ # \ (-) (x2r,y2r) \ # \ / \ # \ / \ # (+)-----------------------------------B (x0,0) # # Input: # mu: shear modulus (in eV/A^3, 1eV/A^3 = 160.22GPa ) # nu: Poisson's ratio # b = [ bx, by, bz ] Burgers vector # rc: cut-off radius # coords = [ x0, x1, y1, x2r, y2r ] # x0 x coord of B (0, infty) (in A) # x1 x coord of A (-infty,infty) (in A) # y1 y coord of A (0, infty) # x2r reduced x coord of P (-0.5, 0.5] # y2r reduced y coord of P (-0.5, 0.5] # trunc = [ xcut, ycut ] truncation paramters # -xcut:xcut times -ycut:ycut number of images will be added # # Output: # Eel: total elastic energy, Eel = Eprm + Eimg # Eprm: primary dipole energy # Eimg: image energy # Esum = [ lnR, xxR2, yyR2, xyR2 ] : intermediate results # lnR : sum of ln(R/rc) # xxR2: sum of x*x/R^2 # yyR2: sum of y*y/R^2 # xyR2: sum of x*y/R^2 # def madsum_iso_non_singular(mu, nu, b, rc, coords, trunc, theta): x0 = coords[0] x1 = coords[1] y1 = coords[2] x2r = coords[3] y2r = coords[4] xcut = trunc[0] ycut = trunc[1] c1 = np.array(np.transpose([x0, 0])) c2 = np.array(np.transpose([x1, y1])) ds = np.array(np.transpose([x2r, y2r])) dr = np.dot([c1, c2], ds) R = np.linalg.norm(dr) dt = dr / R # prm = np.array([np.log(R / rc), dt[0] * dt[0], dt[1] * dt[1], dt[0] * dt[1]]) Efac = -(mu / 4 / pi) * np.array( [ (b[0] * b[0] + b[1] * b[1]) / (1 - nu) + b[2] * b[2], b[1] * b[1] / (1 - nu), b[0] * b[0] / (1 - nu), -2 * b[0] * b[1] / (1 - nu), ] ) Eprm = ( (mu / 2 / pi) * (np.linalg.norm(b) ** 2) * ( ((math.cos(theta) ** 2) + (math.sin(theta) ** 2) / (1 - nu)) * np.log(R / rc) + ((math.cos(theta) ** 2) * 1 / 2) ) ) image = np.array([[0, 0, 2], [dr[0], dr[1], -1], [-dr[0], -dr[1], -1]]) Esum = np.array([0, 0, 0, 0]) for i in range(-xcut, xcut + 1): for j in range(-ycut, ycut + 1): if i != 0 or j != 0: offset = i * c1 + j * c2 for k in range(0, len(image[:, 0])): dR = offset + np.transpose(image[k, 0:2]) R = np.linalg.norm(dR) dt = dR / R Esum = Esum + image[k, 2] * np.array( [np.log(R / rc), dt[0] * dt[0], dt[1] * dt[1], dt[0] * dt[1]] ) ghost = np.array( [ [c1[0] / 2 - dr[0] / 2, c1[1] / 2 - dr[1] / 2, x2r], [-c1[0] / 2 - dr[0] / 2, -c1[1] / 2 - dr[1] / 2, -x2r], [c1[0] / 2 + dr[0] / 2, c1[1] / 2 + dr[1] / 2, -x2r], [-c1[0] / 2 + dr[0] / 2, -c1[1] / 2 + dr[1] / 2, x2r], [c2[0] / 2 - dr[0] / 2, c2[1] / 2 - dr[1] / 2, y2r], [-c2[0] / 2 - dr[0] / 2, -c2[1] / 2 - dr[1] / 2, -y2r], [c2[0] / 2 + dr[0] / 2, c2[1] / 2 + dr[1] / 2, -y2r], [-c2[0] / 2 + dr[0] / 2, -c2[1] / 2 + dr[1] / 2, y2r], ] ) Egst = np.array([0, 0, 0, 0]) for i in range(-xcut, xcut + 1): for j in range(-ycut, ycut + 1): offset = i * c1 + j * c2 for k in range(0, len(ghost[:, 0])): dR = offset + np.transpose(ghost[k, 0:2]) R = np.linalg.norm(dR) dt = dR / R Egst = Egst - ghost[k, 2] * np.array( [np.log(R / rc), dt[0] * dt[0], dt[1] * dt[1], dt[0] * dt[1]] ) Esum = Esum + Egst Eimg = sum(np.multiply(Efac, Esum)) Eel = Eprm + Eimg return [Eel, Eprm, Eimg, Esum] def get_miller_index(inder_string): newstr = "".join((ch if ch in "0123456789.-e" else " ") for ch in inder_string) miller_index = [float(i) for i in newstr.split()] index = ( str(miller_index[0]) + " " + str(miller_index[1]) + " " + str(miller_index[2]) ) return index def get_miller_index_vector(inder_string): newstr = "".join((ch if ch in "0123456789.-e" else " ") for ch in inder_string) miller_index = [float(i) for i in newstr.split()] return miller_index def get_the_dislocation_core_radius(Core_Radius_string): newstr = "".join( (ch if ch in "0123456789.-e" else " ") for ch in Core_Radius_string ) Core_Radius_index = [float(i) for i in newstr.split()] return Core_Radius_index # Read elements, masses, and model from stdin element = input("Elements: ") print(element) mass = input("Masses in g/mol (must match order of elements above): ") print(mass) model = input("Model: ") print(model) lattice_type = input("Lattice type: ") print(lattice_type) dislocation_core_radius_seq = input("dislocation core radius: ") print(dislocation_core_radius_seq) burgers_vector_direction = input("burgers vector direction: ") print(burgers_vector_direction) slip_plane_miller_indices = input("slip plane miller indices: ") print(slip_plane_miller_indices) dislocation_line_direction = input("dislocation line direction: ") print(dislocation_line_direction) print("") # Some directories we need THIS_DIR = os.path.dirname(__file__) MAIN_LAMMPS_INPUT_TEMPLATE = open( os.path.join(THIS_DIR, "min_cg_npt_w.in.template") ).read() DISLOCATION_DIPOLE_TEMPLATE = open( os.path.join(THIS_DIR, "DISLOCATION_DIPOLE.TEMPLATE") ).read() INDIR = os.path.join("output", "lammps_inputs") OUTDIR = os.path.join("output", "lammps_output_log") DUMPDIR = os.path.join("output", "dumpfile") # Ensure the directories we need are created try: os.makedirs(INDIR) except OSError: pass try: os.makedirs(OUTDIR) except OSError: pass try: os.makedirs(DUMPDIR) except OSError: pass # Remove the spece from input. species = " ".join(element.split()) # Calculate the equilibrium lattice constant equilibrium_lattice_constant = get_lattice_constant_cubic( [model], [lattice_type], [species], ["angstrom"] )[0] print("equilibrium lattice constant is", equilibrium_lattice_constant) # Calculate the equilibrium lattice constant cohesive_energy_final = get_cohesive_energy_cubic( [model], [lattice_type], [species], ["eV"] )[0] print("cohesive energy is", cohesive_energy_final) # Calculate the elastic constants elastic_constants = get_elastic_constants_isothermal_cubic( [model], [lattice_type], [species], ["GPa"] ) # The elastic constant transfer. C11 = elastic_constants[0] C12 = elastic_constants[1] C44 = elastic_constants[2] bulk_modulus = 1 / 3 * (C11 + 2 * C12) G_v = (3 * C44 + C11 - C12) / 5 G_r = (C44 * 5 * (C11 - C12)) / (4 * C44 + 3 * (C11 - C12)) shear_modulus = (G_v + G_r) / 2 poission_ratio = (3 * bulk_modulus - 2 * shear_modulus) / ( 2 * (3 * bulk_modulus + shear_modulus) ) # set necessary input parameters in input files latticestructure = get_lattice_type_full_name(lattice_type) length_of_b01 = [10, 20, 30, 40, 50] final_value = [] Main_result = [] burgers_vector = get_miller_index(burgers_vector_direction) burgers_vector_vector = get_miller_index_vector(burgers_vector_direction) print(burgers_vector) slip_plane_direction = get_miller_index(slip_plane_miller_indices) slip_plane_direction_vector = get_miller_index_vector(slip_plane_miller_indices) print(slip_plane_direction) dislocation_line_direction = get_miller_index(dislocation_line_direction) dislocation_line_direction_vector = get_miller_index_vector(dislocation_line_direction) print(dislocation_line_direction) dislocation_core_radius = get_the_dislocation_core_radius(dislocation_core_radius_seq) print(dislocation_core_radius) E_core_final = [] # The loop for generate total energy of supercell containing dipole structure for length_b01 in length_of_b01: Dipole_creation_file = os.path.join( INDIR, "DISLOCATION_DIPOLE.TEMPLATE_" + str(length_b01) + ".scripts" ) with open(Dipole_creation_file, "w") as creation_file: creation_file.write( DISLOCATION_DIPOLE_TEMPLATE.format( DIRNAME=INDIR, length_b01=length_b01, latticestructure=latticestructure, lattice_constant=equilibrium_lattice_constant, poission_ratio=poission_ratio, burgers_vector=burgers_vector, slip_plane_direction=slip_plane_direction, dislocation_line_direction=dislocation_line_direction, ) ) subprocess.call(["md_gpp", Dipole_creation_file]) with gzip.open( os.path.join( INDIR + "/dipole_ref/", "makedp_b01_" + str(length_b01) + ".lammps.gz" ), "rb", ) as f_in: with open( os.path.join( INDIR + "/dipole_ref/", "makedp_b01_" + str(length_b01) + ".lammps" ), "wb", ) as f_out: shutil.copyfileobj(f_in, f_out) makedpfile = open( os.path.join( INDIR + "/dipole_ref/", "makedp_b01_" + str(length_b01) + ".lammps" ) ).read() infile = os.path.join(INDIR, "min_cg_npt_w.in.template") ## xyzfile = os.path.join(INDIR, "makedp_b01_" + str(length_b01) + ".lammps") logfile = os.path.join( OUTDIR, "equilibrate_npt_disloc_dislocation_b01_" + str(length_b01) + ".log" ) dumpfile = os.path.join( DUMPDIR, "makedp_npt_dislocation_b01_" + str(length_b01) + ".dump.gz" ) ### REGULAR EXPRESSIONS FOR MATCHING LAMMPS OUTPUT # Read xyz file, convert the cell to LAMMPS' convention and create a # dump file to read in # NOTE: using a dump file rather than a data file means that it's not # sensitive to the specific atom_style in the event that we're # running with an SM, e.g. if the atom_style of the SM is # 'charge', charges of zero will simply be assigned to each atom # since they're not specified in the dumpfile with open(xyzfile, "w") as xyz_file: xyz_file.write(makedpfile.format()) # Create the LAMMPS input file with open(infile, "w") as in_file: in_file.write( MAIN_LAMMPS_INPUT_TEMPLATE.format( modelname=model, symbol=element, xyzfile=xyzfile, mass=mass, dumpfile=dumpfile, lattice_type=lattice_type, ) ) run_lammps(infile, logfile) ### Now process the output and dumpfile for relevant information # Get Dislocation energy with open(logfile) as outfl: lines = outfl.readlines() try: temp_value = [] for line in lines: line_split = line.split() if len(line_split) <= 10: line_split = [] zzz = [] if line_split: try: for item in line_split: zzz.append(float(item)) temp_value.append(zzz) except: continue zz = 0 for index in range(0, len(temp_value)): if temp_value[index][0] >= temp_value[zz][0]: zz = index Main_result = temp_value[zz] except AttributeError: raise Exception( "Error: Failed to find the dislocation energy in the LAMMPS output log" ) print(Main_result) useful_value = [ length_b01, Main_result[1], -cohesive_energy_final, 0.0000000e00, Main_result[4], Main_result[8], Main_result[9], Main_result[10], Main_result[11], Main_result[12], Main_result[13], ] final_value.append(useful_value) # Calculate the dislocation core energy at different cut-off core radius for core_radius in dislocation_core_radius: rlatt = equilibrium_lattice_constant mu = C44 # mu in units to compute "nu" input vailable nu = poission_ratio # mu and nu shuold be input for the test mu = C44 / 160.22 # mu in units required to MadSum pi = math.pi mod_b = np.linalg.norm([i / 2 * rlatt for i in [1, 1, 1]]) # ## geometry of the simulation box x = [-1, 1, 0] # x: dipole direction y = [-1, -1, 2] # y: z = [1 / 2, 1 / 2, 1 / 2] # z: line direction b = [1 / 2, 1 / 2, 1 / 2] # b: burgers vector # b = z # dislocation dipole # # b01 = 10; # [-1 1 0] direction b02 = 50 # [ -1 -1 2] b03 = 5 # [ 1 1 1] direction # ## Obtaining angle between axes and b = [ 1 1 1]/2 delta = math.acos(np.dot(b, x) / (np.linalg.norm(b) * np.linalg.norm(x))) print("delta = angle(b,x) [degrees]= ", delta * 180 / pi) # ; display(str) # phi = math.acos(np.dot(b, y) / (np.linalg.norm(b) * np.linalg.norm(y))) print("phi = angle(b,y) [degrees]= ", phi * 180 / pi) # ; display(str) theta = math.acos(round(np.dot(b, z) / (np.linalg.norm(b) * np.linalg.norm(z)), 8)) print("theta = angle(b,z) [degrees]= ", theta * 180 / pi) # ; display(str) # ## burgers vector in MADSUM format # b_madsum = np.array( [ i * np.linalg.norm(b) * rlatt for i in [math.cos(delta), math.cos(phi), math.cos(theta)] ] ) alatt = rlatt * math.sqrt(2.0) # <110> direction is the dipole direction -> b01 blatt = rlatt * math.sqrt(6.0) # <112> direction -> b02 clatt = rlatt * math.sqrt(3.0) / 2.0 # 1/2[111] direction is the burgers vector -> b03 N = core_radius rc = N * mod_b cut = np.array([10, 10]) # rc =1; b = b_madsum r_b01 = np.array([]) Eel = np.array([]) Eprm = np.array([]) Eimg = np.array([]) for b01 in [10, 20, 30, 40, 50]: # 10, number of components of b03 coord = np.array([b01 * alatt, 0, b02 * blatt, 0.5, 0]) print(mu, nu, b, rc, coord, cut, theta) ### singular solution # [p1, p2, p3, Esum] = madsum_iso_dislocation(mu, nu, b, rc, coord, cut) # ### non singular solution [p1, p2, p3, Esum] = madsum_iso_non_singular(mu, nu, b, rc, coord, cut, theta) print(p1, p2, p3) r_b01 = np.append(r_b01, b01) Eel = np.append(Eel, p1) Eprm = np.append(Eprm, p2) Eimg = np.append(Eimg, p3) print(r_b01, Eel, Eprm, Eimg, Esum) # fig1 = plt.figure() # plt.plot(r_b01,Eel,'r--',r_b01,Eprm,'g--',r_b01,Eimg,'b--') # plt.legend(['E.elast','Eprm','Eimg'],loc=0) # plt.ylabel('E (eV/A)',fontsize=14) # plt.xlabel('b01 ( units of a_{110})',fontsize=14) # plt.title('W CEA vs4 ; DISLOCATION; rc='+str(N),fontsize=14) # plt.grid(True) # plt.savefig('output/figure1_'+str(core_radius)+'.png') # figure printing scripts------------------------ # --------------------------------------------------- mu = C44 # mu in units to compute "nu" B = bulk_modulus # input vailable ## Residual stress # E_residual=residual_stress_cea(rlatt,mu,nu,b01,b02,b03,alatt,blatt,clatt); # mu = 161; nu = 0.284; mu = mu * 1000 # in MPa...161 GPa lamb = (2 * nu * mu) / (1 - 2 * nu) c_44 = mu c_12 = lamb c_11 = 2 * c_44 + c_12 C = np.zeros([6, 6]) C[0, 0] = c_11 C[1, 1] = c_11 C[2, 2] = c_11 C[0, 1] = c_12 C[0, 2] = c_12 C[1, 0] = c_12 C[1, 2] = c_12 C[2, 0] = c_12 C[2, 1] = c_12 C[3, 3] = c_44 C[4, 4] = c_44 C[5, 5] = c_44 C # GPa S = np.linalg.inv(C) ## E_min_dislocation_npt_cea = np.loadtxt(os.path.join(THIS_DIR , "E_min_dislocation_npt_cea.dat")) E_min_dislocation_npt_cea = np.array(final_value) print(E_min_dislocation_npt_cea) # load E_min_dislocation_npt_cea # Sigma = E_min_dislocation_npt_cea (:,6:11)/10 # for j=1:5 # #Sigma = E_min_dislocation_npt_cea (:,5:11)/10 ; #Pij in LAMMPS (metal units) are in bar=0.1 MPa # E_residual(j)= Sigma(j,:)*S*Sigma(j,:)'; # E_residual in MPa = MJ/m^3 # #Sigma # end E_residual = [] for j in range(0, 5): Sigma = [] for i in range(0, 6): Sigma.append( E_min_dislocation_npt_cea[j, 5 + i] / 10 ) # Pij in LAMMPS (metal units) are in bar=0.1 MPa E_residual.append( np.dot(np.dot(Sigma, S), np.transpose(Sigma)) ) # Sigma*S*Sigma'(original); # E_residual in MPa = MJ/m^3 b01 = [10, 20, 30, 40, 50] Vol = [] for m in [0, 1, 2, 3, 4]: Vol.append( b01[m] * alatt * (10 ** -10) * b02 * blatt * (10 ** -10) * b03 * clatt * (10 ** -10) ) # Vol in m^3 E_residual[m] = E_residual[m] * (10 ** 6) * Vol[m] # E_residual in J E_residual[m] = E_residual[m] / (1.602e-19) # E_residual in eV # fig2 = plt.figure() # plt.plot(b01,E_residual,'r--') # plt.legend(['E.residual'],loc=0) # plt.ylabel('E.residual (eV)',fontsize=14) # plt.xlabel('b01 ( units of a_{110})',fontsize=14) # plt.grid(True) # plt.savefig('output/figure2_'+str(core_radius)+'.png') # -------figure part ------------ # ventana1=figure #open window before creating the plot # set(ventana1,'visible','on') # # ejes1=axes # we create the object axes # p1 = plot(b01,E_residual) # legend('E.residual') # axis([-58.7 58.7 -30 30]) # axis([-40 -20 58.7 -30 30]) # grid on # xlabel('b01 (units of a_{110})','fontsize',14) # ylabel('E.residual (eV)','fontsize',14) # title(['rc=',str(N),'b'],'fontsize',14) E_atm = ( np.array(E_min_dislocation_npt_cea[:, 4]) - np.array( np.multiply( E_min_dislocation_npt_cea[:, 1], E_min_dislocation_npt_cea[:, 2] ) ) - np.array(np.transpose(E_residual)) / 2 ) / b03 # eV/b exp. 5.16 book Wei Cai and Vasily Bulatov E_elastic = np.array(np.transpose(Eel)) # (eV/A) E_core = (E_atm - E_elastic * clatt) / 2 # (eV/b) #### units of (eV/A) E_atm_A = ( np.array(E_min_dislocation_npt_cea[:, 4]) - np.array( np.multiply( E_min_dislocation_npt_cea[:, 1], E_min_dislocation_npt_cea[:, 2] ) ) - np.array(np.transpose(E_residual) / 2) ) / (b03 * clatt) # eV/A exp. 5.16 book Wei Cai and Vasily Bulatov E_core_A = (E_atm_A - E_elastic) / 2 #### # # fig3 = plt.figure() plt.plot( E_min_dislocation_npt_cea[:, 0], E_atm_A, "g--", E_min_dislocation_npt_cea[:, 0], E_elastic, "r--", E_min_dislocation_npt_cea[:, 0], E_core_A, "b--", ) plt.legend(["E.atomistic_A", "Eelastic_A", "Ecore_A"], loc=0) plt.ylabel("E (eV/A)", fontsize=14) plt.xlabel("b01 ([-110])", fontsize=14) plt.title("DISLOCATION; rc= " + str(N) + "b", fontsize=14) plt.grid(True) plt.savefig("output/figure3_" + str(N) + ".png") # Find Dislocation core energy # set a convergence error # print('the dislocation core energy is', E_core_A) # E_diff_core = abs(np.diff(E_core_A)) # E_core_index = np.nonzero(E_diff_core < 0.01) # set a convergence percentage # E_core_final = E_core_A[E_core_index[0]] # print(E_core_final) E_core_final.append(np.polyfit(E_min_dislocation_npt_cea[:, 0], E_core_A, 0)) print(E_core_final) log_core_radius = np.log(dislocation_core_radius) [a_index, b_index] = np.polyfit(log_core_radius, E_core_final, 1) xnew = np.linspace(min(dislocation_core_radius), max(dislocation_core_radius), 300) y = a_index * np.log(xnew) + b_index fig4 = plt.figure() plt.plot(dislocation_core_radius, E_core_final, "*", xnew, y, "g-") plt.legend(["Core energy fitting function"], loc=0) plt.ylabel("E_core[\u03B8=0,a](eV/A)", fontsize=14) plt.xlabel("a (units of b)", fontsize=14) # plt.title('DISLOCATION; rc= '+str(N)+'b',fontsize=14) plt.grid(True) plt.savefig("output/figure4.png") # Create a property instance prop_inst = kim_property_create(1, "dislocation-core-energy-cubic-crystal-npt") # Set all the key-value pairs for this property instance prop_inst = kim_property_modify( prop_inst, 1, "key", "short-name", "source-value", "1", lattice_type, "key", "species", *format_species(element, lattice_type), "key", "a", "source-value", equilibrium_lattice_constant, "source-unit", "angstrom", "key", "basis-atom-coordinates", *format_basis_atom_coordinates(get_basis_atom_coordinates(lattice_type)), "key", "space-group", "source-value", get_space_group(lattice_type), "key", "cauchy-stress", "source-value", "1:6", *[0.0] * 6, "source-unit", "GPa", "key", "slip-plane-miller-indices", "source-value", "1:3", *slip_plane_direction_vector, "key", "dislocation-line-direction", "source-value", "1:3", *dislocation_line_direction_vector, "key", "burgers-vector-direction", "source-value", "1:3", *burgers_vector_vector, "key", "dislocation-core-radius", "source-value", "1:9", *dislocation_core_radius, "key", "core-energy", "source-value", "1:9", *E_core_final, "source-unit", "eV/A", "digits", "5" ) # Dump the results in a file with open("output/results.edn", "w") as fp: kim_property_dump(prop_inst, fp)