#!/usr/bin/env python3 r"""Compute the dislocation core energy of a cubic crystal. The Test Driver computes the dislocation core energy of a cubic crystal at zero temperature and a given stress state for a specified dislocation core cutoff 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) using MD++. - Second, by specifying a dislocation core radius, it will calculate the dislocation core energy. Contributing authors: Yichen Qian (Villanova University), Yaser Afshar (UMN) """ import os from os import makedirs from os.path import join, isfile, dirname, isdir import math import numpy as np import subprocess import gzip import shutil import inspect import json import kim_edn from kim_property import \ kim_property_create, \ kim_property_modify, \ kim_property_dump from ovito.io import import_file from ovito.modifiers import DislocationAnalysisModifier, \ CommonNeighborAnalysisModifier, ExpressionSelectionModifier, \ DeleteSelectedModifier, ClusterAnalysisModifier import convergence as cr PI = math.pi """Global constant, pi.""" verbose = False """Verbose mode, if True it will print useful information.""" class DislocationError(Exception): """Raise an exception. It raises an exception when receives an error message. message (str): The error message. """ def __init__(self, message): """Constuctor.""" _msg = '\nERROR(@' + \ inspect.currentframe().f_back.f_code.co_name + '): ' + message Exception.__init__(self, _msg) self.msg = _msg def __reduce__(self): """Efficient pickling.""" return self.__class__, (self.msg) def __str__(self): """Message string representation.""" return self.msg def print_warning(message): """Print a warning message. Args: msg (str): The warning message. """ _msg = '\nWARNING(@' + \ inspect.currentframe().f_back.f_code.co_name + '): ' + message print(_msg) def check_consistency_of_dislocation_parameters( burgers_vector_direction, dislocation_line_direction, slip_plane_miller_indices): """Verify that the slip plane vector is approximately parallel or anti-parallel to the cross product of the dislocation line direction and the burgers vector. In the special case that the latter cross product is zero, the slip plane can have any non-zero vector value. Args: burgers_vector_direction (1darray): Burgers vector direction dislocation_line_direction (1darray): Dislocation plane index slip_plane_miller_indices (1darray): Slip plane normal direction """ ABS_TOL = 1e-8 norm = np.linalg.norm slip_plane_norm = norm(slip_plane_miller_indices) if slip_plane_norm < ABS_TOL: raise DislocationError("Slip plane vector is approximately zero.") dislocation_cross_burgers = np.cross( dislocation_line_direction, burgers_vector_direction) dislocation_cross_burgers_norm = norm(dislocation_cross_burgers) if dislocation_cross_burgers_norm > ABS_TOL: # Check if the cross product is parallel or anti-parallel to the given # slip plane vector. dislocation_cross_burgers_cross_slip_plane_normal = np.cross( dislocation_cross_burgers, slip_plane_miller_indices ) if norm(dislocation_cross_burgers_cross_slip_plane_normal) > ABS_TOL: emsg = ("Slip plane normal is neither parallel nor anti-parallel " "to the cross product of the dislocation line direction " "and the Burgers vector.") raise DislocationError(emsg) def dxa_extract_position(pipeline, lattice): """Extract dipole positions from an ovito pipeline using DXA. Args: pipeline (:py:class:`~ovito.pipeline.Pipeline`): ovito pipeline which is created for the imported data. lattice (str): Lattice type. Returns: 2darray, 2darray: position, cvec """ # This modifier extracts all dislocations in a crystal and converts # them to continuous line segments. dxa_modifier = DislocationAnalysisModifier() dxa_modifier_lattices = { 'fcc': DislocationAnalysisModifier.Lattice.FCC, 'bcc': DislocationAnalysisModifier.Lattice.BCC, 'diamond': DislocationAnalysisModifier.Lattice.CubicDiamond } try: dxa_modifier.input_crystal_structure = \ dxa_modifier_lattices[lattice] except KeyError: emsg = 'lattice {} not supported in the DXA '.format(lattice) emsg += 'method. Valid lattices are:\n\t- ' emsg += '{}'.format('\n\t- '.join(('fcc', 'bcc', 'diamond'))) raise DislocationError(emsg) pipeline.modifiers.append(dxa_modifier) # access the data data = pipeline.compute() # Print list of dislocation lines: network = data.dislocations if verbose: print('DXA: Found {} dislocation segment(s)'.format( len(network.segments))) position = np.zeros((2, 2)) if len(network.segments) == 2: for segment in network.segments: if verbose: print('Segment {}:\n\tlength={},\n\tBurgers vector={}'.format( segment.id, segment.length, segment.true_burgers_vector)) position[segment.id, 0] = np.mean(segment.points[:, 0]) position[segment.id, 1] = np.mean(segment.points[:, 1]) if verbose: print('\taverage x,y position = {} {}'.format( position[segment.id, 0], position[segment.id, 1])) else: raise DislocationError("Failed to exactly find 2 dislocations.") cvec = data.cell.matrix[0:2, 0:2].T return position, cvec def cna_extract_position(pipeline, lattice): """Extract dipole positions from an ovito pipeline using CNA. Args: pipeline (:py:class:`~ovito.pipeline.Pipeline`): ovito pipeline which is created for the imported data. lattice (str): Lattice type. Returns: 2darray, 2darray: position, cvec """ # This modifier analyzes the local neighborhood of each particle to # identify simple crystalline structures. cna_modifier = CommonNeighborAnalysisModifier() # The structure types recognized by the common neighbor analysis are # FCC, HCP, BCC, and ICO (Icosahedral structure) cna_modifier_lattices = { 'fcc': 'StructureType==1', # 'hcp': 'StructureType==2', 'bcc': 'StructureType==3', # 'ico': 'StructureType==4', } try: select_perfect = ExpressionSelectionModifier( expression=cna_modifier_lattices[lattice]) except KeyError: emsg = 'lattice {} not supported in the CNA '.format(lattice) emsg += 'method. Valid lattices structure types recognized ' emsg += 'by CNA are:\n\t- ' emsg += '{}'.format('\n\t- '.join(('fcc', 'bcc'))) raise DislocationError(emsg) delete_selected = DeleteSelectedModifier() # This modifier groups particles into disconnected clusters based on a # selectable connectivity criterion. cla_modifier = ClusterAnalysisModifier() pipeline.modifiers.append(cna_modifier) pipeline.modifiers.append(select_perfect) pipeline.modifiers.append(delete_selected) pipeline.modifiers.append(cla_modifier) # access the data data = pipeline.compute() cluster_count = data.attributes['ClusterAnalysis.cluster_count'] if verbose: print('CNA: Found {} dislocation segment(s)'.format(cluster_count)) position = np.zeros((2, 2)) if cluster_count == 2: positions = data.particles['Position'] idx = np.digitize(data.particles['Cluster'] - 1, bins=[0, 1]) - 1 for i in range(2): if verbose: print('Segment {}:'.format(i)) # Fold x positions wrt first point in cluster posx = positions[idx == i, 0] lenx = data.cell.matrix[0, 0] posx = posx - np.rint((posx - posx[0]) * 1. / lenx) * lenx position[i, 0] = np.mean(posx) position[i, 1] = np.mean(positions[idx == i, 1]) if verbose: print('\taverage x,y position = {} {}'.format( position[i, 0], position[i, 1])) else: raise DislocationError("Failed to exactly find 2 dislocations.") cvec = data.cell.matrix[0:2, 0:2].T return position, cvec def extract_position(filename, *, method='cna', lattice='bcc'): """Extract dipole positions from a file. Args: filename (str): File name method (str): `dxa`, or `cna` (default: 'cna') CNA: common neighbor analysis DXA: dislocation extraction algorithm. Returns: 2darray, 2darray: position, cvec """ if not isfile(filename): emsg = "The file: \n{}\n does not exist!".format(filename) raise DislocationError(emsg) if method not in ('dxa', 'cna'): emsg = 'method \'{}\' not found. Valid methods '.format(method) emsg += 'are:\n\t- {}'.format('\n\t- '.join(('dxa', 'cna'))) raise DislocationError(emsg) ovito_pipeline = import_file(filename) # Extract dipole position using DXA if method == 'dxa': try: position, cvec = dxa_extract_position(ovito_pipeline, lattice) except: emsg = "Failed to extract dipole position using DXA." raise DislocationError(emsg) # Extract dipole position using CNA if method == 'cna': try: position, cvec = cna_extract_position(ovito_pipeline, lattice) except: emsg = "Failed to extract dipole position using CNA." raise DislocationError(emsg) return position, cvec def run_lammps(input_file, output_file=None, lammps_exe="lmp", num_procs=2): """Run lammps with a given input file and write the output to the output. Args: input_file (str): Input file output_file (str): Output file (Default is None) lammps_exe (str): LAMMPS executable (Default is "lmp") num_procs (int): number of processors (Default is 1) """ if output_file is None: if isinstance(num_procs, int) and num_procs > 1: try: subprocess.check_call(["mpirun", "-np", str(num_procs), lammps_exe, "-in", input_file], shell=False) except subprocess.CalledProcessError: try: extra_info = open("log.lammps").read() except IOError: extra_info = "no log file" emsg = "LAMMPS did not exit properly:\n" + extra_info raise DislocationError(emsg) return try: subprocess.check_call([lammps_exe, "-in", input_file], shell=False) except subprocess.CalledProcessError: try: extra_info = open("log.lammps").read() except IOError: extra_info = "no log file" emsg = "LAMMPS did not exit properly:\n" + extra_info raise DislocationError(emsg) return with open(output_file, "w") as f_out: if isinstance(num_procs, int) and num_procs > 1: try: subprocess.check_call(["mpirun", "-np", str(num_procs), lammps_exe, "-in", input_file], shell=False, stdout=f_out) except subprocess.CalledProcessError: try: extra_info = open("log.lammps").read() except IOError: extra_info = "no log file" emsg = "LAMMPS did not exit properly:\n" + extra_info raise DislocationError(emsg) return try: subprocess.check_call([lammps_exe, "-in", input_file], shell=False, stdout=f_out) except subprocess.CalledProcessError: try: extra_info = open("log.lammps").read() except IOError: extra_info = "no log file" emsg = "LAMMPS did not exit properly:\n" + extra_info raise DislocationError(emsg) # Define the lattice type by full name(will use in generating dipole structure) def get_lattice_type_full_name(lattice): """Get the lattice type by its name based on the input lattice type. Args: lattice (str): Lattice type. One of the "bcc", "fcc", "diamond", and "sc" lattice type. Returns: str: Lattice full name """ if lattice not in ('fcc', 'bcc', 'diamond'): emsg = 'lattice {} is not supported.'.format(lattice) emsg += 'Valid lattices are:\n\t- ' emsg += '{}'.format('\n\t- '.join(('fcc', 'bcc', 'diamond'))) raise DislocationError(emsg) full_names = {"bcc": "body-centered-cubic", "fcc": "face-centered-cubic", "diamond": "diamond-cubic", "sc": "simple-cubic"} return full_names[lattice] def get_basis_atom_coordinates(lattice): """Get the basis atom coordinates based on the input lattice type. Args: lattice (str): Lattice type. One of the "bcc", "fcc", "diamond", and "sc" lattice type. Returns: ndarray: Basis atom coordinates of the lattice type. """ if lattice not in ('fcc', 'bcc', 'diamond'): emsg = 'lattice {} is not supported.'.format(lattice) emsg += 'Valid lattices are:\n\t- ' emsg += '{}'.format('\n\t- '.join(('fcc', 'bcc', 'diamond'))) raise DislocationError(emsg) 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] def format_basis_atom_coordinates(basis_atom_coordinates): """Create a kim_property formatted of atom coordinates. Takes a nested list of basis atom coordinates and creates a list that can be splatted into a kim_property modify command Args: basis_atom_coordinates (str): Basis atom coordinates Returns: str: kim_property formatted string of atom coordinates. """ output = [] for i, atom in enumerate(basis_atom_coordinates): output += ["source-value", str(i + 1), "1:3", *list(map(str, atom))] return output def format_species(species_name, lattice): """Create a kim_property format of species based on the input 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 Args: species_name (str): Species name lattice (str): Lattice type. One of the "bcc", "fcc", "diamond", and "sc" lattice type. Returns: str: kim_property formatted of species. """ if lattice not in ('fcc', 'bcc', 'diamond'): emsg = 'lattice {} is not supported.'.format(lattice) emsg += 'Valid lattices are:\n\t- ' emsg += '{}'.format('\n\t- '.join(('fcc', 'bcc', 'diamond'))) raise DislocationError(emsg) num_basis_atoms = {"fcc": 4, "bcc": 2, "diamond": 8, "sc": 1} n_atoms = num_basis_atoms[lattice] formatted_species = ["source-value", "1:{}".format(n_atoms), *[species_name] * n_atoms] return formatted_species def get_space_group(lattice): """Get the space group based on the input lattice type. Args: lattice (str): Lattice type. One of the "bcc", "fcc", "diamond", and "sc" lattice type. Returns: str: Space group based on the input lattice type. """ if lattice not in ('fcc', 'bcc', 'diamond'): emsg = 'lattice {} is not supported.'.format(lattice) emsg += 'Valid lattices are:\n\t- ' emsg += '{}'.format('\n\t- '.join(('fcc', 'bcc', 'diamond'))) raise DislocationError(emsg) space_groups = {"fcc": "Fm-3m", "bcc": "Im-3m", "diamond": "Fd-3m", "sc": "Pm-3m", } return space_groups[lattice] def madsum_nonsingular(mu, nu, b, rc, dcoord, trunc): r"""Calculate the elastic energy of a dislocation dipole in a supercell under periodic boundary conditions. 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) Args: mu (float): Shear modulus (in eV/A^3, 1eV/A^3 = 160.22GPa). nu (float): Poisson's ratio. b (1darray): [ bx, by, bz ] Burgers vector rc (float): Cut-off radius. dcoord (1darray): [ 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 (1darray): [ xcut, ycut ] truncation paramters -xcut:xcut times -ycut:ycut number of images will be added Returns: float, float, float: Eel, Eprm, Eimg 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 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) Contributors: Wei Cai, caiwei@stanford.edu """ x0 = dcoord[0] x1 = dcoord[1] y1 = dcoord[2] x2r = dcoord[3] y2r = dcoord[4] xcut = trunc[0] ycut = trunc[1] c1 = np.array([x0, 0.0]) c2 = np.array([x1, y1]) ds = np.array([x2r, y2r]) dr = np.matmul(np.vstack((c1, c2)).T, ds) R = np.linalg.norm(dr) Ra2 = R * R + rc * rc Ra = np.sqrt(Ra2) dt = dr / Ra prm = np.array([np.log(Ra / rc), 1.0 / Ra2, 1.0, dt[0] * dt[0], dt[1] * dt[1], dt[0] * dt[1]]) Efac = -mu / 4.0 / PI * np.array([ b[2] * b[2] + 1. / (1. - nu) * (b[0] * b[0] + b[1] * b[1]), -0.5 * rc * rc * b[2] * b[2], 0.5 * b[2] * b[2], 1. / (1. - nu) * b[1] * b[1], 1. / (1. - nu) * b[0] * b[0], -1. / (1. - nu) * 2. * b[0] * b[1] ]) Eprm = np.sum(Efac * prm) * (-2.0) # Add image contribution (naive summation) # image coordinate (x,y) image strength image = np.array([ [0., 0., 2.], [dr[0], dr[1], -1.0], [-dr[0], -dr[1], -1.0] ]) xc = np.arange(-xcut, xcut + 1) yc = np.arange(-ycut, ycut + 1) Xim, Yim = np.meshgrid(xc, yc) Xim = np.tile(Xim.reshape(xc.size, yc.size, 1), (1, 1, 2)) Yim = np.tile(Yim.reshape(xc.size, yc.size, 1), (1, 1, 2)) offset = Xim * c1 + Yim * c2 offset = offset.reshape(-1, 2) Esum = np.zeros(6) offset0 = offset[(offset[:, 0] != 0) | (offset[:, 1] != 0)] for k in range(image.shape[0]): dR = offset0 + image[k, 0:2] R = np.sqrt(np.sum(dR * dR, axis=1, keepdims=True)) Ra2 = R * R + rc * rc Ra = np.sqrt(Ra2) dt = dR/np.tile(Ra, (1, 2)) Esum += image[k, 2] * np.sum(np.hstack([ np.log(Ra / rc), 1. / Ra2, np.ones(Ra.shape), dt * dt, np.prod(dt, axis=1, keepdims=True) ]), axis=0) # Add ghost contribution (correction for conditional convergence) # ghost coordinate (x,y) ghost strength 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.zeros(6) for k in range(ghost.shape[0]): dR = offset + ghost[k, 0:2] R = np.sqrt(np.sum(dR * dR, axis=1, keepdims=True)) Ra2 = R * R + rc * rc Ra = np.sqrt(Ra2) dt = dR / np.tile(Ra, (1, 2)) Egst -= ghost[k, 2] * np.sum(np.hstack([ np.log(Ra / rc), 1. / Ra2, np.ones(Ra.shape), dt * dt, np.prod(dt, axis=1, keepdims=True) ]), axis=0) Esum = Esum + Egst Eimg = np.sum(Efac * Esum) Eel = Eprm + Eimg return Eel, Eprm, Eimg def madsum_iso(mu, nu, b, rc, dcoord, trunc): r"""Calculate the elastic energy of a dislocation dipole in a supercell under periodic boundary conditions This program 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) Args: mu (float): Shear modulus (in eV/A^3, 1eV/A^3 = 160.22GPa). nu (float): Poisson's ratio. b (1darray): [ bx, by, bz ] Burgers vector rc (float): Cut-off radius. dcoord (1darray): [ 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 (1darray): [ xcut, ycut ] truncation paramters -xcut:xcut times -ycut:ycut number of images will be added Returns: float, float, float: Eel, Eprm, Eimg 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 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) Contributors: Nicolas Bertin, bertin1@llnl.gov based on the implemetation by Wei Cai, caiwei@stanford.edu """ x0 = dcoord[0] x1 = dcoord[1] y1 = dcoord[2] x2r = dcoord[3] y2r = dcoord[4] xcut = trunc[0] ycut = trunc[1] c1 = np.array([x0, 0.]) c2 = np.array([x1, y1]) ds = np.array([x2r, y2r]) dr = np.matmul(np.vstack((c1, c2)).T, 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.0/PI)*np.array([ (b[0]**2 + b[1]**2) / (1. - nu) + b[2]**2, b[1]**2 / (1. - nu), b[0]**2 / (1. - nu), -2. * b[0] * b[1] / (1. - nu) ]) Eprm = np.sum(Efac * prm) * (-2.0) # Add image contribution (naive summation) # image coordinate (x,y) image strength image = np.array([ [0., 0., 2.], [dr[0], dr[1], -1.], [-dr[0], -dr[1], -1.] ]) xc = np.arange(-xcut, xcut+1) yc = np.arange(-ycut, ycut+1) Xim, Yim = np.meshgrid(xc, yc) Xim = np.tile(Xim.reshape(xc.size, yc.size, 1), (1, 1, 2)) Yim = np.tile(Yim.reshape(xc.size, yc.size, 1), (1, 1, 2)) offset = Xim*c1 + Yim*c2 offset = offset.reshape(-1, 2) Esum = np.zeros(4,) offset0 = offset[(offset[:, 0] != 0) | (offset[:, 1] != 0)] for k in range(image.shape[0]): dR = offset0 + image[k, 0:2] R = np.sqrt(np.sum(dR**2, axis=1, keepdims=True)) dt = dR/np.tile(R, (1, 2)) Esum += image[k, 2] * np.sum(np.hstack(( np.log(R / rc), dt**2, np.prod(dt, axis=1, keepdims=True) )), axis=0) # Add ghost contribution (correction for conditional convergence) # ghost coordinate (x,y) ghost strength 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.zeros(4,) for k in range(ghost.shape[0]): dR = offset + ghost[k, 0:2] R = np.sqrt(np.sum(dR**2, axis=1, keepdims=True)) dt = dR / np.tile(R, (1, 2)) Egst -= ghost[k, 2] * np.sum(np.hstack(( np.log(R / rc), dt**2, np.prod(dt, axis=1, keepdims=True) )), axis=0) Esum = Esum + Egst Eimg = np.sum(Efac * Esum) Eel = Eprm + Eimg return Eel, Eprm, Eimg def sextic(b, C): """Compute the classical sextic formalism of anisotropic linear elasticity for dislocations [p,A,B,D]=sextic(b,C) Displacement field of dislocation (along z direction) u(k)=Re(-1/(2*pi*i)*Sum_{n=1}^{3}(A(k,n)*D(n)*ln(eta(n)))), where, eta(n)=x+p(n)*y Args: b (1darray): Burger's vector C (ndarray): Elastic Constant Matrix C(3,3,3,3) Returns: ndarray : p, A, B, D, Poly6, a, c0, c1, c2 p(6): Six roots of the sextic polynomial p(n) A(3,3): A(k,n) B{3}(3,3,3): B{n}(i,j,k) D(3): D(n) """ # Construct Poly6 c0 = C[:, 0, :, 0] c1 = C[:, 0, :, 1] + C[:, 1, :, 0] c2 = C[:, 1, :, 1] Poly6 = np.zeros((7,)) Poly6[6] = np.linalg.det(c0) Poly6[5] = \ np.linalg.det(np.vstack([c0[:, 0], c0[:, 1], c1[:, 2]]).T) + \ np.linalg.det(np.vstack([c0[:, 0], c1[:, 1], c0[:, 2]]).T) + \ np.linalg.det(np.vstack([c1[:, 0], c0[:, 1], c0[:, 2]]).T) Poly6[4] = \ np.linalg.det(np.vstack([c0[:, 0], c0[:, 1], c2[:, 2]]).T) + \ np.linalg.det(np.vstack([c0[:, 0], c2[:, 1], c0[:, 2]]).T) + \ np.linalg.det(np.vstack([c2[:, 0], c0[:, 1], c0[:, 2]]).T) + \ np.linalg.det(np.vstack([c0[:, 0], c1[:, 1], c1[:, 2]]).T) + \ np.linalg.det(np.vstack([c1[:, 0], c0[:, 1], c1[:, 2]]).T) + \ np.linalg.det(np.vstack([c1[:, 0], c1[:, 1], c0[:, 2]]).T) Poly6[3] = \ np.linalg.det(np.vstack([c0[:, 0], c1[:, 1], c2[:, 2]]).T) + \ np.linalg.det(np.vstack([c1[:, 0], c0[:, 1], c2[:, 2]]).T) + \ np.linalg.det(np.vstack([c0[:, 0], c2[:, 1], c1[:, 2]]).T) + \ np.linalg.det(np.vstack([c1[:, 0], c2[:, 1], c0[:, 2]]).T) + \ np.linalg.det(np.vstack([c2[:, 0], c0[:, 1], c1[:, 2]]).T) + \ np.linalg.det(np.vstack([c2[:, 0], c1[:, 1], c0[:, 2]]).T) + \ np.linalg.det(np.vstack([c1[:, 0], c1[:, 1], c1[:, 2]]).T) Poly6[2] = \ np.linalg.det(np.vstack([c1[:, 0], c1[:, 1], c2[:, 2]]).T) + \ np.linalg.det(np.vstack([c1[:, 0], c2[:, 1], c1[:, 2]]).T) + \ np.linalg.det(np.vstack([c2[:, 0], c1[:, 1], c1[:, 2]]).T) + \ np.linalg.det(np.vstack([c0[:, 0], c2[:, 1], c2[:, 2]]).T) + \ np.linalg.det(np.vstack([c2[:, 0], c0[:, 1], c2[:, 2]]).T) + \ np.linalg.det(np.vstack([c2[:, 0], c2[:, 1], c0[:, 2]]).T) Poly6[1] = \ np.linalg.det(np.vstack([c1[:, 0], c2[:, 1], c2[:, 2]]).T) + \ np.linalg.det(np.vstack([c2[:, 0], c1[:, 1], c2[:, 2]]).T) + \ np.linalg.det(np.vstack([c2[:, 0], c2[:, 1], c1[:, 2]]).T) Poly6[0] = np.linalg.det(c2) # Solve the Polynomial poly6_roots = np.roots(Poly6) # Sort roots ind = np.argsort(-np.imag(poly6_roots)) p = poly6_roots p[0:3] = poly6_roots[ind[0:3]] p[3:6] = np.conj(p[0:3]) # Construct Matrix A a = [] if (np.abs(p[0] - 1j) < 1e-4) & \ (np.abs(p[1] - 1j) < 1e-4) & (np.abs(p[2] - 1j) < 1e-4): if verbose: print('degeneracy!') for n in range(3): a.append(c0 + c1 * p[n] + c2 * p[n]**2) A = np.array([[1., -1j, 0.], [1j, 1., 0.], [0., 0., 1.]]) else: A = np.zeros((3, 3), dtype=complex) for n in range(3): a.append(c0 + c1 * p[n] + c2 * p[n]**2) w, v = np.linalg.eig(a[n]) t = np.argsort(np.abs(w)) A[:, n] = v[:, t[0]] # Construct Matrix B B = [] for n in range(3): B.append(C[:, :, :, 0] + C[:, :, :, 1] * p[n]) F = np.zeros((3, 3), dtype=complex) for j in range(3): for n in range(3): F[j, n] = B[n][j, 1, 0] * A[0, n] + \ B[n][j, 1, 1] * A[1, n] + B[n][j, 1, 2] * A[2, n] G = np.zeros((6, 6)) G[0:3, 0:3] = np.real(A) G[0:3, 3:6] = -np.imag(A) G[3:6, 0:3] = np.real(F) G[3:6, 3:6] = -np.imag(F) rhs = np.zeros((6,)) rhs[0:3] = b lhs = np.linalg.solve(G, rhs) D = lhs[0:3] + 1j * lhs[3:6] return p, A, B, D, Poly6, a, c0, c1, c2 def madsum_aniso(C0, b0, M, rc, dcoord, trunc): r"""Calculate the elastic energy of a dislocation dipole in a supercell under periodic boundary conditions. This program 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) Args: C0 (1darray): C(i,j,k,l) general elastic constant tensor (in eV/A^3, 1eV/A^3 = 160.22GPa). b0 (1darray): [bx, by, bz] Burgers vector in unit cell frame (in Angstrom), e.g. b0 = [ 1 1 1 ]/2 * 3.1472 (for Molybdenum). M (ndarray): [e1, e2, e3] coordinate transformation matrix, e3 is parallel to dislocation line. rc (float): Cut-off radius. dcoord (1darray): [ 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 (1darray): [ xcut, ycut ] truncation paramters -xcut:xcut times -ycut:ycut number of images will be added Returns: float, float, float: Eel, Eprm, Eimg 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 References: [1] W. Cai, Atomistic and Mesoscale Modeling of Dislocation Mobility, Ph. D. Thesis, Massachusetts Institute of Technology, 2001. (http://asm.mit.edu/caiwei/Download/Thesis/CaiWeiThesis.pdf p.296, for dislocation interaction energy in anisotropic medium.) [2] 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) Contributors: Nicolas Bertin, bertin1@llnl.gov based on the implemetation by Wei Cai, caiwei@stanford.edu """ # Step 1: Solve the energy prefactor in anisotropic elasticity # Transform C Mt = M.T C = np.zeros((3, 3, 3, 3)) for i in range(3): for j in range(3): for k in range(3): for l in range(3): tmp = 0.0 for ip in range(3): for jp in range(3): for kp in range(3): for lp in range(3): tmp += Mt[i, ip] *\ Mt[j, jp] *\ Mt[k, kp] *\ Mt[l, lp] *\ C0[ip, jp, kp, lp] C[i, j, k, l] = tmp # Transform b b = np.matmul(Mt, b0) # Solve the sextic equation in the coordinate specified by M p, A, B, D, poly6, a, c0, c1, c2 = sextic(b, C) del(poly6, a, c0) # Calculate interaction matrix h1b(i,n) # only p and h1b are used in following calculations h1 = np.zeros((3, 3), dtype=complex) for n in range(3): for i in range(3): h1[i, n] = B[n][i, 1, 0] * A[0, n] * D[n]+B[n][i, 1, 1] * \ A[1, n] * D[n]+B[n][i, 1, 2] * A[2, n] * D[n] h1b = np.matmul(h1.T, b) # Step 2: image summation # # (Ref. 1, p. 296) # W(x,y) = sum_{n=1}^{3} Re[ h1b(n) / (2 * pi * I) * ln (x+p_n * y)/rc ] x0 = dcoord[0] x1 = dcoord[1] y1 = dcoord[2] x2r = dcoord[3] y2r = dcoord[4] xcut = trunc[0] ycut = trunc[1] c1 = np.array([x0, 0.]) c2 = np.array([x1, y1]) ds = np.array([x2r, y2r]) dr = np.matmul(np.vstack((c1, c2)).T, ds) prm = np.log((dr[0]+dr[1] * p[0:3])/rc) # New edit # p= p[::-1] # prm = np.log((dr[0] * p[0:3]+dr[1])/rc) Efac = h1b/(2. * PI * 1j) * (-0.5) Eprm = np.sum(np.real(Efac * prm)) * (-2.0) # Add image contribution (naive summation) # image coordinate (x,y) image strength image = np.array([[0., 0., 2.], [dr[0], dr[1], -1.], [-dr[0], -dr[1], -1.]]) xc = np.arange(-xcut, xcut+1) yc = np.arange(-ycut, ycut+1) Xim, Yim = np.meshgrid(xc, yc) Xim = np.tile(Xim.reshape(xc.size, yc.size, 1), (1, 1, 2)) Yim = np.tile(Yim.reshape(xc.size, yc.size, 1), (1, 1, 2)) offset = Xim * c1 + Yim * c2 offset = offset.reshape(-1, 2) Esum = np.zeros((3,), dtype=complex) offset0 = offset[(offset[:, 0] != 0) | (offset[:, 1] != 0)] for k in range(image.shape[0]): dR = offset0 + image[k, 0:2] tmp = np.log((np.tile(dR[:, 0], (3, 1)).T + np.outer(dR[:, 1], p[0:3])) / rc) Esum += image[k, 2] * np.sum(tmp, axis=0) # Add ghost contribution (correction for conditional convergence) # ghost coordinate (x,y) ghost strength 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.zeros((3,), dtype=complex) for k in range(ghost.shape[0]): dR = offset + ghost[k, 0:2] tmp = np.log((np.tile(dR[:, 0], (3, 1)).T + np.outer(dR[:, 1], p[0:3])) / rc) Egst -= ghost[k, 2] * np.sum(tmp, axis=0) Esum = Esum + Egst Eimg = np.sum(np.real(Efac * Esum)) Eel = Eprm + Eimg return Eel, Eprm, Eimg def madsum_cubic(C, b0, M, rc, dcoord, trunc): r"""Calculate the elastic energy of a dislocation dipole in a supercell under periodic boundary conditions This program 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) Args: C (1darray): [c11, c12, c44] cubic elastic constants (in eV/A^3, 1eV/A^3 = 160.22GPa). b0 (1darray): [bx, by, bz] Burgers vector in unit cell frame (in Angstrom), e.g. b0 = [ 1 1 1 ]/2 * 3.1472 (for Molybdenum). M (ndarray): [e1, e2, e3] coordinate transformation matrix, e3 is parallel to dislocation line. rc (float): Cut-off radius. dcoord (1darray): [ 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 (1darray): [ xcut, ycut ] truncation paramters -xcut:xcut times -ycut:ycut number of images will be added Returns: float, float, float: Eel, Eprm, Eimg 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 References: [1] W. Cai, Atomistic and Mesoscale Modeling of Dislocation Mobility, Ph. D. Thesis, Massachusetts Institute of Technology, 2001. (http://asm.mit.edu/caiwei/Download/Thesis/CaiWeiThesis.pdf p.296, for dislocation interaction energy in anisotropic medium.) [2] 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) Contributors: Nicolas Bertin, bertin1@llnl.gov based on the implementation by Wei Cai, caiwei@stanford.edu """ c11 = C[0] c12 = C[1] c44 = C[2] # Voigt elasticity matrix C66 = np.array([[c11, c12, c12, 0.0, 0.0, 0.0], [c12, c11, c12, 0.0, 0.0, 0.0], [c12, c12, c11, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, c44, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, c44, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, c44]]) # Voigt notation cind = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]]) C0 = np.zeros((3, 3, 3, 3)) for i in range(3): for j in range(3): for k in range(3): for l in range(3): C0[i, j, k, l] = C66[cind[i, j], cind[k, l]] return madsum_aniso(C0, b0, M, rc, dcoord, trunc) def generate_dipole_parameters(burgers_vector, dislocation_line_direction, slip_plane_direction): """Generate a supercell parameters. Args: burgers_vector (1darray): Burgers vector dislocation_line_direction (1darray): Dislocation plane index slip_plane_direction (1darray): Slip plane direction Returns: 1darray, 1darray, 1darray, float: cm1, cm2, cm3, angle x direction, dipole direction, dipole line direction, and supercell angle """ # Burgers vector direction bv = np.array(burgers_vector, copy=False) bv_norm = np.linalg.norm(bv) # dislocation plane index c3 = np.array(dislocation_line_direction, copy=True) c2 = np.array(slip_plane_direction, copy=False) cm2 = np.array(c2, copy=False) y0 = np.cross(c2, bv) y0 = y0.astype(int) my = np.gcd(y0[0], np.gcd(y0[1], y0[2])) y0 = y0 / my y0_norm = np.linalg.norm(y0) x = bv / bv_norm y = y0 / y0_norm # maximum Miller index of repeat vectors allowed to # generate supercells of various character angles c3x = np.dot(c3, x) c3y = np.dot(c3, y) angle = np.arctan2(c3y, c3x) * 180.0 / PI # determine supercell size # compute complementary supercell repeat vector c1 = np.cross(c3, c2) c1 = c1.astype(int) m1 = np.gcd(c1[0], np.gcd(c1[1], c1[2])) cm1 = c1 / m1 c3 = c3.astype(int) m3 = np.gcd(c3[0], np.gcd(c3[1], c3[2])) cm3 = c3 / m3 return cm1, cm2, cm3, angle def parse_coeffs(line): """Parse a line of coeffs. Args: line (str): input string. Returns: 1darray: Array of values """ words = "".join((ch if ch in "0123456789.-e" else " ") for ch in line) try: coeffs = [float(word) for word in words.split()] except: msg = "The input line contains a non floating-point character." raise DislocationError(msg) return coeffs def dump_property(property_name, lattice, species, lattice_constant, slip_plane_miller_indices, dislocation_line_direction, burgers_vector_direction, dislocation_core_radius_, core_energy_nonsingular, core_energy_nonsingular_std, core_energy_isotropic, core_energy_isotropic_std, core_energy_anisotropic, core_energy_anisotropic_std): """Create and dump the KIM property instance. Create a KIM property instance and take the generated instance validate against the property definition and writes the instance out to a file in edn format. Args: property_name {str} -- - A string containing the property name or - unique ID of the property, or - a path-like object giving the pathname (absolute or relative to the current working directory) of the file to be opened lattice (str): Lattice type species (str): Species name lattice_constant (float): Equilibrium lattice constant slip_plane_miller_indices (1darray): Slip plane direction dislocation_line_direction (1darray): Dislocation plane index burgers_vector_direction (1darray): Burgers vector direction dislocation_core_radius_ (1darray): Dislocation core radii core_energy_nonsingular (1darray): [description] core_energy_nonsingular_std (1darray): [description] core_energy_isotropic (1darray): [description] core_energy_isotropic_std (1darray): [description] core_energy_anisotropic (1darray): [description] core_energy_anisotropic_std (1darray): [description] """ # Create a property instance prop_inst = kim_property_create(1, property_name) num_dislocation_core_radii = len(dislocation_core_radius_) if num_dislocation_core_radii < 1: raise DislocationError("There is no dislocation core radius.") if len(core_energy_nonsingular) != num_dislocation_core_radii or \ len(core_energy_nonsingular_std) != num_dislocation_core_radii or \ len(core_energy_isotropic) != num_dislocation_core_radii or \ len(core_energy_isotropic_std) != num_dislocation_core_radii or \ len(core_energy_anisotropic) != num_dislocation_core_radii or \ len(core_energy_anisotropic_std) != num_dislocation_core_radii: emsg = "The input length do not match with each other." raise DislocationError(emsg) ndcr_str = "1:{}".format(num_dislocation_core_radii) # 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, "key", "species", *format_species(species, lattice), "key", "a", "source-value", lattice_constant, "source-unit", "angstrom", "key", "basis-atom-coordinates", *format_basis_atom_coordinates(get_basis_atom_coordinates(lattice)), "key", "space-group", "source-value", get_space_group(lattice), "key", "cauchy-stress", "source-value", "1:6", *[0.0] * 6, "source-unit", "GPa", "key", "slip-plane-miller-indices", "source-value", "1:3", *slip_plane_miller_indices, "key", "dislocation-line-direction", "source-value", "1:3", *dislocation_line_direction, "key", "burgers-vector-direction", "source-value", "1:3", *burgers_vector_direction, "key", "dislocation-core-radius", "source-value", ndcr_str, *dislocation_core_radius_, "key", "core-energy-nonsingular", "source-value", ndcr_str, *core_energy_nonsingular, "source-std-uncert-value", ndcr_str, *core_energy_nonsingular_std, "source-unit", "eV/A", "key", "core-energy-isotropic", "source-value", ndcr_str, *core_energy_isotropic, "source-std-uncert-value", ndcr_str, *core_energy_isotropic_std, "source-unit", "eV/A", "key", "core-energy-anisotropic", "source-value", ndcr_str, *core_energy_anisotropic, "source-std-uncert-value", ndcr_str, *core_energy_anisotropic_std, "source-unit", "eV/A", ) # Dump the results in a file with open("output/results.edn", "w") as fp: kim_property_dump(prop_inst, fp) def get_shear_modulus_and_poisson_ratio(elastic_constants): """Get the shear modulus and Poisson's ratio. The shear modulus describes the material's response to shear stress. The Poisson's ratio describes the response in the directions orthogonal to uniaxial stress. Args: elastic_constants (array like): [c11, c12, c44] The isothermal elastic constants of a cubic crystal array containing three floats. Returns: float, float: shear_modulus, poisson_ratio """ c11 = elastic_constants[0] c12 = elastic_constants[1] poisson_ratio = c12 / (c11 + c12) Ey = c11 - 2.0 * c12**2 / (c11 + c12) shear_modulus = Ey / 2.0 / (1.0 + poisson_ratio) return shear_modulus, poisson_ratio def find_burgers_coord(angle, x_direction, x_length, dislocation_line_direction, z_length): """Compute Burgers vector in scaled coordinates. Args: angle (float): supercell angle x_direction (1darray): x direction vector x_length (float): supercell size along the x-direction dislocation_line_direction (1darray): dipole line direction z_length (float): supercell size along the z-direction Returns: 1darray: Burgers vector in scaled coordinates """ # compute Burgers vector in scaled coordinates x_vec = x_length * x_direction z_vec = z_length * dislocation_line_direction det = x_vec[0] * z_vec[1] - x_vec[1] * z_vec[0] if np.abs(det) > 1e20: a = 0.5 * (z_vec[1] - z_vec[0]) / det b = 0.5 * (x_vec[0] - x_vec[1]) / det else: det = x_vec[2] * z_vec[1] - x_vec[1] * z_vec[2] a = 0.5 * (z_vec[1] - z_vec[2]) / det b = 0.5 * (x_vec[2] - x_vec[1]) / det # make sure always remove atom, not insert atom if angle < 0: a = -a b = -b return np.array([a, 0.0, b], dtype=np.float64) def get_supercell_info(x_direction, dipole_direction, dislocation_line_direction, angle, mult, aspect_ratio, x_length, z_length): """Get the supercell size in x, y, z directions and Burgers vector. Args: x_direction (1darray): x direction vector dipole_direction (1darray): dipole direction dislocation_line_direction (1darray): dislocation line direction angle (float): supercell angle mult (float): multiplication factor aspect_ratio (float): aspect ratio x/y x_length (float): supercell size along the x-direction z_length (float): supercell size along the z-direction Returns: float, float, float, 1darray: x_length, y_length, z_length, bs supercell size in x, y, z directions and Burgers vector in scaled coordinates """ l1 = np.linalg.norm(x_direction) l2 = np.linalg.norm(dipole_direction) l3 = np.linalg.norm(dislocation_line_direction) x_length = np.ceil(mult * x_length / l1) z_length = np.ceil(mult * z_length / l3) y_length = round(aspect_ratio * x_length * l1 / l2) burgers_vector_in_scaled_coordinates = \ find_burgers_coord(angle, x_direction, x_length, dislocation_line_direction, z_length) return x_length, y_length, z_length, burgers_vector_in_scaled_coordinates def get_number_of_processors_to_use(mult, max_procs): """Get the number of processors to use based on the multiplication factor. Args: mult (int or float): multiplication factor max_procs (int): maximum number of processors Returns: int: number of processors to use """ # make a guess on the number of processors to use mult = round(mult) if mult in range(5): return min(2, max_procs) elif mult in range(5, 9): return min(4, max_procs) elif mult in range(9, 11): return min(8, max_procs) elif mult in range(11, 13): return min(16, max_procs) return min(32, max_procs) class SuperCell: """Seuper cell class.""" def __init__(self, species, mass, kim_model, lattice, lattice_constant, cohesive_energy, max_procs): """Initialize the Seuper cell object. Args: species (str): species mass (str): masses in g/mol (must match the order of species above) kim_model (str): KIM model name lattice (str): lattice type lattice_constant (float): equilibrium lattice constant max_procs (int): maximum number of processors """ self.species = species self.mass = mass self.kim_model = kim_model self.lattice_type = lattice self.lattice_constant = lattice_constant self.cohesive_energy = cohesive_energy self.max_procs = max_procs self.init() def init(self): """Initialize.""" # Current directory template_directory = dirname(__file__) # Fisr step minimization try: self.min_cg_1_template = \ open(join(template_directory, "in.min_cg_1.template")).read() except: emsg = "Failed to read the LAMMPS input template file for the " emsg += "1st relaxation." raise DislocationError(emsg) try: self.min_cg_npt_template = \ open(join(template_directory, "in.min_cg_npt.template")).read() except: emsg = "Failed to read the LAMMPS input template file." raise DislocationError(emsg) try: self.dislocation_dipole_template = \ open(join(template_directory, "dislocation_dipole.template")).read() except: emsg = "Failed to read the dislocation dipole template file." raise DislocationError(emsg) self.lammps_input_directory = join("output", "lammps", "in") self.lammps_output_directory = join("output", "lammps", "out") self.dumping_directory = join("output", "dump") # Ensure the directories we need are created try: makedirs(self.lammps_input_directory) except OSError: pass if not isdir(self.lammps_input_directory): emsg = "Failed to create a lammps input directory." raise DislocationError(emsg) try: makedirs(self.lammps_output_directory) except OSError: pass if not isdir(self.lammps_output_directory): emsg = "Failed to create a lammps output directory." raise DislocationError(emsg) try: makedirs(self.dumping_directory) except OSError: pass if not isdir(self.dumping_directory): raise DislocationError("Failed to create a dumping directory.") def create_dipole(self, n1, n2, n3, bs, x_direction, dipole_direction, dislocation_line_direction, mult, nu): """Crate dipole input using MD++. Args: n1 (float): supercell size in the x direction n2 (float): supercell size in the y direction n3 (float): supercell size in the z direction bs (1darray): Burgers vector in scaled coordinates dipole_direction (1darray): dipole direction x_direction (1darray): x direction dislocation_line_direction (1darray): dislocation line direction mult (float): multiplication factor nu (float): poisson ratio """ dipole_creation_file = \ join(self.lammps_input_directory, "dislocation_dipole_creation_{}.mdpp".format(round(mult, 3))) # set necessary input parameters in input files lattice_full_name = get_lattice_type_full_name(self.lattice_type) with open(dipole_creation_file, "w") as f_out: f_out.write( self.dislocation_dipole_template.format( input_directory=self.lammps_input_directory, lattice_structure=lattice_full_name, lattice_constant=self.lattice_constant, poisson_ratio=nu, x_direction_0=x_direction[0], x_direction_1=x_direction[1], x_direction_2=x_direction[2], dipole_direction_0=dipole_direction[0], dipole_direction_1=dipole_direction[1], dipole_direction_2=dipole_direction[2], dislocation_line_direction_0=dislocation_line_direction[0], dislocation_line_direction_1=dislocation_line_direction[1], dislocation_line_direction_2=dislocation_line_direction[2], length_b01=n1, length_b02=n2, length_b03=n3, bs0=bs[0], bs1=bs[1], bs2=bs[2], mult=round(mult, 3), ) ) try: # Run the MD++ with the created input file subprocess.call(["md_gpp", dipole_creation_file]) # subprocess.check_call(["md_mac", dipole_creation_file]) except subprocess.CalledProcessError: raise DislocationError("Failed to run MD++.") def cg_min_with_lj_potential(self, mult): """Apply a cg minimization algorithm. Apply a cg minimization algorithm using an ad-hoc lj potential to relax the system after dislocation dipole creation from the MD++. Args: mult (float): multiplication factor Returns: str: xyz_file_relaxed name of the xyz_file after relaxation """ dipole_lammps_ref_file = join( self.lammps_input_directory, "dipole_ref", "dipole_{}.lammps".format(round(mult, 3))) if not os.path.isfile("{}.gz".format(dipole_lammps_ref_file)): emsg = "The file:\n{}.gz".format(dipole_lammps_ref_file) emsg += "\ndoes not exist." raise DislocationError(emsg) with gzip.open("{}.gz".format(dipole_lammps_ref_file), "rb",) as f_in: with open(dipole_lammps_ref_file, "wb",) as f_out: # Copy the contents f_in to f_out. shutil.copyfileobj(f_in, f_out) xyz_file = join(self.lammps_input_directory, "dipole_b01_{}.lammps".format(round(mult, 3))) # 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(xyz_file, "w") as f_out: d_file = open(dipole_lammps_ref_file).read() f_out.write(d_file.format()) xyz_file_relaxed = \ join(self.lammps_input_directory, "dipole_b01_{}_relaxed.lammps".format(round(mult, 3))) in_lammps = join(self.lammps_input_directory, "in.min_cg_1") # Create the LAMMPS input file with open(in_lammps, "w") as f_out: try: epsilon, sigma, cutoff = \ get_lj_parameters(self.lattice_type, self.lattice_constant, self.cohesive_energy) except: raise DislocationError("Failed to get the lj parameters.") f_out.write( self.min_cg_1_template.format( xyzfile=xyz_file, mass=self.mass, cutoff=cutoff, epsilon=epsilon, sigma=sigma, writefile=xyz_file_relaxed, ) ) log_lammps = \ join(self.lammps_output_directory, "dipole_b01_{}_relaxed.lammps.log".format(round(mult, 3))) # Running the lammps script num_procs = get_number_of_processors_to_use(mult, self.max_procs) try: run_lammps(in_lammps, log_lammps, lammps_exe="lmp", num_procs=num_procs) except: raise DislocationError("Failed to run lammps.") return xyz_file_relaxed def extract_result(self, log_lammps, n1): """Process the lammps log file and extract the relevant information. Process the output dumpfile for relevant information. Get the Dislocation energy. Args: log_lammps (str): LAMMPS log file n1 (float): supercell size in the x direction Returns: list: useful_value """ main_result = None try: lines = open(log_lammps).readlines() for line in lines: line_split = line.split() if len(line_split) <= 10: continue zzz = [] try: for item in line_split: zzz.append(float(item)) except: continue if main_result is None: main_result = zzz elif zzz[0] >= main_result[0]: main_result = zzz except: emsg = "Failed to find the dislocation " emsg += "energy in the LAMMPS output log file:\n" emsg += log_lammps raise DislocationError(emsg) if main_result is None: emsg = "Failed to find the dislocation " emsg += "energy in the LAMMPS output log file" raise DislocationError(emsg) if verbose: print("main_result=\n{}".format(main_result)) useful_value = [ n1, # supercell size along the x-direction main_result[1], # number of atoms -self.cohesive_energy, 0.0, main_result[4], # PotEng main_result[8], # Pxx main_result[9], # Pyy main_result[10], # Pzz main_result[11], # Pxy main_result[12], # Pyz main_result[13], # Pxz ] return useful_value def compute_info_and_energy(self, x_length, y_length, z_length, burgers_vector_in_scaled_coordinates, x_direction, dipole_direction, dislocation_line_direction, mult, nu): """Compute the total energy of supercell containing dipole structure and the information. Args: x_length (float): supercell size in the x direction y_length (float): supercell size in the y direction z_length (float): supercell size in the z direction burgers_vector_in_scaled_coordinates (1darray): Burgers vector in scaled coordinates dipole_direction (1darray): dipole direction x_direction (1darray): x direction dislocation_line_direction (1darray): dislocation line direction angle (float) : supercell angle mult (float): multiplication factor nu (float) : poisson ratio Returns: 1darray, list: coord, useful_value coord (1darray): [ x0, x1, y1, x2r, y2r ] useful_value (list) """ # Create dipole using MD++ try: self.create_dipole(x_length, y_length, z_length, burgers_vector_in_scaled_coordinates, x_direction, dipole_direction, dislocation_line_direction, mult, nu) except: raise DislocationError("Failed to create dipole using MD++.") if verbose: print("Finished creating dipole using MD++") # Apply a cg minimization to the MD++ output using an ad-hoc # lj potential. try: xyz_file = self.cg_min_with_lj_potential(mult) except: raise DislocationError("Failed to minimize the system.") if verbose: print("Finished lj relaxation using LAMMPS") dump_file = \ join(self.dumping_directory, "dislocation_dipole_npt_{}.dump.gz".format(round(mult, 3))) in_lammps = join(self.lammps_input_directory, "in.min_cg_npt") # Create the LAMMPS input file with open(in_lammps, "w") as f_out: f_out.write( self.min_cg_npt_template.format( xyzfile=xyz_file, modelname=self.kim_model, symbol=self.species, lattice_type=self.lattice_type, dumpfile=dump_file, ) ) log_lammps = \ join(self.lammps_output_directory, "equilibrated_dislocation_dipole_npt_" "{}.log".format(round(mult, 3))) # Running the lammps script num_procs = get_number_of_processors_to_use(mult, self.max_procs) try: run_lammps(in_lammps, log_lammps, lammps_exe="lmp", num_procs=num_procs) except: raise DislocationError("Failed to run lammps.") if verbose: print("Finished running LAMMPS") try: method = 'cna' if self.lattice_type in ('fcc', 'bcc') else 'dxa' pos, cvec = extract_position(dump_file, method=method, lattice=self.lattice_type) except: raise DislocationError("Failed to extract position.") # Supercell orientation # Dipole position x0 = cvec[0][0] x1 = cvec[1][0] y1 = cvec[1][1] r2 = pos[1] - pos[0] if pos[0][0] < pos[1][0] else pos[0] - pos[1] S = np.array([[x0, x1], [0., y1]]) S_inv = np.linalg.inv(S) dr = np.matmul(S_inv, r2) if verbose: print("S = {}, r2 = {}, dr = {}".format(S, r2, dr)) del S, S_inv, r2 x2r = dr[0] y2r = dr[1] coord = [x0, x1, y1, x2r, y2r] try: useful_value = self.extract_result(log_lammps, x_length) except: emsg = "Failed to extract results from the lammps log file." raise DislocationError(emsg) return coord, useful_value def num_nearest_neighbors_in_reference_structure(lattice): """Get the number of nearest neighbors in the reference structure. Args: lattice (str): lattice type Returns: int: number of nearest neighbors in the reference structure """ if lattice not in ('fcc', 'bcc', 'diamond'): emsg = 'Lattice {} is not supported.'.format(lattice) emsg += 'Valid lattices are:\n\t- ' emsg += '{}'.format('\n\t- '.join(('fcc', 'bcc', 'diamond'))) raise DislocationError(emsg) num_nearest_neighbors = { 'fcc': 12, 'bcc': 8, 'diamond': 4, } return num_nearest_neighbors[lattice] def crystal_structure_radius(lattice, lattice_constant, shell_number): """Get the radius for shells of the `fcc, bcc, diamond` crystal strucrues. Args: lattice (str): lattice type lattice_constant (float): equilibrium lattice constant shell_number (int): shell number Returns: float: radius """ if lattice not in ('fcc', 'bcc', 'diamond'): emsg = 'Lattice {} is not supported.'.format(lattice) emsg += 'Valid lattices are:\n\t- ' emsg += '{}'.format('\n\t- '.join(('fcc', 'bcc', 'diamond'))) raise DislocationError(emsg) if lattice_constant <= 0: emsg = "Invalid lattice_constant = {}.".format(lattice_constant) raise DislocationError(emsg) if shell_number not in (1, 2, 3, 4): raise DislocationError("Invalid shell number.") CST = { 'fcc': { 1: lattice_constant / np.sqrt(2), 2: lattice_constant, 3: lattice_constant * np.sqrt(3 / 2), 4: lattice_constant * np.sqrt(2)}, 'bcc': { 1: lattice_constant * np.sqrt(3 / 2), 2: lattice_constant, 3: lattice_constant * np.sqrt(2), 4: lattice_constant * np.sqrt(11 / 4)}, 'diamond': { 1: lattice_constant * np.sqrt(3 / 4), 2: lattice_constant / np.sqrt(2), 3: lattice_constant * np.sqrt(11 / 4), 4: lattice_constant}} return CST[lattice][shell_number] def get_lj_parameters(lattice, lattice_constant, cohesive_energy): """Get the parameters of the standard 12/6 Lennard-Jones potential. Args: lattice (str): lattice type lattice_constant (float): equilibrium lattice constant cohesive_energy (float): cohesive energy Returns: float, float, float: spsilon, sigma, cutoff spsilon: energy units. sigma: distance units. `sigma` is defined in the LJ formula as the zero-crossing distance for the potential cutoff: distance units. """ # Approximate the epsilon nn1 = num_nearest_neighbors_in_reference_structure(lattice) epsilon = cohesive_energy * 2.0 / nn1 # Approximate the sigma csr_1 = crystal_structure_radius(lattice, lattice_constant, 1) sigma = csr_1 / 2**(1/6) # Approximate the cutoff csr_2 = crystal_structure_radius(lattice, lattice_constant, 2) cutoff = 0.5 * (csr_1 + csr_2) return epsilon, sigma, cutoff def create_multiplication_factor_list(start, increment, num_multiplication_factor, increment_max=0.5, multiplication_limit=2.0): """Create a list of multiplication factors. Args: start (float): starting value in the multiplication factors array. increment (float): incremental value num_multiplication_factor (int): number of multiplication factors increment_max (float, optional): maximum incremental value. (default: {0.5}). multiplication_limit (float, optional): multiplication limit, where from that point we use the maximum incremental value. (default: {2.0}). Returns: 1darray, float: multiplication factors, incremental value """ mfv = [start, *[increment] * num_multiplication_factor] mults = np.add.accumulate(mfv) for m in range(mults.size): mults[m] = round(mults[m], 3) if increment < increment_max: increment = increment_max mults_indices = np.where(mults > multiplication_limit)[0] if len(mults_indices) > 0: mfv = np.array(mfv) mfv[mults_indices] = increment mults = np.add.accumulate(mfv) for m in range(mults.size): mults[m] = round(mults[m], 3) return mults, increment # initial global variables multiplication_factor_start = 1.0 multiplication_factor_increment = 0.1 supercell_size = [np.zeros(3)] # Indicating to use the convergence algorithm or not # # Note: # If it is set to False you only need to set the initial run length to be # used as the number of multiplication factor starting from # multiplication_factor_start and increamenting with # multiplication_factor_increment # # Note: # There is a default limit on multiplication_factor_increment which means by # default multiplication_factor_increment is set to 0.5 after reaching the # multiplication_limit = 2. # These are set in create_multiplication_factor_list routine. If different # behavior is expected please update these values! use_convergence = True # Initial run length initial_run_length = 15 # The maximum run length represents a cost constraint. maximum_run_length = 21 # The maximum number of steps as an equilibration hard limit. If the # algorithm finds equilibration_step greater than this limit it will fail. # For the default None, the function is using `maximum_run_length // 2` as # the maximum equilibration step. maximum_equilibration_step = 20 # A relative half-width requirement or the accuracy parameter. Target value # for the ratio of halfwidth to sample mean. If n_variables > 1, # relative_accuracy can be a scalar to be used for all variables or a 1darray # of values of size n_variables. relative_accuracy = 0.1 if __name__ == "__main__": # Read elements, masses, and model from stdin element = input("Elements: ") if verbose: print("element = {}".format(element)) # Remove the space from input. species = " ".join(element.split()) mass = input("Masses in g/mol (must match the order of elements above): ") if verbose: print("mass = {}".format(mass)) model = input("Model: ") if verbose: print("model = {}".format(model)) lattice_type = input("Lattice type: ") if verbose: print("lattice_type = {}".format(lattice_type)) # Get the equilibrium lattice constant and cohesive energy equilibrium_lattice_constant, zero_temperature_cohesive_energy = [ float(val) for val in json.loads(input("Lattice constant (m) and cohesive energy (J): "))] # Convert lattice constant from m to angstroms equilibrium_lattice_constant *= 1e10 # Convert cohesive energy from J to eV zero_temperature_cohesive_energy *= 6.2415091e18 if verbose: print("equilibrium_lattice_constant = {} angstroms".format( equilibrium_lattice_constant)) print("cohesive_energy = {} eV".format( zero_temperature_cohesive_energy)) # Get the isothermal cubic elastic constants for this lattice isothermal_elastic_constants = [ float(val) for val in json.loads(input("Cubic isothermal elastic constants [C11, C12, C44] (Pa): "))] # Convert elastic constants from Pa to GPa isothermal_elastic_constants = [ val * 1e-9 for val in isothermal_elastic_constants] if verbose: print("isothermal elastic constants = {} GPa ".format( isothermal_elastic_constants)) dislocation_core_radius_str = input("dislocation core radius: ") if verbose: print("dislocation_core_radius_str = {}".format( dislocation_core_radius_str)) dislocation_core_radius = parse_coeffs(dislocation_core_radius_str) dislocation_core_radius_size = len(dislocation_core_radius) # The Burgers vector of the dislocation given as a vector of three real # numbers relative to the lattice parameter, e.g. [0.5, 0.5, 0] corresponds # to a Burgers vectors of [a/2, a/2, 0] burgers_vector_direction_str = input("burgers vector direction: ") if verbose: print("burgers_vector_direction_str = {}".format( burgers_vector_direction_str)) burgers_vector_vector = parse_coeffs(burgers_vector_direction_str) slip_plane_miller_indices_str = input("slip plane miller indices: ") if verbose: print("slip_plane_miller_indices_str = {}".format( slip_plane_miller_indices_str)) slip_plane_direction_vector = parse_coeffs(slip_plane_miller_indices_str) dislocation_line_direction_str = input("dislocation line direction: ") if verbose: print("dislocation_line_direction_str = {}".format( dislocation_line_direction_str)) dislocation_line_direction_vector = parse_coeffs( dislocation_line_direction_str) max_num_procs = max(1, int(input("max num processors: "))) if verbose: print("max_num_procs = {}".format(max_num_procs)) ###################### # Postprocess inputs ###################### # Check consistency of burgers vector, dislocation line direction, and slip # plane miller indices check_consistency_of_dislocation_parameters( burgers_vector_vector, dislocation_line_direction_vector, slip_plane_direction_vector) shear_modulus, poisson_ratio = get_shear_modulus_and_poisson_ratio( isothermal_elastic_constants) # cubic elastic constants in units required to MadSum elastic_constants_in_madsum_units = np.array( isothermal_elastic_constants) / 160.22 # Shear modulus in units required to MadSum shear_modulus_in_madsum_units = shear_modulus / 160.22 # dipole_c1: x direction # dipole_c2: dipole direction # dipole_c3: line direction dipole_c1, dipole_c2, dipole_c3, supercell_angle = \ generate_dipole_parameters(burgers_vector_vector, dislocation_line_direction_vector, slip_plane_direction_vector) dipole_c1_norm = np.linalg.norm(dipole_c1) dipole_c2_norm = np.linalg.norm(dipole_c2) dipole_c3_norm = np.linalg.norm(dipole_c3) alatt = equilibrium_lattice_constant * dipole_c1_norm # convert alat to m alatt *= 1e-10 # <110> direction is the dipole direction -> b01 blatt = equilibrium_lattice_constant * dipole_c2_norm # convert alat to m blatt *= 1e-10 # <112> direction -> b02 clatt = equilibrium_lattice_constant * dipole_c3_norm # convert alat to m clatt *= 1e-10 # vol in m^3 abc_vol = alatt * blatt * clatt # geometry of the simulation box # burgers vector burgers_vector = 0.5 * np.array(burgers_vector_vector) burgers_vector_norm = np.linalg.norm(burgers_vector) burgers_vector_in_unit_cell = equilibrium_lattice_constant * burgers_vector burgers_vector_in_unit_cell_magnitude = np.linalg.norm( burgers_vector_in_unit_cell) # M: Coordinate transformation matrix coordinate_transformation_matrix = \ np.vstack((dipole_c1 / dipole_c1_norm, dipole_c2 / dipole_c2_norm, dipole_c3 / dipole_c3_norm)).T b_madsum = np.matmul(coordinate_transformation_matrix.T, burgers_vector_in_unit_cell) # Obtaining angle between axes and b = [1 1 1] / 2 delta = math.acos( round(np.dot(burgers_vector, dipole_c1) / (burgers_vector_norm * dipole_c1_norm), 8)) if verbose: print("delta = angle(b, x) [degrees] = {} ".format(delta * 180 / PI)) phi = math.acos( round(np.dot(burgers_vector, dipole_c2) / (burgers_vector_norm * dipole_c2_norm), 8)) if verbose: print("phi = angle(b, y) [degrees] = {} ".format(phi * 180 / PI)) theta = math.acos( round(np.dot(burgers_vector, dipole_c3) / (burgers_vector_norm * dipole_c3_norm), 8)) if verbose: print("theta = angle(b, z) [degrees] = {} ".format(theta * 180 / PI)) # burgers vector in MADSUM format burgers_vector_madsum = burgers_vector_norm * \ equilibrium_lattice_constant * \ np.array([math.cos(delta), math.cos(phi), math.cos(theta)]) # C in GPa c11 = isothermal_elastic_constants[0] * 1000 c12 = isothermal_elastic_constants[1] * 1000 c44 = isothermal_elastic_constants[2] * 1000 # C: Voigt elasticity matrix voigt_elasticity_matrix = np.array([[c11, c12, c12, 0.0, 0.0, 0.0], [c12, c11, c12, 0.0, 0.0, 0.0], [c12, c12, c11, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, c44, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, c44, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, c44]]) # The Voigt elasticity matrix inverse voigt_elasticity_matrix_inv = np.linalg.inv(voigt_elasticity_matrix) del c11, c12, c44, voigt_elasticity_matrix E_core_final_nonsingular = [] E_core_final_iso = [] E_core_final_cub = [] super_cell = SuperCell(species, mass, model, lattice_type, equilibrium_lattice_constant, zero_temperature_cohesive_energy, max_num_procs) def get_trajectory(step): # multiplication_factor global multiplication_factor_start global multiplication_factor_increment global supercell_size supercell_size = [supercell_size[-1]] multiplication_factor_list = [] multiplication_factor_list_size = 0 final_value = [] coords = [] length_b03 = [] # Create a list of multiplication factors mults, multiplication_factor_increment = \ create_multiplication_factor_list(multiplication_factor_start, multiplication_factor_increment, step) # Update the start factor for the next call multiplication_factor_start = mults[-1] aspect_ratio = 1.5 # aspect ratio x_length / y_length x_length = 10.0 # supercell size along the x-direction z_length = 3.0 # supercell size along the z-direction while True: for mult_ind, mult in enumerate(mults[:-1]): n1, n2, n3, bs = get_supercell_info(dipole_c1, dipole_c2, dipole_c3, supercell_angle, mult, aspect_ratio, x_length, z_length) # The same cell size if np.array_equal([n1, n2, n3], supercell_size[-1]): if verbose: wmsg = "The created super cell has the same size of " wmsg += "({}, {}, {}), as the ".format(n1, n2, n3) wmsg += "previous cell." print_warning(wmsg) continue try: coord, useful_value = \ super_cell.compute_info_and_energy( n1, n2, n3, bs, dipole_c1, dipole_c2, dipole_c3, mult, poisson_ratio) except: if verbose: print_warning("Failed to compute info and energy.") continue multiplication_factor_list.append(mult) supercell_size.append(np.array([n1, n2, n3])) length_b03.append(n3) coords.append(coord) final_value.append(useful_value) multiplication_factor_list_size = \ len(multiplication_factor_list) # If we have enough samples, we need to update the # starting index and the corresponding increment if multiplication_factor_list_size == step: multiplication_factor_start = mults[mult_ind + 1] break # We do not have enough samples to stop if multiplication_factor_list_size < step: # Create a list of multiplication factors mults, multiplication_factor_increment = \ create_multiplication_factor_list( multiplication_factor_start, multiplication_factor_increment, step) # Update the start factor for the next call multiplication_factor_start = mults[-1] continue break trajectory = np.empty((3 * dislocation_core_radius_size, multiplication_factor_list_size), dtype=np.float64) for core_radius_ind, core_radius in enumerate(dislocation_core_radius): # cut-off radius cutoff_radius = core_radius * burgers_vector_in_unit_cell_magnitude # truncation paramters truncation = (10, 10) E_min_dislocation_npt_cea = np.array(final_value, copy=True) Eel_nonsingular = np.empty(multiplication_factor_list_size) Eprm_nonsingular = np.empty(multiplication_factor_list_size) Eimg_nonsingular = np.empty(multiplication_factor_list_size) Eel_iso = np.empty(multiplication_factor_list_size) Eprm_iso = np.empty(multiplication_factor_list_size) Eimg_iso = np.empty(multiplication_factor_list_size) Eel_cub = np.empty(multiplication_factor_list_size) Eprm_cub = np.empty(multiplication_factor_list_size) Eimg_cub = np.empty(multiplication_factor_list_size) E_residual = np.empty(multiplication_factor_list_size) E_atom_A = np.empty(multiplication_factor_list_size) if verbose: print("-" * 100) print('core radius = {}'.format(core_radius)) for mult in range(multiplication_factor_list_size): coord = coords[mult] total_elastic_energy_nonsingular, \ primary_dipole_energy_nonsingular, \ image_energy_nonsingular = madsum_nonsingular( shear_modulus_in_madsum_units, poisson_ratio, b_madsum, cutoff_radius, coord, truncation) Eel_nonsingular[mult] = total_elastic_energy_nonsingular Eprm_nonsingular[mult] = primary_dipole_energy_nonsingular Eimg_nonsingular[mult] = image_energy_nonsingular if verbose: print('multiplication factor = {}'.format(mult)) print("Eel_nonsingular = {}".format( Eel_nonsingular[:mult + 1])) total_elastic_energy_iso, \ primary_dipole_energy_iso, \ image_energy_iso = madsum_iso( shear_modulus_in_madsum_units, poisson_ratio, b_madsum, cutoff_radius, coord, truncation) Eel_iso[mult] = total_elastic_energy_iso Eprm_iso[mult] = primary_dipole_energy_iso Eimg_iso[mult] = image_energy_iso if verbose: print("Eel_iso = {}".format(Eel_iso[:mult + 1])) total_elastic_energy_cubic, \ primary_dipole_energy_cubic, \ image_energy_cubic = madsum_cubic( elastic_constants_in_madsum_units, burgers_vector_in_unit_cell, coordinate_transformation_matrix, cutoff_radius, coord, truncation) Eel_cub[mult] = total_elastic_energy_cubic Eprm_cub[mult] = primary_dipole_energy_cubic Eimg_cub[mult] = image_energy_cubic if verbose: print("Eel_cub = {}".format(Eel_cub[:mult + 1])) print("CC={}".format(elastic_constants_in_madsum_units)) print("b0={}".format(burgers_vector_in_unit_cell)) print("M={}".format(coordinate_transformation_matrix)) print("rc={}".format(cutoff_radius)) print("b_madsum={}".format(b_madsum)) print("shear modulus={}".format( shear_modulus_in_madsum_units)) print("poisson ratio={}".format(poisson_ratio)) print("coord={}".format(coord)) print("truncation={}".format(truncation)) print("E_min_dislocation_npt_cea = {}".format( E_min_dislocation_npt_cea)) # Pij in LAMMPS (metal units) are in bar=0.1 MPa sigma = \ [E_min_dislocation_npt_cea[mult, 5 + i] / 10 for i in range(6)] # vol in m^3 vol = np.prod(supercell_size[mult]) * abc_vol # sigma * voigt_elasticity_matrix_inv * sigma' (original) # E_residual in MPa = MJ/m^3 E_residual_temp = np.dot( np.dot(sigma, voigt_elasticity_matrix_inv), np.transpose(sigma)) E_residual_temp = E_residual_temp * \ (10 ** 6) * vol / (1.602e-19) if verbose: print("E_residual_temp = {}".format(E_residual_temp)) # E_residual in J E_residual[mult] = E_residual_temp # E_residual in eV E_atm_A_temp = ( np.array(E_min_dislocation_npt_cea[mult, 4]) - np.array(np.multiply(E_min_dislocation_npt_cea[mult, 1], E_min_dislocation_npt_cea[mult, 2])) - np.transpose(E_residual_temp) / 2 ) / (length_b03[mult] * clatt*1e10) if verbose: print("E_atm_A_temp = {}".format(E_atm_A_temp)) E_atom_A[mult] = E_atm_A_temp if verbose: print("-" * 50) print("E_atom_A = {}".format(E_atom_A)) print("E_residual = {}".format(E_residual)) E_elastic_nonsingular = np.transpose(Eel_nonsingular) E_elastic_iso = np.transpose(Eel_iso) E_elastic_cub = np.transpose(Eel_cub) E_core_A_nonsingular = (E_atom_A - E_elastic_nonsingular) / 2 if verbose: print("E_core_A_nonsingular = {}".format(E_core_A_nonsingular)) trajectory[core_radius_ind * 3, :] = E_core_A_nonsingular[:] E_core_A_iso = (E_atom_A - E_elastic_iso) / 2 if verbose: print("E_core_A_iso = {}".format(E_core_A_iso)) trajectory[core_radius_ind * 3 + 1, :] = E_core_A_iso[:] E_core_A_cub = (E_atom_A - E_elastic_cub) / 2 if verbose: print("E_core_A_cub = {}".format(E_core_A_cub)) trajectory[core_radius_ind * 3 + 2, :] = E_core_A_cub[:] return trajectory if use_convergence: try: msg = cr.run_length_control( get_trajectory=get_trajectory, n_variables=3 * dislocation_core_radius_size, initial_run_length=initial_run_length, run_length_factor=1, maximum_run_length=maximum_run_length, maximum_equilibration_step=maximum_equilibration_step, relative_accuracy=relative_accuracy, confidence_interval_approximation_method='subsample', batch_size=2, fp='return', fp_format='edn') except Exception as e: msg = '{}'.format(e) raise cr.CVGError(msg) else: _traj = get_trajectory(step=initial_run_length) if not np.all(np.isfinite(_traj)): msg = 'there is/are value/s in the input which is/are ' msg += 'non-finite or not number.' raise DislocationError(msg) traj = np.array(_traj, dtype=np.float64) msg = {"converged": False} for i in range(3 * dislocation_core_radius_size): msg[i] = { "mean": np.mean(traj[i, :]), "standard_deviation": np.std(traj[i, :]), } msg = kim_edn.dumps(msg) kim_obj = kim_edn.loads(msg) E_core_energy_nonsingular = [] E_core_energy_nonsingular_std = [] E_core_energy_isotropic = [] E_core_energy_isotropic_std = [] E_core_energy_anisotropic = [] E_core_energy_anisotropic_std = [] for ind in range(dislocation_core_radius_size): key = '{}'.format(ind * 3) E_core_energy_nonsingular.append(kim_obj[key]["mean"]) E_core_energy_nonsingular_std.append( kim_obj[key]["standard_deviation"]) key = '{}'.format(ind * 3 + 1) E_core_energy_isotropic.append(kim_obj[key]["mean"]) E_core_energy_isotropic_std.append( kim_obj[key]["standard_deviation"]) key = '{}'.format(ind * 3 + 2) E_core_energy_anisotropic.append(kim_obj[key]["mean"]) E_core_energy_anisotropic_std.append( kim_obj[key]["standard_deviation"]) kim_property_name = "dislocation-core-energy-cubic-crystal-npt" dump_property(kim_property_name, lattice_type, species, equilibrium_lattice_constant, slip_plane_direction_vector, dislocation_line_direction_vector, burgers_vector_vector, dislocation_core_radius, E_core_energy_nonsingular, E_core_energy_nonsingular_std, E_core_energy_isotropic, E_core_energy_isotropic_std, E_core_energy_anisotropic, E_core_energy_anisotropic_std)