Commit c9fcabb5 authored by Peter Eastman's avatar Peter Eastman
Browse files

More updates to dispersion PME

parent c274f08d
...@@ -412,7 +412,8 @@ public: ...@@ -412,7 +412,8 @@ public:
bool usesPeriodicBoundaryConditions() const { bool usesPeriodicBoundaryConditions() const {
return nonbondedMethod == NonbondedForce::CutoffPeriodic || return nonbondedMethod == NonbondedForce::CutoffPeriodic ||
nonbondedMethod == NonbondedForce::Ewald || nonbondedMethod == NonbondedForce::Ewald ||
nonbondedMethod == NonbondedForce::PME; nonbondedMethod == NonbondedForce::PME ||
nonbondedMethod == NonbondedForce::LJPME;
} }
protected: protected:
ForceImpl* createImpl() const; ForceImpl* createImpl() const;
......
...@@ -90,7 +90,7 @@ void NonbondedForceImpl::initialize(ContextImpl& context) { ...@@ -90,7 +90,7 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
exceptions[particle1].insert(particle2); exceptions[particle1].insert(particle2);
exceptions[particle2].insert(particle1); exceptions[particle2].insert(particle1);
} }
if (owner.getNonbondedMethod() != NonbondedForce::NoCutoff) { if (owner.getNonbondedMethod() != NonbondedForce::NoCutoff && owner.getNonbondedMethod() != NonbondedForce::CutoffNonPeriodic) {
Vec3 boxVectors[3]; Vec3 boxVectors[3];
system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]); system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
double cutoff = owner.getCutoffDistance(); double cutoff = owner.getCutoffDistance();
......
...@@ -42,7 +42,7 @@ NonbondedForceProxy::NonbondedForceProxy() : SerializationProxy("NonbondedForce" ...@@ -42,7 +42,7 @@ NonbondedForceProxy::NonbondedForceProxy() : SerializationProxy("NonbondedForce"
} }
void NonbondedForceProxy::serialize(const void* object, SerializationNode& node) const { void NonbondedForceProxy::serialize(const void* object, SerializationNode& node) const {
node.setIntProperty("version", 1); node.setIntProperty("version", 2);
const NonbondedForce& force = *reinterpret_cast<const NonbondedForce*>(object); const NonbondedForce& force = *reinterpret_cast<const NonbondedForce*>(object);
node.setIntProperty("forceGroup", force.getForceGroup()); node.setIntProperty("forceGroup", force.getForceGroup());
node.setIntProperty("method", (int) force.getNonbondedMethod()); node.setIntProperty("method", (int) force.getNonbondedMethod());
...@@ -59,6 +59,11 @@ void NonbondedForceProxy::serialize(const void* object, SerializationNode& node) ...@@ -59,6 +59,11 @@ void NonbondedForceProxy::serialize(const void* object, SerializationNode& node)
node.setIntProperty("nx", nx); node.setIntProperty("nx", nx);
node.setIntProperty("ny", ny); node.setIntProperty("ny", ny);
node.setIntProperty("nz", nz); node.setIntProperty("nz", nz);
force.getLJPMEParameters(alpha, nx, ny, nz);
node.setDoubleProperty("ljAlpha", alpha);
node.setIntProperty("ljnx", nx);
node.setIntProperty("ljny", ny);
node.setIntProperty("ljnz", nz);
node.setIntProperty("recipForceGroup", force.getReciprocalSpaceForceGroup()); node.setIntProperty("recipForceGroup", force.getReciprocalSpaceForceGroup());
SerializationNode& particles = node.createChildNode("Particles"); SerializationNode& particles = node.createChildNode("Particles");
for (int i = 0; i < force.getNumParticles(); i++) { for (int i = 0; i < force.getNumParticles(); i++) {
...@@ -76,7 +81,8 @@ void NonbondedForceProxy::serialize(const void* object, SerializationNode& node) ...@@ -76,7 +81,8 @@ void NonbondedForceProxy::serialize(const void* object, SerializationNode& node)
} }
void* NonbondedForceProxy::deserialize(const SerializationNode& node) const { void* NonbondedForceProxy::deserialize(const SerializationNode& node) const {
if (node.getIntProperty("version") != 1) int version = node.getIntProperty("version");
if (version < 1 || version > 2)
throw OpenMMException("Unsupported version number"); throw OpenMMException("Unsupported version number");
NonbondedForce* force = new NonbondedForce(); NonbondedForce* force = new NonbondedForce();
try { try {
...@@ -93,6 +99,13 @@ void* NonbondedForceProxy::deserialize(const SerializationNode& node) const { ...@@ -93,6 +99,13 @@ void* NonbondedForceProxy::deserialize(const SerializationNode& node) const {
int ny = node.getIntProperty("ny", 0); int ny = node.getIntProperty("ny", 0);
int nz = node.getIntProperty("nz", 0); int nz = node.getIntProperty("nz", 0);
force->setPMEParameters(alpha, nx, ny, nz); force->setPMEParameters(alpha, nx, ny, nz);
if (version >= 2) {
alpha = node.getDoubleProperty("ljAlpha", 0.0);
nx = node.getIntProperty("ljnx", 0);
ny = node.getIntProperty("ljny", 0);
nz = node.getIntProperty("ljnz", 0);
force->setLJPMEParameters(alpha, nx, ny, nz);
}
force->setReciprocalSpaceForceGroup(node.getIntProperty("recipForceGroup", -1)); force->setReciprocalSpaceForceGroup(node.getIntProperty("recipForceGroup", -1));
const SerializationNode& particles = node.getChildNode("Particles"); const SerializationNode& particles = node.getChildNode("Particles");
for (int i = 0; i < (int) particles.getChildren().size(); i++) { for (int i = 0; i < (int) particles.getChildren().size(); i++) {
......
...@@ -53,6 +53,9 @@ void testSerialization() { ...@@ -53,6 +53,9 @@ void testSerialization() {
double alpha = 0.5; double alpha = 0.5;
int nx = 3, ny = 5, nz = 7; int nx = 3, ny = 5, nz = 7;
force.setPMEParameters(alpha, nx, ny, nz); force.setPMEParameters(alpha, nx, ny, nz);
double dalpha = 0.8;
int dnx = 4, dny = 6, dnz = 7;
force.setLJPMEParameters(dalpha, dnx, dny, dnz);
force.addParticle(1, 0.1, 0.01); force.addParticle(1, 0.1, 0.01);
force.addParticle(0.5, 0.2, 0.02); force.addParticle(0.5, 0.2, 0.02);
force.addParticle(-0.5, 0.3, 0.03); force.addParticle(-0.5, 0.3, 0.03);
...@@ -84,6 +87,13 @@ void testSerialization() { ...@@ -84,6 +87,13 @@ void testSerialization() {
ASSERT_EQUAL(nx, nx2); ASSERT_EQUAL(nx, nx2);
ASSERT_EQUAL(ny, ny2); ASSERT_EQUAL(ny, ny2);
ASSERT_EQUAL(nz, nz2); ASSERT_EQUAL(nz, nz2);
double dalpha2;
int dnx2, dny2, dnz2;
force2.getLJPMEParameters(dalpha2, dnx2, dny2, dnz2);
ASSERT_EQUAL(dalpha, dalpha2);
ASSERT_EQUAL(dnx, dnx2);
ASSERT_EQUAL(dny, dny2);
ASSERT_EQUAL(dnz, dnz2);
for (int i = 0; i < force.getNumParticles(); i++) { for (int i = 0; i < force.getNumParticles(); i++) {
double charge1, sigma1, epsilon1; double charge1, sigma1, epsilon1;
double charge2, sigma2, epsilon2; double charge2, sigma2, epsilon2;
......
...@@ -40,6 +40,7 @@ CutoffNonPeriodic = forcefield.CutoffNonPeriodic ...@@ -40,6 +40,7 @@ CutoffNonPeriodic = forcefield.CutoffNonPeriodic
CutoffPeriodic = forcefield.CutoffPeriodic CutoffPeriodic = forcefield.CutoffPeriodic
Ewald = forcefield.Ewald Ewald = forcefield.Ewald
PME = forcefield.PME PME = forcefield.PME
LJPME = forcefield.LJPME
HBonds = forcefield.HBonds HBonds = forcefield.HBonds
AllBonds = forcefield.AllBonds AllBonds = forcefield.AllBonds
......
...@@ -169,7 +169,7 @@ class AmberPrmtopFile(object): ...@@ -169,7 +169,7 @@ class AmberPrmtopFile(object):
---------- ----------
nonbondedMethod : object=NoCutoff nonbondedMethod : object=NoCutoff
The method to use for nonbonded interactions. Allowed values are The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, or LJPME.
nonbondedCutoff : distance=1*nanometer nonbondedCutoff : distance=1*nanometer
The cutoff distance to use for nonbonded interactions The cutoff distance to use for nonbonded interactions
constraints : object=None constraints : object=None
...@@ -202,7 +202,7 @@ class AmberPrmtopFile(object): ...@@ -202,7 +202,7 @@ class AmberPrmtopFile(object):
added to a hydrogen is subtracted from the heavy atom to keep their added to a hydrogen is subtracted from the heavy atom to keep their
total mass the same. total mass the same.
ewaldErrorTolerance : float=0.0005 ewaldErrorTolerance : float=0.0005
The error tolerance to use if nonbondedMethod is Ewald or PME. The error tolerance to use if nonbondedMethod is Ewald, PME, or LJPME.
switchDistance : float=0*nanometers switchDistance : float=0*nanometers
The distance at which the potential energy switching function is The distance at which the potential energy switching function is
turned on for Lennard-Jones interactions. If the switchDistance is 0 turned on for Lennard-Jones interactions. If the switchDistance is 0
...@@ -222,10 +222,11 @@ class AmberPrmtopFile(object): ...@@ -222,10 +222,11 @@ class AmberPrmtopFile(object):
ff.CutoffNonPeriodic:'CutoffNonPeriodic', ff.CutoffNonPeriodic:'CutoffNonPeriodic',
ff.CutoffPeriodic:'CutoffPeriodic', ff.CutoffPeriodic:'CutoffPeriodic',
ff.Ewald:'Ewald', ff.Ewald:'Ewald',
ff.PME:'PME'} ff.PME:'PME',
ff.LJPME:'LJPME'}
if nonbondedMethod not in methodMap: if nonbondedMethod not in methodMap:
raise ValueError('Illegal value for nonbonded method') raise ValueError('Illegal value for nonbonded method')
if not self._prmtop.getIfBox() and nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME): if not self._prmtop.getIfBox() and nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME, ff.LJPME):
raise ValueError('Illegal nonbonded method for a non-periodic system') raise ValueError('Illegal nonbonded method for a non-periodic system')
constraintMap = {None:None, constraintMap = {None:None,
ff.HBonds:'h-bonds', ff.HBonds:'h-bonds',
......
...@@ -690,7 +690,7 @@ class CharmmPsfFile(object): ...@@ -690,7 +690,7 @@ class CharmmPsfFile(object):
The parameter set to use to parametrize this molecule The parameter set to use to parametrize this molecule
nonbondedMethod : object=NoCutoff nonbondedMethod : object=NoCutoff
The method to use for nonbonded interactions. Allowed values are The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, or LJPME.
nonbondedCutoff : distance=1*nanometer nonbondedCutoff : distance=1*nanometer
The cutoff distance to use for nonbonded interactions. The cutoff distance to use for nonbonded interactions.
switchDistance : distance=0*nanometer switchDistance : distance=0*nanometer
...@@ -728,7 +728,7 @@ class CharmmPsfFile(object): ...@@ -728,7 +728,7 @@ class CharmmPsfFile(object):
added to a hydrogen is subtracted from the heavy atom to keep their added to a hydrogen is subtracted from the heavy atom to keep their
total mass the same. total mass the same.
ewaldErrorTolerance : float=0.0005 ewaldErrorTolerance : float=0.0005
The error tolerance to use if the nonbonded method is Ewald or PME. The error tolerance to use if the nonbonded method is Ewald, PME, or LJPME.
flexibleConstraints : bool=True flexibleConstraints : bool=True
Are our constraints flexible or not? Are our constraints flexible or not?
verbose : bool=False verbose : bool=False
...@@ -746,10 +746,10 @@ class CharmmPsfFile(object): ...@@ -746,10 +746,10 @@ class CharmmPsfFile(object):
cutoff = cutoff.value_in_unit(u.nanometers) cutoff = cutoff.value_in_unit(u.nanometers)
if nonbondedMethod not in (ff.NoCutoff, ff.CutoffNonPeriodic, if nonbondedMethod not in (ff.NoCutoff, ff.CutoffNonPeriodic,
ff.CutoffPeriodic, ff.Ewald, ff.PME): ff.CutoffPeriodic, ff.Ewald, ff.PME, ff.LJPME):
raise ValueError('Illegal value for nonbonded method') raise ValueError('Illegal value for nonbonded method')
if not hasbox and nonbondedMethod in (ff.CutoffPeriodic, if not hasbox and nonbondedMethod in (ff.CutoffPeriodic,
ff.Ewald, ff.PME): ff.Ewald, ff.PME, ff.LJPME):
raise ValueError('Illegal nonbonded method for a ' raise ValueError('Illegal nonbonded method for a '
'non-periodic system') 'non-periodic system')
if implicitSolvent not in (HCT, OBC1, OBC2, GBn, GBn2, None): if implicitSolvent not in (HCT, OBC1, OBC2, GBn, GBn2, None):
...@@ -1009,6 +1009,8 @@ class CharmmPsfFile(object): ...@@ -1009,6 +1009,8 @@ class CharmmPsfFile(object):
force.setNonbondedMethod(mm.NonbondedForce.Ewald) force.setNonbondedMethod(mm.NonbondedForce.Ewald)
elif nonbondedMethod is ff.PME: elif nonbondedMethod is ff.PME:
force.setNonbondedMethod(mm.NonbondedForce.PME) force.setNonbondedMethod(mm.NonbondedForce.PME)
elif nonbondedMethod is ff.LJPME:
force.setNonbondedMethod(mm.NonbondedForce.LJPME)
else: else:
raise ValueError('Cutoff method is not understood') raise ValueError('Cutoff method is not understood')
...@@ -1088,8 +1090,7 @@ class CharmmPsfFile(object): ...@@ -1088,8 +1090,7 @@ class CharmmPsfFile(object):
mm.Discrete2DFunction(num_lj_types, num_lj_types, bcoef)) mm.Discrete2DFunction(num_lj_types, num_lj_types, bcoef))
cforce.addPerParticleParameter('type') cforce.addPerParticleParameter('type')
cforce.setForceGroup(self.NONBONDED_FORCE_GROUP) cforce.setForceGroup(self.NONBONDED_FORCE_GROUP)
if (nonbondedMethod is ff.PME or nonbondedMethod is ff.Ewald or if (nonbondedMethod in (ff.PME, ff.LJPME, ff.Ewald, ff.CutoffPeriodic)):
nonbondedMethod is ff.CutoffPeriodic):
cforce.setNonbondedMethod(cforce.CutoffPeriodic) cforce.setNonbondedMethod(cforce.CutoffPeriodic)
cforce.setCutoffDistance(nonbondedCutoff) cforce.setCutoffDistance(nonbondedCutoff)
cforce.setUseLongRangeCorrection(True) cforce.setUseLongRangeCorrection(True)
......
...@@ -165,11 +165,11 @@ class DesmondDMSFile(object): ...@@ -165,11 +165,11 @@ class DesmondDMSFile(object):
---------- ----------
nonbondedMethod : object=NoCutoff nonbondedMethod : object=NoCutoff
The method to use for nonbonded interactions. Allowed values are The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, or LJPME.
nonbondedCutoff : distance=1*nanometer nonbondedCutoff : distance=1*nanometer
The cutoff distance to use for nonbonded interactions The cutoff distance to use for nonbonded interactions
ewaldErrorTolerance : float=0.0005 ewaldErrorTolerance : float=0.0005
The error tolerance to use if nonbondedMethod is Ewald or PME. The error tolerance to use if nonbondedMethod is Ewald, PME, or LJPME.
removeCMMotion : boolean=True removeCMMotion : boolean=True
If true, a CMMotionRemover will be added to the System If true, a CMMotionRemover will be added to the System
hydrogenMass : mass=None hydrogenMass : mass=None
...@@ -185,7 +185,7 @@ class DesmondDMSFile(object): ...@@ -185,7 +185,7 @@ class DesmondDMSFile(object):
boxSize = self.topology.getUnitCellDimensions() boxSize = self.topology.getUnitCellDimensions()
if boxSize is not None: if boxSize is not None:
sys.setDefaultPeriodicBoxVectors((boxSize[0], 0, 0), (0, boxSize[1], 0), (0, 0, boxSize[2])) sys.setDefaultPeriodicBoxVectors((boxSize[0], 0, 0), (0, boxSize[1], 0), (0, 0, boxSize[2]))
elif nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME): elif nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME, ff.LJPME):
raise ValueError('Illegal nonbonded method for a non-periodic system') raise ValueError('Illegal nonbonded method for a non-periodic system')
# Create all of the particles # Create all of the particles
...@@ -207,7 +207,8 @@ class DesmondDMSFile(object): ...@@ -207,7 +207,8 @@ class DesmondDMSFile(object):
ff.CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic, ff.CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
ff.CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic, ff.CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic,
ff.Ewald:mm.NonbondedForce.Ewald, ff.Ewald:mm.NonbondedForce.Ewald,
ff.PME:mm.NonbondedForce.PME} ff.PME:mm.NonbondedForce.PME,
ff.LJPME:mm.NonbondedForce.LJPME}
nb.setNonbondedMethod(methodMap[nonbondedMethod]) nb.setNonbondedMethod(methodMap[nonbondedMethod])
nb.setCutoffDistance(nonbondedCutoff) nb.setCutoffDistance(nonbondedCutoff)
nb.setEwaldErrorTolerance(ewaldErrorTolerance) nb.setEwaldErrorTolerance(ewaldErrorTolerance)
......
...@@ -115,6 +115,11 @@ class PME(Singleton): ...@@ -115,6 +115,11 @@ class PME(Singleton):
return 'PME' return 'PME'
PME = PME() PME = PME()
class LJPME(Singleton):
def __repr__(self):
return 'LJPME'
LJPME = LJPME()
# Enumerated values for constraint type # Enumerated values for constraint type
class HBonds(Singleton): class HBonds(Singleton):
...@@ -971,7 +976,7 @@ class ForceField(object): ...@@ -971,7 +976,7 @@ class ForceField(object):
The Topology for which to create a System The Topology for which to create a System
nonbondedMethod : object=NoCutoff nonbondedMethod : object=NoCutoff
The method to use for nonbonded interactions. Allowed values are The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, or LJPME.
nonbondedCutoff : distance=1*nanometer nonbondedCutoff : distance=1*nanometer
The cutoff distance to use for nonbonded interactions The cutoff distance to use for nonbonded interactions
constraints : object=None constraints : object=None
...@@ -2195,7 +2200,8 @@ class NonbondedGenerator(object): ...@@ -2195,7 +2200,8 @@ class NonbondedGenerator(object):
CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic, CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic, CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic,
Ewald:mm.NonbondedForce.Ewald, Ewald:mm.NonbondedForce.Ewald,
PME:mm.NonbondedForce.PME} PME:mm.NonbondedForce.PME,
LJPME:mm.NonbondedForce.LJPME}
if nonbondedMethod not in methodMap: if nonbondedMethod not in methodMap:
raise ValueError('Illegal nonbonded method for NonbondedForce') raise ValueError('Illegal nonbonded method for NonbondedForce')
force = mm.NonbondedForce() force = mm.NonbondedForce()
...@@ -2307,7 +2313,7 @@ class LennardJonesGenerator(object): ...@@ -2307,7 +2313,7 @@ class LennardJonesGenerator(object):
self.force.addTabulatedFunction('acoef', mm.Discrete2DFunction(numLjTypes, numLjTypes, acoef)) self.force.addTabulatedFunction('acoef', mm.Discrete2DFunction(numLjTypes, numLjTypes, acoef))
self.force.addTabulatedFunction('bcoef', mm.Discrete2DFunction(numLjTypes, numLjTypes, bcoef)) self.force.addTabulatedFunction('bcoef', mm.Discrete2DFunction(numLjTypes, numLjTypes, bcoef))
self.force.addPerParticleParameter('type') self.force.addPerParticleParameter('type')
if nonbondedMethod in [CutoffPeriodic, Ewald, PME]: if nonbondedMethod in [CutoffPeriodic, Ewald, PME, LJPME]:
self.force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic) self.force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic)
elif nonbondedMethod is NoCutoff: elif nonbondedMethod is NoCutoff:
self.force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff) self.force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff)
......
...@@ -552,7 +552,7 @@ class GromacsTopFile(object): ...@@ -552,7 +552,7 @@ class GromacsTopFile(object):
---------- ----------
nonbondedMethod : object=NoCutoff nonbondedMethod : object=NoCutoff
The method to use for nonbonded interactions. Allowed values are The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, or LJPME.
nonbondedCutoff : distance=1*nanometer nonbondedCutoff : distance=1*nanometer
The cutoff distance to use for nonbonded interactions The cutoff distance to use for nonbonded interactions
constraints : object=None constraints : object=None
...@@ -570,7 +570,7 @@ class GromacsTopFile(object): ...@@ -570,7 +570,7 @@ class GromacsTopFile(object):
The solvent dielectric constant to use in the implicit solvent The solvent dielectric constant to use in the implicit solvent
model. model.
ewaldErrorTolerance : float=0.0005 ewaldErrorTolerance : float=0.0005
The error tolerance to use if nonbondedMethod is Ewald or PME. The error tolerance to use if nonbondedMethod is Ewald, PME, or LJPME.
removeCMMotion : boolean=True removeCMMotion : boolean=True
If true, a CMMotionRemover will be added to the System If true, a CMMotionRemover will be added to the System
hydrogenMass : mass=None hydrogenMass : mass=None
...@@ -589,7 +589,7 @@ class GromacsTopFile(object): ...@@ -589,7 +589,7 @@ class GromacsTopFile(object):
boxVectors = self.topology.getPeriodicBoxVectors() boxVectors = self.topology.getPeriodicBoxVectors()
if boxVectors is not None: if boxVectors is not None:
sys.setDefaultPeriodicBoxVectors(*boxVectors) sys.setDefaultPeriodicBoxVectors(*boxVectors)
elif nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME): elif nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME, ff.LJPME):
raise ValueError('Illegal nonbonded method for a non-periodic system') raise ValueError('Illegal nonbonded method for a non-periodic system')
nb = mm.NonbondedForce() nb = mm.NonbondedForce()
sys.addForce(nb) sys.addForce(nb)
...@@ -877,7 +877,8 @@ class GromacsTopFile(object): ...@@ -877,7 +877,8 @@ class GromacsTopFile(object):
ff.CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic, ff.CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
ff.CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic, ff.CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic,
ff.Ewald:mm.NonbondedForce.Ewald, ff.Ewald:mm.NonbondedForce.Ewald,
ff.PME:mm.NonbondedForce.PME} ff.PME:mm.NonbondedForce.PME,
ff.LJPME:mm.NonbondedForce.LJPME}
nb.setNonbondedMethod(methodMap[nonbondedMethod]) nb.setNonbondedMethod(methodMap[nonbondedMethod])
nb.setCutoffDistance(nonbondedCutoff) nb.setCutoffDistance(nonbondedCutoff)
nb.setEwaldErrorTolerance(ewaldErrorTolerance) nb.setEwaldErrorTolerance(ewaldErrorTolerance)
......
...@@ -776,6 +776,8 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No ...@@ -776,6 +776,8 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
force.setNonbondedMethod(mm.NonbondedForce.Ewald) force.setNonbondedMethod(mm.NonbondedForce.Ewald)
elif nonbondedMethod == 'PME': elif nonbondedMethod == 'PME':
force.setNonbondedMethod(mm.NonbondedForce.PME) force.setNonbondedMethod(mm.NonbondedForce.PME)
elif nonbondedMethod == 'LJPME':
force.setNonbondedMethod(mm.NonbondedForce.LJPME)
else: else:
raise Exception("Cutoff method not understood.") raise Exception("Cutoff method not understood.")
...@@ -885,7 +887,7 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No ...@@ -885,7 +887,7 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
ii, jj, chg, sig, eps = force.getExceptionParameters(i) ii, jj, chg, sig, eps = force.getExceptionParameters(i)
cforce.addExclusion(ii, jj) cforce.addExclusion(ii, jj)
# Now set the various properties based on the NonbondedForce object # Now set the various properties based on the NonbondedForce object
if nonbondedMethod in ('PME', 'Ewald', 'CutoffPeriodic'): if nonbondedMethod in ('PME', 'LJPME', 'Ewald', 'CutoffPeriodic'):
cforce.setNonbondedMethod(cforce.CutoffPeriodic) cforce.setNonbondedMethod(cforce.CutoffPeriodic)
cforce.setCutoffDistance(nonbondedCutoff) cforce.setCutoffDistance(nonbondedCutoff)
cforce.setUseLongRangeCorrection(True) cforce.setUseLongRangeCorrection(True)
......
...@@ -186,6 +186,8 @@ class TestAPIUnits(unittest.TestCase): ...@@ -186,6 +186,8 @@ class TestAPIUnits(unittest.TestCase):
self.assertTrue(force.usesPeriodicBoundaryConditions()) self.assertTrue(force.usesPeriodicBoundaryConditions())
force.setNonbondedMethod(NonbondedForce.PME) force.setNonbondedMethod(NonbondedForce.PME)
self.assertTrue(force.usesPeriodicBoundaryConditions()) self.assertTrue(force.usesPeriodicBoundaryConditions())
force.setNonbondedMethod(NonbondedForce.LJPME)
self.assertTrue(force.usesPeriodicBoundaryConditions())
def testCmapForce(self): def testCmapForce(self):
""" Tests the CMAPTorsionForce API features """ """ Tests the CMAPTorsionForce API features """
......
...@@ -21,12 +21,14 @@ class TestAmberPrmtopFile(unittest.TestCase): ...@@ -21,12 +21,14 @@ class TestAmberPrmtopFile(unittest.TestCase):
"""Test the AmberPrmtopFile.createSystem() method.""" """Test the AmberPrmtopFile.createSystem() method."""
def test_NonbondedMethod(self): def test_NonbondedMethod(self):
"""Test all five options for the nonbondedMethod parameter.""" """Test all six options for the nonbondedMethod parameter."""
methodMap = {NoCutoff:NonbondedForce.NoCutoff, methodMap = {NoCutoff:NonbondedForce.NoCutoff,
CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic, CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic,
CutoffPeriodic:NonbondedForce.CutoffPeriodic, CutoffPeriodic:NonbondedForce.CutoffPeriodic,
Ewald:NonbondedForce.Ewald, PME: NonbondedForce.PME} Ewald:NonbondedForce.Ewald,
PME:NonbondedForce.PME,
LJPME:NonbondedForce.LJPME}
for method in methodMap: for method in methodMap:
system = prmtop1.createSystem(nonbondedMethod=method) system = prmtop1.createSystem(nonbondedMethod=method)
forces = system.getForces() forces = system.getForces()
...@@ -37,7 +39,7 @@ class TestAmberPrmtopFile(unittest.TestCase): ...@@ -37,7 +39,7 @@ class TestAmberPrmtopFile(unittest.TestCase):
def test_Cutoff(self): def test_Cutoff(self):
"""Test to make sure the nonbondedCutoff parameter is passed correctly.""" """Test to make sure the nonbondedCutoff parameter is passed correctly."""
for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME]: for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, LJPME]:
system = prmtop1.createSystem(nonbondedMethod=method, system = prmtop1.createSystem(nonbondedMethod=method,
nonbondedCutoff=2*nanometer, nonbondedCutoff=2*nanometer,
constraints=HBonds) constraints=HBonds)
...@@ -51,7 +53,7 @@ class TestAmberPrmtopFile(unittest.TestCase): ...@@ -51,7 +53,7 @@ class TestAmberPrmtopFile(unittest.TestCase):
def test_EwaldErrorTolerance(self): def test_EwaldErrorTolerance(self):
"""Test to make sure the ewaldErrorTolerance parameter is passed correctly.""" """Test to make sure the ewaldErrorTolerance parameter is passed correctly."""
for method in [Ewald, PME]: for method in [Ewald, PME, LJPME]:
system = prmtop1.createSystem(nonbondedMethod=method, system = prmtop1.createSystem(nonbondedMethod=method,
ewaldErrorTolerance=1e-6, ewaldErrorTolerance=1e-6,
constraints=HBonds) constraints=HBonds)
......
...@@ -14,12 +14,14 @@ class TestDesmondDMSFile(unittest.TestCase): ...@@ -14,12 +14,14 @@ class TestDesmondDMSFile(unittest.TestCase):
self.dms = DesmondDMSFile(path) self.dms = DesmondDMSFile(path)
def test_NonbondedMethod(self): def test_NonbondedMethod(self):
"""Test all five options for the nonbondedMethod parameter.""" """Test all six options for the nonbondedMethod parameter."""
methodMap = {NoCutoff:NonbondedForce.NoCutoff, methodMap = {NoCutoff:NonbondedForce.NoCutoff,
CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic, CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic,
CutoffPeriodic:NonbondedForce.CutoffPeriodic, CutoffPeriodic:NonbondedForce.CutoffPeriodic,
Ewald:NonbondedForce.Ewald, PME: NonbondedForce.PME} Ewald:NonbondedForce.Ewald,
PME:NonbondedForce.PME,
LJPME:NonbondedForce.LJPME}
for method in methodMap: for method in methodMap:
system = self.dms.createSystem(nonbondedMethod=method) system = self.dms.createSystem(nonbondedMethod=method)
forces = system.getForces() forces = system.getForces()
...@@ -30,7 +32,7 @@ class TestDesmondDMSFile(unittest.TestCase): ...@@ -30,7 +32,7 @@ class TestDesmondDMSFile(unittest.TestCase):
def test_Cutoff(self): def test_Cutoff(self):
"""Test to make sure the nonbondedCutoff parameter is passed correctly.""" """Test to make sure the nonbondedCutoff parameter is passed correctly."""
for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME]: for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, LJPME]:
system = self.dms.createSystem(nonbondedMethod=method, system = self.dms.createSystem(nonbondedMethod=method,
nonbondedCutoff=2*nanometer) nonbondedCutoff=2*nanometer)
cutoff_distance = 0.0*nanometer cutoff_distance = 0.0*nanometer
...@@ -43,7 +45,7 @@ class TestDesmondDMSFile(unittest.TestCase): ...@@ -43,7 +45,7 @@ class TestDesmondDMSFile(unittest.TestCase):
def test_EwaldErrorTolerance(self): def test_EwaldErrorTolerance(self):
"""Test to make sure the ewaldErrorTolerance parameter is passed correctly.""" """Test to make sure the ewaldErrorTolerance parameter is passed correctly."""
for method in [Ewald, PME]: for method in [Ewald, PME, LJPME]:
system = self.dms.createSystem(nonbondedMethod=method, system = self.dms.createSystem(nonbondedMethod=method,
ewaldErrorTolerance=1e-6) ewaldErrorTolerance=1e-6)
tolerance = 0 tolerance = 0
......
...@@ -33,12 +33,14 @@ class TestForceField(unittest.TestCase): ...@@ -33,12 +33,14 @@ class TestForceField(unittest.TestCase):
def test_NonbondedMethod(self): def test_NonbondedMethod(self):
"""Test all five options for the nonbondedMethod parameter.""" """Test all six options for the nonbondedMethod parameter."""
methodMap = {NoCutoff:NonbondedForce.NoCutoff, methodMap = {NoCutoff:NonbondedForce.NoCutoff,
CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic, CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic,
CutoffPeriodic:NonbondedForce.CutoffPeriodic, CutoffPeriodic:NonbondedForce.CutoffPeriodic,
Ewald:NonbondedForce.Ewald, PME: NonbondedForce.PME} Ewald:NonbondedForce.Ewald,
PME:NonbondedForce.PME,
LJPME:NonbondedForce.LJPME}
for method in methodMap: for method in methodMap:
system = self.forcefield1.createSystem(self.pdb1.topology, system = self.forcefield1.createSystem(self.pdb1.topology,
nonbondedMethod=method) nonbondedMethod=method)
...@@ -62,7 +64,7 @@ class TestForceField(unittest.TestCase): ...@@ -62,7 +64,7 @@ class TestForceField(unittest.TestCase):
def test_Cutoff(self): def test_Cutoff(self):
"""Test to make sure the nonbondedCutoff parameter is passed correctly.""" """Test to make sure the nonbondedCutoff parameter is passed correctly."""
for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME]: for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, LJPME]:
system = self.forcefield1.createSystem(self.pdb1.topology, system = self.forcefield1.createSystem(self.pdb1.topology,
nonbondedMethod=method, nonbondedMethod=method,
nonbondedCutoff=2*nanometer, nonbondedCutoff=2*nanometer,
...@@ -776,7 +778,7 @@ class AmoebaTestForceField(unittest.TestCase): ...@@ -776,7 +778,7 @@ class AmoebaTestForceField(unittest.TestCase):
def test_NonbondedMethod(self): def test_NonbondedMethod(self):
"""Test all five options for the nonbondedMethod parameter.""" """Test both options for the nonbondedMethod parameter."""
methodMap = {NoCutoff:AmoebaMultipoleForce.NoCutoff, methodMap = {NoCutoff:AmoebaMultipoleForce.NoCutoff,
PME:AmoebaMultipoleForce.PME} PME:AmoebaMultipoleForce.PME}
......
...@@ -22,12 +22,14 @@ class TestGromacsTopFile(unittest.TestCase): ...@@ -22,12 +22,14 @@ class TestGromacsTopFile(unittest.TestCase):
self.top2 = GromacsTopFile('systems/implicit.top') self.top2 = GromacsTopFile('systems/implicit.top')
def test_NonbondedMethod(self): def test_NonbondedMethod(self):
"""Test all five options for the nonbondedMethod parameter.""" """Test all six options for the nonbondedMethod parameter."""
methodMap = {NoCutoff:NonbondedForce.NoCutoff, methodMap = {NoCutoff:NonbondedForce.NoCutoff,
CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic, CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic,
CutoffPeriodic:NonbondedForce.CutoffPeriodic, CutoffPeriodic:NonbondedForce.CutoffPeriodic,
Ewald:NonbondedForce.Ewald, PME: NonbondedForce.PME} Ewald:NonbondedForce.Ewald,
PME:NonbondedForce.PME,
LJPME:NonbondedForce.LJPME}
for method in methodMap: for method in methodMap:
system = self.top1.createSystem(nonbondedMethod=method) system = self.top1.createSystem(nonbondedMethod=method)
forces = system.getForces() forces = system.getForces()
...@@ -52,7 +54,7 @@ class TestGromacsTopFile(unittest.TestCase): ...@@ -52,7 +54,7 @@ class TestGromacsTopFile(unittest.TestCase):
def test_Cutoff(self): def test_Cutoff(self):
"""Test to make sure the nonbondedCutoff parameter is passed correctly.""" """Test to make sure the nonbondedCutoff parameter is passed correctly."""
for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME]: for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME, LJPME]:
system = self.top1.createSystem(nonbondedMethod=method, system = self.top1.createSystem(nonbondedMethod=method,
nonbondedCutoff=2*nanometer, nonbondedCutoff=2*nanometer,
constraints=HBonds) constraints=HBonds)
...@@ -66,7 +68,7 @@ class TestGromacsTopFile(unittest.TestCase): ...@@ -66,7 +68,7 @@ class TestGromacsTopFile(unittest.TestCase):
def test_EwaldErrorTolerance(self): def test_EwaldErrorTolerance(self):
"""Test to make sure the ewaldErrorTolerance parameter is passed correctly.""" """Test to make sure the ewaldErrorTolerance parameter is passed correctly."""
for method in [Ewald, PME]: for method in [Ewald, PME, LJPME]:
system = self.top1.createSystem(nonbondedMethod=method, system = self.top1.createSystem(nonbondedMethod=method,
ewaldErrorTolerance=1e-6, ewaldErrorTolerance=1e-6,
constraints=HBonds) constraints=HBonds)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment