Source code for ase.calculators.lj

from __future__ import division

import numpy as np

from ase.neighborlist import NeighborList
from ase.calculators.calculator import Calculator, all_changes
from ase.calculators.calculator import PropertyNotImplementedError


[docs]class LennardJones(Calculator): implemented_properties = ['energy', 'forces', 'stress'] default_parameters = {'epsilon': 1.0, 'sigma': 1.0, 'rc': None} nolabel = True def __init__(self, **kwargs): Calculator.__init__(self, **kwargs) def calculate(self, atoms=None, properties=['energy'], system_changes=all_changes): Calculator.calculate(self, atoms, properties, system_changes) natoms = len(self.atoms) sigma = self.parameters.sigma epsilon = self.parameters.epsilon rc = self.parameters.rc if rc is None: rc = 3 * sigma if 'numbers' in system_changes: self.nl = NeighborList([rc / 2] * natoms, self_interaction=False) self.nl.update(self.atoms) positions = self.atoms.positions cell = self.atoms.cell e0 = 4 * epsilon * ((sigma / rc)**12 - (sigma / rc)**6) energy = 0.0 forces = np.zeros((natoms, 3)) stress = np.zeros((3, 3)) for a1 in range(natoms): neighbors, offsets = self.nl.get_neighbors(a1) cells = np.dot(offsets, cell) d = positions[neighbors] + cells - positions[a1] r2 = (d**2).sum(1) c6 = (sigma**2 / r2)**3 c6[r2 > rc**2] = 0.0 energy -= e0 * (c6 != 0.0).sum() c12 = c6**2 energy += 4 * epsilon * (c12 - c6).sum() f = (24 * epsilon * (2 * c12 - c6) / r2)[:, np.newaxis] * d forces[a1] -= f.sum(axis=0) for a2, f2 in zip(neighbors, f): forces[a2] += f2 stress += np.dot(f.T, d) if 'stress' in properties: if self.atoms.number_of_lattice_vectors == 3: stress += stress.T.copy() stress *= -0.5 / self.atoms.get_volume() self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]] else: raise PropertyNotImplementedError self.results['energy'] = energy self.results['free_energy'] = energy self.results['forces'] = forces