#!/usr/bin/env python """ Vacancy Formation Energy (VFE) and Relaxation Volume (VRV) Test Driver Works for both cubic and hcp Crystals, at 0 K and 0 GPa This version speeds up the calculation by using more potential energy evaluations to decrease the steps required in FIRE relaxation Date: 2015/09/02 Author: Junhao Li , Cornell University Updated 2018/04/02 Daniel S. Karls, University of Minnesota Last Update: 2023/06/23 Eric Fuemmeler, University of Minnesota """ # ASE Modules from ase.build import bulk from ase.optimize import FIRE # KIM Modules from ase.calculators.kim.kim import KIM # Python Modules import os from scipy.optimize import fmin import numpy as np import sys import re import json import math from collections import OrderedDict # Parameters for Production FIRE_LOG = 'fire.log' FIRE_MAX_STEPS = 50 FIRE_UNCERT_STEPS = 20 FIRE_TOL = 1e-3 # absolute FMIN_FTOL = 1e-6 # relative FMIN_XTOL = 1e-10 # relative VFE_TOL = 1e-5 # absolute MAX_LOOPS = 20 CELL_SIZE_MIN = 3 CELL_SIZE_MAX = 5 COLLAPSE_CRITERIA_VOLUME = 0.1 COLLAPSE_CRITERIA_ENERGY = 0.1 DYNAMIC_CELL_SIZE = True # Increase Cell Size According to lattice structure EPS = 1e-3 # Parameters for Debugging # FIRE_MAX_STEPS = 200 # FIRE_TOL = 1e-3 # absolute # FMIN_FTOL = 1e-3 # relative # FMIN_XTOL = 1e-5 # relative # Extrapolation Parameters FITS_CNT = [2, 3, 3, 3, 3] # Number of data points used for each fitting FITS_ORDERS = [ [0, 3], [0, 3], [0, 3, 4], [0, 3, 5], [0, 3, 6], ] # Number of orders included in each fitting # Fit Results Used (Corresponding to the above) FITS_VFE_VALUE = 0 # Vacancy Formation Energy FITS_VFE_UNCERT = [1, 2] FITS_VRV_VALUE = 0 # Vacancy Relaxation Volume FITS_VRV_UNCERT = [1, 2] # Strings for Output KEY_SOURCE_VALUE = 'source-value' KEY_SOURCE_UNIT = 'source-unit' KEY_SOURCE_UNCERT = 'source-std-uncert-value' UNIT_ENERGY = 'eV' UNIT_LENGTH = 'angstrom' UNIT_ANGLE = 'degree' UNIT_PRESSURE = 'GPa' UNIT_VOLUME = UNIT_LENGTH + '^3' SPACE_GROUPS = { 'fcc': 'Fm-3m', 'bcc': 'Im-3m', 'sc': 'Pm-3m', 'diamond': 'Fd-3m', 'hcp': 'P63/mmc', } WYCKOFF_CODES = { 'fcc': ['4a'], 'bcc': ['2a'], 'sc': ['1a'], 'diamond': ['8a'], 'hcp': ['2d'], } WYCKOFF_SITES = { 'fcc': [[0.0, 0.0, 0.0]], 'bcc': [[0.0, 0.0, 0.0]], 'sc': [[0.0, 0.0, 0.0]], 'diamond': [[0.0, 0.0, 0.0]], 'hcp': [[2.0 / 3.0, 1.0 / 3.0, 0.25]], } VUFE_PROP_ID = 'tag:staff@noreply.openkim.org,2015-07-28:property/monovacancy-neutral-unrelaxed-formation-potential-energy-crystal-npt' VFE_PROP_ID = 'tag:staff@noreply.openkim.org,2015-07-28:property/monovacancy-neutral-relaxed-formation-potential-energy-crystal-npt' VRV_PROP_ID = 'tag:staff@noreply.openkim.org,2015-07-28:property/monovacancy-neutral-relaxation-volume-crystal-npt' class VacancyCalculation(object): # Class for calculating vacancy formation energy and relaxation volume def __init__(self, calc, elem, model, lattice, latticeConsts): self.calc = calc self.elem = elem self.model = model self.lattice = lattice self.latticeConsts = latticeConsts self.VFEUncert = 0 self.VRVUncert = 0 if lattice == 'hcp': atoms = bulk( elem, a = latticeConsts[0], c = latticeConsts[1], crystalstructure = 'hcp', ) else: atoms = bulk( elem, a = latticeConsts[0], crystalstructure = lattice, cubic = True, ) atoms.set_calculator(calc) self.atoms = atoms if DYNAMIC_CELL_SIZE == True: numAtoms = atoms.get_number_of_atoms() factor = math.pow(8 / numAtoms, 0.333) global CELL_SIZE_MIN, CELL_SIZE_MAX CELL_SIZE_MIN = int(math.ceil(factor * CELL_SIZE_MIN)) CELL_SIZE_MAX = CELL_SIZE_MIN + 2 print('CELL_SIZE_MIN:', CELL_SIZE_MIN) print('CELL_SIZE_MAX:', CELL_SIZE_MAX) print('Smallest System Size:', numAtoms * CELL_SIZE_MIN**3) print('Largest System Size:', numAtoms * CELL_SIZE_MAX**3) def _createSupercell(self, size): atoms = self.atoms.copy() atoms.set_calculator(self.calc) atoms *= (size, size, size) return atoms def _cellVector2Cell(self, cellVector): # Reconstruct cell From cellVector cell = [ [cellVector[0], 0, 0], [cellVector[1], cellVector[2], 0], [cellVector[3], cellVector[4], cellVector[5]] ] return cell def _cell2CellVector(self, cell): # Extract cellVector From cell # For reducing degree of freedom during relaxation cellVector = [ cell[0, 0], cell[1, 0], cell[1, 1], cell[2, 0], cell[2, 1], cell[2, 2], ] return cellVector def _getVFE(self, cellVector, atoms, enAtoms, numAtoms): newCell = self._cellVector2Cell(cellVector) atoms.set_cell(newCell, scale_atoms = True) enAtomsWithVacancy = atoms.get_potential_energy() enVacancy = enAtomsWithVacancy - enAtoms * (numAtoms - 1) / numAtoms return enVacancy def _getResultsForSize(self, size): # Setup Environment unrelaxedCell = self.atoms.get_cell() * size unrelaxedCellVector = self._cell2CellVector(unrelaxedCell) atoms = self._createSupercell(size) numAtoms = atoms.get_number_of_atoms() enAtoms = atoms.get_potential_energy() unrelaxedCellEnergy = enAtoms unrelaxedCellVolume = np.abs(np.linalg.det(unrelaxedCell)) print('\nSupercell Size:\n', size) print('Unrelaxed Cell:\n', unrelaxedCell) print('Unrelaxed Cell Vector:\n', unrelaxedCellVector) print('Unrelaxed Cell Energy:\n', unrelaxedCellEnergy) # Create Vacancy del atoms[0] enAtomsWithVacancy = atoms.get_potential_energy() print('Energy of Unrelaxed Cell With Vacancy:\n', unrelaxedCellEnergy) enVacancyUnrelaxed = enAtomsWithVacancy - enAtoms * (numAtoms - 1) / numAtoms # Self Consistent Relaxation enVacancy = 0 relaxedCellVector = unrelaxedCellVector loop = 0 while 1: # Position Relaxation print('==========') print('Loop:', loop) print('Position Relaxation...') dyn = FIRE(atoms) # dyn = FIRE(atoms, logfile = FIRE_LOG) dyn.run(fmax = FIRE_TOL, steps = FIRE_MAX_STEPS) numSteps = dyn.get_number_of_steps() if numSteps >= FIRE_MAX_STEPS: print('WARNING: Max number of steps exceeded. Structure may be unstable.') # sys.exit(0) print('Relaxation Completed. Steps:', numSteps) # Cell Size Relaxation print('Cell Size Relaxation...') tmpCellVector, tmpEnVacancy = fmin( self._getVFE, relaxedCellVector, args = (atoms, enAtoms, numAtoms), ftol = FMIN_FTOL, xtol = FMIN_XTOL, full_output = True, )[:2] # Convergence Requirement Satisfied if abs(tmpEnVacancy - enVacancy) < VFE_TOL and dyn.get_number_of_steps() < 1: dyn.run(fmax = FIRE_TOL * EPS, steps = FIRE_UNCERT_STEPS) tmpCellVector, tmpEnVacancy = fmin( self._getVFE, relaxedCellVector, args = (atoms, enAtoms, numAtoms), ftol = FMIN_FTOL * EPS, xtol = FMIN_XTOL * EPS, full_output = True, )[:2] self.VFEUncert = np.abs(tmpEnVacancy - enVacancy) enVacancy = tmpEnVacancy oldVolume = np.linalg.det(self._cellVector2Cell(relaxedCellVector)) newVolume = np.linalg.det(self._cellVector2Cell(tmpCellVector.tolist())) self.VRVUncert = np.abs(newVolume - oldVolume) relaxedCellVector = tmpCellVector.tolist() break enVacancy = tmpEnVacancy relaxedCellVector = tmpCellVector.tolist() # Check Loop Limit loop += 1 if loop > MAX_LOOPS: print('Loops Limit Exceeded. Structure Unstable.') sys.exit(0) # Output Temporary Result relaxedCell = self._cellVector2Cell(relaxedCellVector) relaxedCellVolume = np.abs(np.linalg.det(relaxedCell)) relaxationVolume = unrelaxedCellVolume - relaxedCellVolume print('Current VFE:', enVacancy) print('Energy of Supercell:', enAtoms) print('Unrelaxed Cell Volume:', unrelaxedCellVolume) print('Current Relaxed Cell Volume:', relaxedCellVolume) print('Current Relaxation Volume:', relaxationVolume) print('Current Cell:\n', np.array(self._cellVector2Cell(relaxedCellVector))) # Determine Collapse if np.abs(relaxationVolume) > COLLAPSE_CRITERIA_VOLUME * unrelaxedCellVolume: print('System Collapsed. Volume significantly changed.') sys.exit(0) if np.abs(enVacancy) > COLLAPSE_CRITERIA_ENERGY * np.abs(enAtoms): print('System Collapsed. System Energy significantly changed.') sys.exit(0) # Print Summary print('---------------') print('Calculation Completed.') print('Number Of Atoms in Supercell:', numAtoms) print('Vacancy Formation Energy (relaxed):', enVacancy) print('Vacancy Formation Energy (unrelaxed):', enVacancyUnrelaxed) print('Unrelaxed Cell Volume:', unrelaxedCellVolume) print('Relaxed Cell Volume:', relaxedCellVolume) print('Relaxation Volume:', relaxationVolume) print('Relaxed Cell Vector:\n', relaxedCellVector) print('Unrelaxed Cell Vector:\n', unrelaxedCellVector) print('Relaxed Cell:\n', np.array(self._cellVector2Cell(relaxedCellVector))) print('Unrelaxed Cell:\n', np.array(self._cellVector2Cell(unrelaxedCellVector))) return enVacancyUnrelaxed, relaxedCellVector, enVacancy, relaxationVolume def _getUnitVector(self, vec): return vec / np.linalg.norm(vec) def _getAngle(self, vec1, vec2): # Get acute angle between two vectors in degrees (always between 0 - 90) vec1Unit = self._getUnitVector(vec1) vec2Unit = self._getUnitVector(vec2) angle = np.arccos(np.dot(vec1Unit, vec2Unit)) if np.isnan(angle): return 0.0 angle = angle * 180.0 / np.pi # if angle < 0: # return 180.0 + angle return angle def _getFit(self, xdata, ydata, orders): # Polynomial Fitting with Specific Orders A = [] print('\nFit with Size:', xdata) print('Orders:', orders) for order in orders: A.append(np.power(xdata * 1.0, -order)) A = np.vstack(A).T print('Matrix A (Ax = y):\n', A) print('Data for Fitting:', ydata) res = np.linalg.lstsq(A, ydata, rcond=None) print('Fitting Results:', res) return res[0] def _getValueUncert(self, valueFitId, uncertFitIds, systematicUncert, maxSizeId, dataSource): # Get sourceValue and sourceUncert use only certain size and fits # Get source value valueFitCnt = FITS_CNT[valueFitId] sourceValue = dataSource[valueFitId][maxSizeId - valueFitCnt + 1] # Get source uncertainty (statistical) sourceUncert = 0 for uncertFitId in uncertFitIds: uncertFitCnt = FITS_CNT[uncertFitId] uncertValue = dataSource[uncertFitId][maxSizeId - uncertFitCnt + 1] sourceUncert = max([abs(uncertValue - sourceValue), sourceUncert]) # Include systematic error, assuming independent of statistical errors sourceUncert = math.sqrt(sourceUncert**2 + systematicUncert**2) return sourceValue, sourceUncert def getResults(self): unitBulk = self.atoms unitCell = unitBulk.get_cell() # Calculate VFE and VRV for Each Size sizes = [] unrelaxedformationEnergyBySize = [] formationEnergyBySize = [] relaxationVolumeBySize = [] print('\n[Calculation]') for size in range(CELL_SIZE_MIN, CELL_SIZE_MAX + 1): unrelaxedFormationEnergy, relaxedCellVector, relaxedFormationEnergy, relaxationVolume = self._getResultsForSize(size) sizes.append(size) unrelaxedformationEnergyBySize.append(unrelaxedFormationEnergy) formationEnergyBySize.append(relaxedFormationEnergy) relaxationVolumeBySize.append(relaxationVolume) print('\n[Calculation Results Summary]') print('Sizes:', sizes) print('Unrelaxed Formation Energy By Size:\n', unrelaxedformationEnergyBySize) print('Formation Energy By Size:\n', formationEnergyBySize) print('Relaxation Volume By Size:\n', relaxationVolumeBySize) # Data for skipping computation when debugging extrapolation and output # sizes = [3, 4, 5, 6, 7, 8, 9] # formationEnergyBySize = [ # 0.6721479768766585 , # 0.67372899358906579, # 0.67440913973746319, # 0.6747228089247983 , # 0.67488432759455463, # 0.6749755557248136 , # 0.67503091578691965, # ] # relaxationVolumeBySize = [ # 8.2664887840680876, # 8.2145358736270282, # 8.2008345712674782, # 8.1943833508903481, # 8.1916426682910242, # 8.1898981954873307, # 8.1889297673697001, # ] # Extrapolate for VFE and VRV of Infinite Size print('\n[Extrapolation]') naSizes = np.array(sizes) naUnrelaxedFormationEnergyBySize = np.array(unrelaxedformationEnergyBySize) naFormationEnergyBySize = np.array(formationEnergyBySize) naRelaxationVolumeBySize = np.array(relaxationVolumeBySize) unrelaxedformationEnergyFitsBySize = [] formationEnergyFitsBySize = [] relaxationVolumeFitsBySize = [] for i in range(0, len(FITS_CNT)): cnt = FITS_CNT[i] # Num of Data Points Used orders = FITS_ORDERS[i] # Orders Included print('Fitting with', cnt, 'points, including orders', orders) unrelaxedformationEnergyFits = [] formationEnergyFits = [] relaxationVolumeFits = [] for j in range(0, len(sizes) - cnt + 1): print('Fit with data beginning', j) xdata = naSizes[j:(j + cnt)] unrelaxedformationEnergyFits.append(self._getFit( xdata, naUnrelaxedFormationEnergyBySize[j:(j + cnt)], orders )[0]) formationEnergyFits.append(self._getFit( xdata, naFormationEnergyBySize[j:(j + cnt)], orders )[0]) relaxationVolumeFits.append(self._getFit( xdata, naRelaxationVolumeBySize[j:(j + cnt)], orders )[0]) unrelaxedformationEnergyFitsBySize.append(unrelaxedformationEnergyFits) formationEnergyFitsBySize.append(formationEnergyFits) relaxationVolumeFitsBySize.append(relaxationVolumeFits) # Output Fitting Results print('\n[Fitting Results Summary]') print('Sizes:', sizes) print('Data Points Used:', FITS_CNT) print('Orders Included:\n', FITS_ORDERS) print('Unrelaxed Formation Energy Fits By Size:\n', unrelaxedformationEnergyFitsBySize) print('Formation Energy Fits By Size:\n', formationEnergyFitsBySize) print('Relaxation Volume Fits By Size:\n', relaxationVolumeFitsBySize) # Obtain Extrapolated Value and Uncertainty unrelaxedformationEnergy, unrelaxedformationEnergyUncert = self._getValueUncert( FITS_VFE_VALUE, FITS_VFE_UNCERT, # FMIN_FTOL * formationEnergyBySize[-1], self.VFEUncert, 2, unrelaxedformationEnergyFitsBySize, ) formationEnergy, formationEnergyUncert = self._getValueUncert( FITS_VFE_VALUE, FITS_VFE_UNCERT, # FMIN_FTOL * formationEnergyBySize[-1], self.VFEUncert, 2, formationEnergyFitsBySize, ) relaxationVolume, relaxationVolumeUncert = self._getValueUncert( FITS_VRV_VALUE, FITS_VRV_UNCERT, # FMIN_XTOL * (self.latticeConsts[0] * CELL_SIZE_MAX)**3, self.VRVUncert, 2, relaxationVolumeFitsBySize, ) # Construct Results Dictionary unrelaxedformationEnergyResult = OrderedDict([ ('property-id', VUFE_PROP_ID), ('instance-id', 1), ('unrelaxed-formation-potential-energy', V(unrelaxedformationEnergy, UNIT_ENERGY, unrelaxedformationEnergyUncert)), ]) formationEnergyResult = OrderedDict([ ('property-id', VFE_PROP_ID), ('instance-id', 2), ('relaxed-formation-potential-energy', V(formationEnergy, UNIT_ENERGY, formationEnergyUncert)), ]) relaxationVolumeResult = OrderedDict([ ('property-id', VRV_PROP_ID), ('instance-id', 3), ('relaxation-volume', V(relaxationVolume, UNIT_VOLUME, relaxationVolumeUncert)), ]) hostInfo = OrderedDict([ ('host-cauchy-stress', V([0, 0, 0, 0, 0, 0], UNIT_PRESSURE)), ('host-removed-atom', V(1)), ('host-short-name', V([self.lattice])), ('host-a', V(np.linalg.norm(unitCell[0]), UNIT_LENGTH)), ('host-b', V(np.linalg.norm(unitCell[1]), UNIT_LENGTH)), ('host-c', V(np.linalg.norm(unitCell[2]), UNIT_LENGTH)), ('host-alpha', V(self._getAngle(unitCell[1], unitCell[2]), UNIT_ANGLE)), ('host-beta', V(self._getAngle(unitCell[2], unitCell[0]), UNIT_ANGLE)), ('host-gamma', V(self._getAngle(unitCell[0], unitCell[1]), UNIT_ANGLE)), ('host-space-group', V(SPACE_GROUPS[self.lattice])), ('host-wyckoff-multiplicity-and-letter', V(WYCKOFF_CODES[self.lattice])), ('host-wyckoff-coordinates', V(WYCKOFF_SITES[self.lattice])), ('host-wyckoff-species', V([self.elem] * len(WYCKOFF_CODES[self.lattice]))), ]) reservoirInfo = OrderedDict([ ('reservoir-cohesive-potential-energy', V(-unitBulk.get_potential_energy()/unitBulk.get_number_of_atoms(), UNIT_ENERGY)), ('reservoir-short-name', V([self.lattice])), ('reservoir-cauchy-stress', V([0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UNIT_PRESSURE)), ('reservoir-a', V(np.linalg.norm(unitCell[0]), UNIT_LENGTH)), ('reservoir-b', V(np.linalg.norm(unitCell[1]), UNIT_LENGTH)), ('reservoir-c', V(np.linalg.norm(unitCell[2]), UNIT_LENGTH)), ('reservoir-alpha', V(self._getAngle(unitCell[1], unitCell[2]), UNIT_ANGLE)), ('reservoir-beta', V(self._getAngle(unitCell[2], unitCell[0]), UNIT_ANGLE)), ('reservoir-gamma', V(self._getAngle(unitCell[0], unitCell[1]), UNIT_ANGLE)), ('reservoir-space-group', V(SPACE_GROUPS[self.lattice])), ('reservoir-wyckoff-multiplicity-and-letter', V(WYCKOFF_CODES[self.lattice])), ('reservoir-wyckoff-coordinates', V(WYCKOFF_SITES[self.lattice])), ('reservoir-wyckoff-species', V([self.elem] * len(WYCKOFF_CODES[self.lattice]))), ]) unrelaxedformationEnergyResult.update(hostInfo) unrelaxedformationEnergyResult.update(reservoirInfo) formationEnergyResult.update(hostInfo) formationEnergyResult.update(reservoirInfo) relaxationVolumeResult.update(hostInfo) results = [unrelaxedformationEnergyResult, formationEnergyResult, relaxationVolumeResult] return results def V(value, unit = '', uncert = ''): # Generate OrderedDict for JSON Dump res = OrderedDict([ (KEY_SOURCE_VALUE, value), ]) if unit != '': res.update(OrderedDict([ (KEY_SOURCE_UNIT, unit), ])) if uncert != '': res.update(OrderedDict([ (KEY_SOURCE_UNCERT, uncert) ])) return res def main(): # Create Environment and Common Variables # If it's hcp crystal, process both a and c # Otherwise, process a, ignore c # Obtain Inputs model = input("Enter the name of the KIM Model you wish to perform calculations for:\n") elem = input("Enter the name of the species you wish to simulate:\n") lattice = input("Enter the lattice type of the crystal ('bcc', 'diamond', 'fcc', 'hcp', or 'sc'):\n") latticeConstA = input("Enter the lattice constant 'a' in meters:\n") latticeConstC = input("If the lattice type is hcp, enter the lattice constant 'c' in meters (if you are not using hcp as the lattice, simply put any value here, as it will be ignored):\n") # Cast lattice to lower case lattice = lattice.lower() # Convert the lattice constant values from the query from meters to Angstroms if lattice == 'hcp': latticeConstA = float(latticeConstA)*1e10 latticeConstC = float(latticeConstC)*1e10 latticeConsts = [float(latticeConstA), float(latticeConstC)] else: latticeConstA = float(latticeConstA)*1e10 latticeConsts = [float(latticeConstA)] calc = KIM(model) print('Parameters Input:', elem, lattice, model, latticeConsts) # Create Instance instance = VacancyCalculation(calc, elem, model, lattice, latticeConsts) # Obtain Results res = instance.getResults() # Output results print('\n[Final Results]') resStr = json.dumps(res, separators = (' ',' '), indent = 4) print(resStr) with open(os.path.abspath('output/results.edn'), 'w') as f: f.write(resStr) if __name__ == '__main__': main()