#!/usr/bin/env python3 """ Author: Brandon Runnels Date: June 2019 Institution: Department of Mechanical and Aerospace Engineering University of Colorado Colorado Springs Email: brunnels@uccs.edu Description: This file contains a single function that computes the grain boundary (GB) energy for a range of tilt angles and an arbitrarily selected tilt axis. The user specifies the initial orientations and the range of tilt angles, as well as a maximum repeat- ing unit cell size. The script then creates a LAMMPS input file for all possible tilt boundaries within that set of parameters, runs LAMMPS, collects output, and returns the theta ranges, en- ergies, sigma values, and minimum seperation dist- ances to the calling function. For details on the input parameters, please see README.txt Copyright: 2016,2017,2019 University of Colorado Colorado Springs """ import os from math import ceil from fractions import Fraction import pickle import textwrap import numpy from scipy import linalg def compute_gb_tilt_energy( output_dir, model_name, species, a1, a2, a3, theta_min, theta_max, max_denominator, lattice_style, lattice_constant, a_curve, cohesive_energy_curve, offset_grid_frac, x_repeat=1, z_repeat=1, min_cell_height=10, offsets_top=-1, offsets_bottom=-1, num_mpi_processors=1, cohesive_energy=0, mass=0, isolated_atom_energy=0, gzip=False, verbose=False, ): # Check if the user has specified a custom grid of offset translation vectors to # loop over for each tilt angle grid_search = True if offsets_top != -1 and offsets_bottom != -1: grid_search = False elif (offsets_top == -1 and offsets_bottom != -1) or ( offsets_top != -1 and offsets_bottom == -1 ): raise ValueError( "Error: you must specify BOTH offsets_top and offsets_bottom or NEITHER offsets_top or offsets_bottom" ) theta_res = 1000 # use cohesive energy curve to estimate cutoff distance if cohesive_energy_curve[-1] > 0 : # did not compress far enough to reach zero energy point # do linear approximation of last 2 points instead aneg=a_curve[-1] * 1e10 # convert from m to angstroms eneg=cohesive_energy_curve[-1]/1.602177e-19 #convert from J to eV apos=a_curve[-2] * 1e10 # convert from m to angstroms epos=cohesive_energy_curve[-2]/1.602177e-19 #convert from J to eV # solve for the x interecept of the linear approximation of the cohesive energy curve cutoff_distance_si = -epos * ((aneg - apos)/(eneg - epos)) + apos cutoff_distance = cutoff_distance_si / lattice_constant # convert to units of lattice constant else : # Find the separation where cohesive energy changes sign from compressive stress for count,energy in enumerate(cohesive_energy_curve) : energy=energy / 1.602177e-19 # convert from J to eV dist = a_curve[count]/lattice_constant if energy < 0 and dist < 1: # make sure we're at the compressive zero energy point at separations < a aneg=a_curve[count] * 1e10 # convert from m to angstroms eneg=energy apos=a_curve[count-1] * 1e10 # convert from m to angstroms epos=cohesive_energy_curve[count-1]/ 1.602177e-19 # convert from J to eV # solve for the x interecept of the linear approximation of the cohesive energy curve cutoff_distance_si = -epos * ((aneg - apos)/(eneg - epos)) + apos cutoff_distance = cutoff_distance_si / lattice_constant # convert to units of lattice constant break # unit vectors corresponding to the un-tilted configuration a1 = numpy.asarray(a1) a2 = numpy.asarray(a2) a3 = numpy.asarray(a3) e1 = a1 / linalg.norm(a1) e2 = a2 / linalg.norm(a2) e3 = a3 / linalg.norm(a3) # check for orthogonality if abs(a1.dot(a2)) > 1e-8 or abs(a2.dot(a3)) > 1e-8 or abs(a3.dot(a1)) > 1e-8: raise ValueError("Vectors (a1, a2, a3) do not form an orthogonal set") # check for right hand rule if ( numpy.cross(e1, e2).dot(e3) < 0 or numpy.cross(e2, e3).dot(e1) < 0 or numpy.cross(e3, e1).dot(e2) < 0 ): raise ValueError("Vectors (a1, a2, a3) do not obey right hand rule") # initialize output arrays thetas = [] gb_energies = [] minimum_distances = [] sigmas = [] # Integers to remember the last CSL configuration to avoid repeating calculations mold = 0 nold = 0 pold = 0 qold = 0 aux_results_info = {} # Iterate thorough all thetas in the range for th in numpy.linspace(theta_min, theta_max, theta_res): # Compute the repeating distance in the x direction and snap to the nearest finite integer rationalization if ( -1 < (numpy.tan(numpy.radians(th / 2)) * linalg.norm(a1) / linalg.norm(a2)) < 1 ): frac = Fraction( numpy.tan(numpy.radians(th / 2)) * linalg.norm(a1) / linalg.norm(a2) ).limit_denominator(max_denominator) m = int(frac.denominator) n = int(frac.numerator) else: frac = Fraction( linalg.norm(a2) / numpy.tan(numpy.radians(th / 2)) / linalg.norm(a1) ).limit_denominator(max_denominator) m = int(frac.numerator) n = int(frac.denominator) # Compute the repeating distance in the y direction and snap to the nearest finite integer rationalization if ( -1 < (numpy.tan(numpy.radians(th / 2)) * linalg.norm(a2) / linalg.norm(a1)) < 1 ): frac = Fraction( numpy.tan(numpy.radians(th / 2)) * linalg.norm(a2) / linalg.norm(a1) ).limit_denominator(max_denominator) p = int(frac.numerator) q = int(frac.denominator) else: frac = Fraction( linalg.norm(a1) / numpy.tan(numpy.radians(th / 2)) / linalg.norm(a2) ).limit_denominator(max_denominator) p = int(frac.denominator) q = int(frac.numerator) # Avoid repeating previously found rationalizations if [m, n, p, q] == [mold, nold, pold, qold]: continue mold = m nold = n pold = p qold = q # Find the actual theta for the rationalized angle thetaA = numpy.degrees( 2 * numpy.arctan2(n * linalg.norm(a2), m * linalg.norm(a1)) ) thetaB = numpy.degrees( 2 * numpy.arctan2(p * linalg.norm(a1), q * linalg.norm(a2)) ) # Sometimes the top and bottom crystals return different angles; in this case, skip if abs(thetaA - thetaB) > 1e-8: continue theta = thetaA if verbose: print("\n========================================") print("theta = " + str(theta)) print("========================================\n") # Compute integral dimensions of the unit cell lambda1 = m * linalg.norm(a1) * numpy.cos( numpy.radians(theta / 2) ) + n * linalg.norm(a2) * numpy.sin(numpy.radians(theta / 2)) lambda2 = p * linalg.norm(a1) * numpy.sin( numpy.radians(theta / 2) ) + q * linalg.norm(a2) * numpy.cos(numpy.radians(theta / 2)) lambda3 = linalg.norm(a3) # Convert integral dimensions to actual lengths x_cell = lambda1 * lattice_constant y_cell = lambda2 * lattice_constant z_cell = lambda3 * lattice_constant # Compute the minimum number of repeats in the y direction to satisfy the minimum height requirement y_repeat = int(ceil(min_cell_height * lattice_constant / y_cell)) # Compute coordinates of the box xmin = 0.0 xmax = x_cell * x_repeat ymin = -y_cell * y_repeat - 0.000000001 * lattice_constant ymax = y_cell * y_repeat + 0.000000001 * lattice_constant zmin = 0.0 zmax = z_cell * z_repeat # Compute rotated direction vectors X1 = m * a1 - n * a2 X2 = m * a1 + n * a2 Y1 = p * a1 + q * a2 Y2 = -p * a1 + q * a2 Z1 = a3 Z2 = a3 if grid_search: # Form grid over x_cell and z_cell (recall that offset_grid_frac is # expressed as a fraction of the lattice constant). We subtract # 0.05*grid_mesh from the upper limit of the grid extent so as to # ensure that the shifts corresponding to lambda1 and lambda3 # themselves are not included since they're equivalent to zero # shift by symmetry grid_mesh = offset_grid_frac * lattice_constant x_shifts = numpy.arange( 0, x_cell - 0.05 * grid_mesh, grid_mesh, dtype=numpy.float64 ) z_shifts = numpy.arange( 0, z_cell - 0.05 * grid_mesh, grid_mesh, dtype=numpy.float64 ) offsets_top = [[0.5 * x, 0, 0.5 * z] for x in x_shifts for z in z_shifts] offsets_bottom = [ [-0.5 * x, 0, -0.5 * z] for x in x_shifts for z in z_shifts ] minimum_energy = float("inf") # Create the directory structure for dumping the output directory = os.path.join(output_dir, "dump_{:0>8.4f}".format(theta)) if not os.path.exists(directory): os.makedirs(directory) lammps_input_file = os.path.join(directory, "generated_input.gb") lammps_output_file = os.path.join(directory, "run_output.gb") # Files written by LAMMPS to communicate relevant results energy_file = os.path.join(directory, "energy.out") numatoms_file = os.path.join(directory, "numatoms.out") mindistance_file = os.path.join(directory, "mindistance.out") dump_file = os.path.join(directory, "dump_post.cfg" + (".gz " if gzip else "")) success_file = os.path.join(directory, "success.out") # Container for storing extra info about the results. This isn't used in generating # the property instances that will go in results.edn, but is rather used to # create a pickle file in each dump directory for bookkeeping aux_results_info = [] for offset_top, offset_bottom in zip(offsets_top, offsets_bottom): aux_results_info_entry = {} LAMMPS_INPUT_TEMPLATE = textwrap.dedent( r""" clear kim init {model_name} metal unit_conversion_mode dimension 3 boundary p p p # Adjust input variables for possible units changes with # Simulator Models variable lattice_constant_converted equal {lattice_constant}*${{_u_distance}} variable xmin_converted equal {xmin}*${{_u_distance}} variable xmax_converted equal {xmax}*${{_u_distance}} variable ymin_converted equal {ymin}*${{_u_distance}} variable ymax_converted equal {ymax}*${{_u_distance}} variable zmin_converted equal {zmin}*${{_u_distance}} variable zmax_converted equal {zmax}*${{_u_distance}} # Simulation domain lattice {lattice_style} ${{lattice_constant_converted}} region whole block ${{xmin_converted}} ${{xmax_converted}} & ${{ymin_converted}} ${{ymax_converted}} & ${{zmin_converted}} ${{zmax_converted}} & units box create_box 2 whole # Upper region region upper block INF INF 0.0 ${{ymax_converted}} INF INF units box lattice {lattice_style} ${{lattice_constant_converted}} & orient x {lattice_orient_top_x} & orient y {lattice_orient_top_y} & orient z {lattice_orient_top_z} create_atoms 1 region upper group upper type 1 variable mass1_converted equal {mass}*${{_u_mass}} mass 1 ${{mass1_converted}} # Move atoms in upper region displace_atoms upper move {displace_upper} units box # Lower region region lower block INF INF ${{ymin_converted}} 0.0 INF INF units box lattice {lattice_style} ${{lattice_constant_converted}} & orient x {lattice_orient_bottom_x} & orient y {lattice_orient_bottom_y} & orient z {lattice_orient_bottom_z} create_atoms 2 region lower group lower type 2 # Move atoms in lower region displace_atoms lower move {displace_lower} units box variable mass2_converted equal {mass}*${{_u_mass}} mass 2 ${{mass2_converted}} # Interatomic potential and neighbor settings kim interactions {species} {species} # Set up neighbor list method variable neigh_skin_converted equal 2.0*${{_u_distance}} neighbor ${{neigh_skin_converted}} bin neigh_modify delay 10 check yes # Delete overlapping atoms variable delete_distance equal {offset_grid_frac}*${{lattice_constant_converted}} delete_atoms overlap ${{delete_distance}} all all # Solver min_style cg minimize 1e-6 1e-6 5000 10000 fix 1 all box/relax x 0.0 y 0.0 z 0.0 couple none vmax 0.001 minimize 1e-6 1e-6 1000 10000 unfix 1 # Variables used to rescale the positions and energies so that # the quantities in the dumpfile are in the original metal units # (angstrom and eV) even if we're running with a Simulator Model # that uses different units variable pe_metal equal "c_thermo_pe/v__u_energy" variable lx_metal equal lx/${{_u_distance}} variable ly_metal equal ly/${{_u_distance}} variable lz_metal equal lz/${{_u_distance}} variable press_metal equal "c_thermo_press/v__u_pressure" variable pxx_metal equal pxx/${{_u_pressure}} variable pyy_metal equal pyy/${{_u_pressure}} variable pzz_metal equal pzz/${{_u_pressure}} variable mass_metal atom mass/${{_u_mass}} # Compute centrosymmetry, particle energy, and min dist variables # Note that if the Model does not define particle energy, this quantity will # be zero. compute csym all centro/atom {lattice_style} compute particle_eng all pe/atom compute particle_engsum all reduce sum c_particle_eng compute csymsum all reduce sum c_csym compute distance all pair/local dist compute mindist all reduce min c_distance # Isolated atom energy in eV (no unit conversion necessary) variable isolated_atom_energy equal {isolated_atom_energy} variable particle_eng_metal atom "c_particle_eng/v__u_energy - v_isolated_atom_energy" variable csym_metal atom c_csym/(${{_u_distance}}*${{_u_distance}}) variable csymsum_metal equal c_csymsum/(${{_u_distance}}*${{_u_distance}}) variable mindist_metal equal c_mindist/${{_u_distance}} reset_timestep 0 run 0 thermo 0 thermo_style custom step v_pe_metal c_particle_engsum & v_lx_metal v_ly_metal v_lz_metal & press v_press_metal v_pxx_metal v_pyy_metal v_pzz_metal & v_mindist_metal v_csymsum_metal run 0 # Spit out number of atoms variable num_atoms equal "count(all)" print "${{num_atoms}}" file {numatoms_file} # Spit out total energy relative to sum of isolated energies, converted back to our # original metal units variable adjusted_pe_metal equal ${{pe_metal}}-${{num_atoms}}*${{isolated_atom_energy}} print "${{adjusted_pe_metal}} eV" file {energy_file} # Spit out distances, converted back to our original metal units print "${{mindist_metal}} Angstroms" file {mindistance_file} # Write out dump file write_dump all cfg{gzip} {dump_file} & mass type xs ys zs id v_mass_metal v_csym_metal v_particle_eng_metal & modify element {species} {species} print "This indicates that LAMMPS ran successfully" file {success_file} """ ) lammps_input = LAMMPS_INPUT_TEMPLATE.format( model_name=model_name, lattice_constant=lattice_constant, lattice_style=lattice_style, species=species, mass=mass, isolated_atom_energy=isolated_atom_energy, displace_upper=" ".join(map(str, offset_top)), displace_lower=" ".join(map(str, offset_bottom)), xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, offset_grid_frac=offset_grid_frac, lattice_orient_top_x=" ".join(map(str, X1)), lattice_orient_top_y=" ".join(map(str, Y1)), lattice_orient_top_z=" ".join(map(str, Z1)), lattice_orient_bottom_x=" ".join(map(str, X2)), lattice_orient_bottom_y=" ".join(map(str, Y2)), lattice_orient_bottom_z=" ".join(map(str, Z2)), gzip="/gz" if gzip else "", dump_file=dump_file, numatoms_file=numatoms_file, energy_file=energy_file, success_file=success_file, mindistance_file=mindistance_file, ) with open(lammps_input_file, "w") as lammps_input_fl: lammps_input_fl.write(textwrap.dedent(lammps_input)) # Run LAMMPS if num_mpi_processors > 1: os.system( "mpirun -np {} lammps -i {} > {}".format( num_mpi_processors, lammps_input_file, lammps_output_file ) ) elif num_mpi_processors == 1: os.system( "lammps -i {} > {}".format(lammps_input_file, lammps_output_file) ) else: raise ValueError("Error: invalid number of mpi processors") if not os.path.isfile(dump_file): raise RuntimeError( "Error: LAMMPS did not generate a dump file -- something probably " "went wrong." ) if not os.path.isfile(success_file): raise RuntimeError( "Error: there was an error in the LAMMPS run -- check the LAMMPS " "log file for more detail" ) # Collect output from LAMMPS that was recorded in files and delete them with open(energy_file) as f: energy = float(f.read().split()[0]) with open(numatoms_file) as f: numatoms = int(f.read()) with open(mindistance_file) as f: mindistance = float(f.read().split()[0]) # Compute grain boundary energy if theta == 0 and abs(cohesive_energy) < 1e-16: cohesive_energy = energy / numatoms excess_energy = energy - numatoms * (-abs(cohesive_energy)) gb_energy_ev_angstrom2 = ( excess_energy / (x_repeat * (xmax - xmin)) / (z_repeat * (zmax - zmin)) / 2 ) gb_energy_J_m2 = gb_energy_ev_angstrom2 * 16.0218 print("Energy: {}".format(gb_energy_J_m2)) # Record auxiliary information for posterity aux_results_info_entry["offset"] = (offset_top, offset_bottom) aux_results_info_entry["relaxed_energy_J_m2"] = gb_energy_J_m2 aux_results_info_entry["min_distance_angstrom"] = mindistance aux_results_info_entry["numatoms"] = numatoms aux_results_info.append(aux_results_info_entry) if gb_energy_J_m2 < minimum_energy: minimum_energy = gb_energy_J_m2 # Rename files for this run so that they're not deleted os.replace(energy_file, os.path.join(directory, "energy_min.out")) os.replace(numatoms_file, os.path.join(directory, "numatoms_min.out")) os.replace( mindistance_file, os.path.join(directory, "mindistance_min.out") ) os.replace( lammps_input_file, os.path.join(directory, "generated_input_min.gb") ) os.replace( lammps_output_file, os.path.join(directory, "run_output_min.gb") ) os.replace(dump_file, os.path.join(directory, "dump_post_min.cfg")) # Clean up files written by LAMMPS for fl in ( energy_file, numatoms_file, mindistance_file, lammps_input_file, lammps_output_file, dump_file, success_file, ): if os.path.isfile(fl): os.remove(fl) with open(os.path.join(directory, "aux_results_info.pkl"), "wb") as f: pickle.dump(aux_results_info, f, protocol=3) # Compute sigma value sigma = int(lambda1 * lambda2 * lambda3) while sigma % 2 == 0: sigma = round(sigma / 2) # Store results thetas.append(theta) if verbose: print("") print("minimum gb energy = {} eV".format(minimum_energy)) gb_energies.append(minimum_energy) if verbose: print("minimum distance = {} Angstroms".format(mindistance)) minimum_distances.append(mindistance) if verbose: print("sigma value = {}".format(sigma)) sigmas.append(sigma) print() return [thetas, gb_energies, minimum_distances, sigmas]