Source code for moldesign.models.base

# Copyright 2016 Autodesk Inc.
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.
import itertools
import numpy as np

import moldesign as mdt
from moldesign import units as u
from moldesign.method import Method

[docs]class EnergyModelBase(Method): """ Base class for all energy models """ # TODO:should add some architecture to check implementations, e.g. # TODO: ensure prep gets called when necessary, make sure all parameters can be consumed DEFAULT_PROPERTIES = ['potential_energy', 'forces'] """List[str]: list of the properties that are always calculated by this method""" ALL_PROPERTIES = DEFAULT_PROPERTIES """List[str]: List of all the properties that this model can calculate""" PARAMETERS = []
[docs] def calculate(self, requests): """Calculate the the default properties and any additiona requests Arguments: requests (List[str]): the requested properties to calculate Returns: utils.DotDict: A dict of calculated properties (or a job object that will return them) """ self.prep() raise NotImplementedError('EnergyModelBase is an abstract base class')
[docs] def minimize(self, method='bfgs', **kwargs): """ If the energy model provides its own minimizer, it should be hooked up here """ raise NotImplementedError()
[docs] def get_formal_charge(self): """Determine the formal charge of the molecular system. This can be set either as a molecular attribute OR in the parameters of the energy model. Returns: u.Scalar[charge]: the formal charge used for this model """ if 'charge' in self.params and self.params.charge is not None: return self.params.charge.value_in(u.q_e) elif hasattr(self.mol, 'charge'): return self.mol.charge.value_in(u.q_e) else: return 0
[docs] def finite_difference_force(self, direction=0, stepsize=0.025 * u.angstrom): """ Compute force using a finite difference with the given step size. Args: direction (int): EITHER +1, -1, (for one-sided finite differences) or 0 (for central difference - better but twice as expensive) step (u.Scalar[lenght]): step size to take in each direction Returns: u.Vector[force]: force vector, len= `self.mol.ndims` """ # TODO: this should totally be parallelized - how do we request/configure/control this? forces = np.zeros(self.mol.ndims) * u.default.force properties = [] if direction in (-1, 1): e0 = self.mol.calc_potential_energy() else: assert direction == 0, 'Finite difference direction must be -1, 0, or 1' for iatom, idim in itertools.product(xrange(self.mol.num_atoms), xrange(3)): print '\rFinite differencing %s for atom %d/%d'%('xyz'[iatom], iatom+1, self.mol.num_atoms), if direction == 0: self.mol.positions[iatom, idim] += stepsize / 2.0 eplus = self.mol.calc_potential_energy() pplus = self.mol.positions[iatom, idim] -= stepsize eminus = self.mol.calc_potential_energy() pminus = self.mol.positions[iatom, idim] += stepsize / 2.0 # resets the position properties.append((pminus, pplus)) forces[iatom, idim] = (eminus - eplus) / stepsize elif direction in (-1, 1): self.mol.positions[iatom, idim] += direction * stepsize enew = self.mol.calc_potential_energy() properties.append( forces[iatom, idim] = (e0 - enew) / (direction * stepsize) self.mol.positions[iatom, idim] -= direction * stepsize return forces, properties
[docs] def prep(self): """ Prepare to run. Possibly do a test to ensure that the model is ready. """ raise NotImplementedError('EnergyModelBase is an abstract base class')
[docs]class MMBase(EnergyModelBase): """Common interface for molecular mechanics""" PARAMETERS = (EnergyModelBase.PARAMETERS + mdt.parameters.mm_model_parameters.values()) def __init__(self, *args, **kwargs): super(MMBase, self).__init__(*args, **kwargs) self.mdtforcefield = None
[docs]class QMBase(EnergyModelBase): """Common interface for quantum mechanics""" PARAMETERS = mdt.parameters.qm_model_parameters.values() DEFAULT_PROPERTIES = ['potential_energy', 'nuclear_repulsion', 'dipole_moment', 'orbitals', 'orbital_energies'] ALL_PROPERTIES = DEFAULT_PROPERTIES # properties will be a pretty long list for most packages
[docs] def set_wfn_guess(self): raise NotImplementedError
[docs]class QMMMBase(EnergyModelBase): DEFAULT_PROPERTIES = ['potential_energy', 'qm_energy', 'mm_energy', 'interaction_energy' 'qm_dipole_moment', 'orbitals', 'orbital_energies'] ALL_PROPERTIES = DEFAULT_PROPERTIES PARAMETERS = MMBase.PARAMETERS + QMBase.PARAMETERS