Source code for MMTK.ForceFields.ForceFieldTest
# This module implements test functions.
#
# Written by Konrad Hinsen
#
"""
Force field consistency tests
"""
__docformat__ = 'restructuredtext'
from MMTK import Utility
from Scientific.Geometry import Vector, ex, ey, ez
from Scientific import N
import itertools
#
# Check consistency of energies and gradients.
#
[docs]def gradientTest(universe, atoms = None, delta = 0.0001):
"""
Test gradients by comparing to numerical derivatives of the energy.
:param universe: the universe on which the test is performed
:type universe: :class:`~MMTK.Universe.Universe`
:param atoms: the atoms of the universe for which the gradient
is tested (default: all atoms)
:type atoms: list
:param delta: the step size used in calculating the numerical derivatives
:type delta: float
"""
e0, grad = universe.energyAndGradients()
print 'Energy: ', e0
if atoms is None:
atoms = universe.atomList()
for a in atoms:
print a
print grad[a]
num_grad = []
for v in [ex, ey, ez]:
x = a.position()
a.setPosition(x+delta*v)
eplus = universe.energy()
a.setPosition(x-delta*v)
eminus = universe.energy()
a.setPosition(x)
num_grad.append(0.5*(eplus-eminus)/delta)
print Vector(num_grad)
#
# Check consistency of gradients and force constants.
#
[docs]def forceConstantTest(universe, atoms = None, delta = 0.0001):
"""
Test force constants by comparing to the numerical derivatives
of the gradients.
:param universe: the universe on which the test is performed
:type universe: :class:`~MMTK.Universe.Universe`
:param atoms: the atoms of the universe for which the gradient
is tested (default: all atoms)
:type atoms: list
:param delta: the step size used in calculating the numerical derivatives
:type delta: float
"""
e0, grad0, fc = universe.energyGradientsAndForceConstants()
if atoms is None:
atoms = universe.atomList()
for a1, a2 in itertools.chain(itertools.izip(atoms, atoms),
Utility.pairs(atoms)):
print a1, a2
print fc[a1, a2]
num_fc = []
for v in [ex, ey, ez]:
x = a1.position()
a1.setPosition(x+delta*v)
e_plus, grad_plus = universe.energyAndGradients()
a1.setPosition(x-delta*v)
e_minus, grad_minus = universe.energyAndGradients()
a1.setPosition(x)
num_fc.append(0.5*(grad_plus[a2]-grad_minus[a2])/delta)
print N.array(map(lambda a: a.array, num_fc))
#
# Check consistency of gradients and virial
#
[docs]def virialTest(universe):
"""
Test the virial by comparing to an explicit computation from
positions and gradients.
:param universe: the universe on which the test is performed
:type universe: :class:`~MMTK.Universe.Universe`
"""
ev = universe.energyEvaluator()
e, grad = ev(gradients = True)
virial = ev.lastVirial()
conf = universe.configuration()
print virial, -(conf*grad).sumOverParticles()
if __name__ == '__main__':
from MMTK import *
from MMTK.ForceFields import Amber94ForceField
delta = 0.001
world = InfiniteUniverse(Amber94ForceField())
m = Molecule('water')
m.O.translateBy(Vector(0.,0.,0.01))
m.H1.translateBy(Vector(0.01,0.,0.))
atoms = None
world.addObject(m)
gradientTest(world, atoms, delta)
forceConstantTest(world, atoms, delta)
virialTest(world)