#!/usr/bin/env python """ Lattice constant Test Driver Computes the lattice constant for any material and any cubic crystal structure by minimizing the energy using a simplex minimization Date: 2012/09/03 Author: Alex Alemi & Matt Bierbaum Last Update: 2016/02/26 Junhao Li 2017/08/21 Ellad Tadmor 2018/02/17 Dan Karls 2019/02/24 Dan Karls (added logic to check if lattice is unstable, handle very small lattice constants for Toy Models) """ from ase.lattice.cubic import BodyCenteredCubic, FaceCenteredCubic, Diamond, SimpleCubic import scipy.optimize as opt from ase.calculators.kim.kim import KIM import numpy as np import jinja2 import os import json import sys def create_cubic_lattice(symbol, lattice, a, ncells_per_side): """Note that we use ase.lattice instead of ase.build.bulk so that we can construct finite clusters of unit cells if necessary.""" dimensions = (ncells_per_side, ncells_per_side, ncells_per_side) if lattice == "bcc": slab = BodyCenteredCubic(symbol, latticeconstant=a, size=dimensions) elif lattice == "fcc": slab = FaceCenteredCubic(symbol, latticeconstant=a, size=dimensions) elif lattice == "diamond": slab = Diamond(symbol, latticeconstant=a, size=dimensions) elif lattice == "sc": slab = SimpleCubic(symbol, latticeconstant=a, size=dimensions) return slab def energy(a, slab, cell, positions): """ Compute the energy of a lattice given the lattice constant. Multiply the positions for every lattice type since in this test only diamond has multiple sites, all other positions are zero. """ slab.set_positions( positions * a ) slab.set_cell( cell * a ) energy = slab.get_potential_energy() return energy def rescale_to_get_nonzero_forces(atoms, ftol): ''' If the given configuration has forces smaller than 'ftol' (presumably because the distance between atoms is too large), rescale it making it smaller until the maximum force is nonzero. In a perfect crystal, the crystal is rescaled until the atoms on the surface reach the minimum value (internal atoms will have zero force). Note that any periodicity is turned off for the rescaling and then restored at the end. ''' if len(atoms)<2: raise Exception('ERROR: Invalid configuration. Must have at least 2 atoms. Number of atoms = {}'.format(len(atoms))) # Temporarily turn off any periodicity pbc_save = atoms.get_pbc() cell = atoms.get_cell() atoms.set_pbc([False, False, False]) # Rescale cell and atoms forces = atoms.get_forces() fmax = max(abs(forces.min()), abs(forces.max())) # find max in abs value if fmax np.finfo(delpmin).tiny: atoms.positions *= 0.5 # make configuration half the size cell *= 0.5 # make cell half the size delpmin *= 0.5 forces = atoms.get_forces() # get max force fmax = max(abs(forces.min()), abs(forces.max())) if fmax >= ftol: # Restore periodicity atoms.set_pbc(pbc_save) atoms.set_cell(cell) return # success! raise Exception('ERROR: Unable to scale configuration down to nonzero forces.') def get_lattice_constant(model, symbol, lattice): """ Here, we start with an initial guess of 2.5 Angstroms for the lattice constant. We then attempt to minimize the potential energy using this initial value and, if the optimization fails for any reason, we increment the initial guess by 0.5 Angstroms and repeat the calculation. This process repeats until either the calculation succeeds or until a maximum initial guess is exceeded. """ calc = KIM(model) repeat = 0 slab = create_cubic_lattice(symbol, lattice, a=1, ncells_per_side=1) cell = slab.get_cell() positions = slab.get_positions() slab.set_calculator(calc) particles = len(slab) init_guess = 2.5 max_lattice_constant = 10.0 # First, compute the isolated energy of the cell of atoms. This is assumed # to be equivalent to the energy at a lattice constant of 1000 Angstroms energy_at_inf = energy(1000.0, slab, cell, positions) # Now check if the energy at a=init_guess is approximately equal to the isolated energy. # This is unlikely, but can occur for Toy Models. On the other hand, many # non-toy Models may produce an error because init_guess is too small for # them to handle and, in this case, we simply proceed to the loop below, # which will catch the associated exceptions as it tries larger initial # guesses. init_guess_too_large = False try: energy_at_init_guess = energy(init_guess, slab, cell, positions) if abs(energy_at_init_guess - energy_at_inf) < 1e-8: init_guess_too_large = True except: pass if not init_guess_too_large: # Start at init_guess and increment up to max_lattice_constant before giving up done = False while not done: print("\nAttempting to perform relaxation using initial lattice constant guess of {} Angstroms\n".format(init_guess, '.1f')) try: aopt_arr, eopt, iterations, funccalls, warnflag = opt.fmin(energy, init_guess, args=(slab,cell,positions), full_output=True, xtol=1e-8) done = True except: print("\nFailed to perform relaxation using initial lattice constant guess of {} Angstroms\n".format(init_guess, '.1f')) init_guess += 0.5 if init_guess > max_lattice_constant: done = True raise Exception("Exceeded maximum of {} for initial lattice constant guess for relaxation".format(max_lattice_constant, '.1f')) else: print("\nDetected that energy at a lattice constant of {} Angstroms was equal to isolated energy (i.e. no interactions). Shrinking box until significant forces are present...\n".format(init_guess, '.1f')) # In this case, we construct a 2x2x2 cell, turn off periodic boundary # and start shrinking the box until we get non-zero forces. # This lattice constant is then used as the initial guess. ncells_per_side=2 dummylattice = create_cubic_lattice(symbol, lattice, a=init_guess, ncells_per_side=ncells_per_side) dummylattice.set_calculator(calc) rescale_to_get_nonzero_forces(dummylattice, ftol=0.01) new_init_guess = dummylattice.get_cell()[0,0]/ncells_per_side min_lattice_constant = 0.5*new_init_guess decrement_value = 0.1*new_init_guess # Try new_init_guess, 0.9*new_init_guess, 0.8*new_init_guess, ..., min_lattice_constant before giving up done = False while not done: print("\nAttempting to perform relaxation using initial lattice constant guess of {} Angstroms\n".format(new_init_guess, '.4f')) try: aopt_arr, eopt, iterations, funccalls, warnflag = opt.fmin(energy, new_init_guess, args=(slab,cell,positions), full_output=True, xtol=1e-8) done = True except: print("\nFailed to perform relaxation using initial lattice constant guess of {} Angstroms\n".format(init_guess, '.4f')) new_init_guess -= decrement_value if new_init_guess < min_lattice_constant: done = True raise Exception("Fell below minimum of {} for initial lattice constant guess for relaxation.".format(min_lattice_constant, '.1f')) # Check to see if the lattice constant we converged to is actually just # the one large enough so that particles are not interacting, i.e. this # lattice is unstable if abs(eopt-energy_at_inf) < 1e-8: raise Exception("Calculation converged to lattice constant ({}) for which particles do not interact. " "This may indicate the lattice is unstable.".format(aopt_arr[0])) info = {"iterations": iterations, "func_calls": funccalls, "warnflag": warnflag, "repeat": repeat} return aopt_arr[0], -eopt/particles, info if __name__ == "__main__": #grab from stdin (or a file) symbol = raw_input("element=") lattice = raw_input("lattice type=") model = raw_input("modelname=") print symbol, lattice, model aopt, eopt, info = get_lattice_constant(model=model, symbol=symbol, lattice=lattice) space_groups = {"fcc": "Fm-3m", "bcc": "Im-3m", "sc": "Pm-3m", "diamond": "Fd-3m"} wyckoff_codes = {"fcc": "4a", "bcc": "2a", "sc": "1a", "diamond": "8a"} dummylattice = create_cubic_lattice(symbol, lattice, a=1, ncells_per_side=1) normed_basis = { lattice: json.dumps(dummylattice.positions.tolist(), separators=(' ', ' ')) for lattice in space_groups.keys() } count = dummylattice.positions.shape[0] # count atoms in a unit cell # pack the results in a dictionary results = {"lattice_constant": aopt, "cohesive_energy": eopt, "element": symbol, "species": '" "'.join([symbol] * count), # repeat symbols to match normed_basis "crystal_structure": lattice, "space_group": space_groups[lattice], "wyckoff_code": wyckoff_codes[lattice], "basis_atoms": normed_basis[lattice]} results.update(info) print results template_environment = jinja2.Environment( loader=jinja2.FileSystemLoader('/'), block_start_string='@[', block_end_string=']@', variable_start_string='@<', variable_end_string='>@', comment_start_string='@#', comment_end_string='#@', undefined=jinja2.StrictUndefined, ) #template the EDN output with open(os.path.abspath("output/results.edn"), "w") as f: template = template_environment.get_template(os.path.abspath("results.edn.tpl")) f.write(template.render(**results))