"vscode:/vscode.git/clone" did not exist on "a16af40531f6c5dabfd7f0b22b835ecfe3e43464"
Commit f16a609b authored by Jason Swails's avatar Jason Swails
Browse files

Change implicitSolentKappa to implicitSolventSaltConc and convert to

implicitSolventKappa the same (black-magic) way that Amber does. I have never
been able to reproduce the prefactor used to convert to kappa in Amber... sigh.
parent 901108a0
...@@ -31,6 +31,7 @@ USE OR OTHER DEALINGS IN THE SOFTWARE. ...@@ -31,6 +31,7 @@ USE OR OTHER DEALINGS IN THE SOFTWARE.
__author__ = "Peter Eastman" __author__ = "Peter Eastman"
__version__ = "1.0" __version__ = "1.0"
from math import sqrt
from simtk.openmm.app import Topology from simtk.openmm.app import Topology
from simtk.openmm.app import PDBFile from simtk.openmm.app import PDBFile
from simtk.openmm.app.internal import amber_file_parser from simtk.openmm.app.internal import amber_file_parser
...@@ -142,7 +143,7 @@ class AmberPrmtopFile(object): ...@@ -142,7 +143,7 @@ class AmberPrmtopFile(object):
def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer, def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
constraints=None, rigidWater=True, implicitSolvent=None, constraints=None, rigidWater=True, implicitSolvent=None,
implicitSolventKappa=0.0*(1/unit.nanometer), implicitSolventSaltConc=0.0*(unit.moles/unit.liter),
soluteDielectric=1.0, solventDielectric=78.5, soluteDielectric=1.0, solventDielectric=78.5,
removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005): removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005):
"""Construct an OpenMM System representing the topology described by this prmtop file. """Construct an OpenMM System representing the topology described by this prmtop file.
...@@ -155,7 +156,8 @@ class AmberPrmtopFile(object): ...@@ -155,7 +156,8 @@ class AmberPrmtopFile(object):
Allowed values are None, HBonds, AllBonds, or HAngles. Allowed values are None, HBonds, AllBonds, or HAngles.
- rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument - rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument
- implicitSolvent (object=None) If not None, the implicit solvent model to use. Allowed values are HCT, OBC1, OBC2, GBn, or GBn2. - implicitSolvent (object=None) If not None, the implicit solvent model to use. Allowed values are HCT, OBC1, OBC2, GBn, or GBn2.
- implicitSolventKappa (float=0.0*1/unit.nanometer) The Debye-screening parameter corresponding to ionic strength used for implicit solvent - implicitSolventSaltConc (float=0.0*unit.moles/unit.liter) The salt concentration for GB
calculations (modelled as a debye screening parameter)
- soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model. - soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model.
- solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model. - solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model.
- removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System - removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
...@@ -197,6 +199,16 @@ class AmberPrmtopFile(object): ...@@ -197,6 +199,16 @@ class AmberPrmtopFile(object):
implicitString = 'GBn2' implicitString = 'GBn2'
else: else:
raise ValueError('Illegal value for implicit solvent model') raise ValueError('Illegal value for implicit solvent model')
# Use the conversion of ionic strength to kappa used by Amber (external
# dielectric of 78.5 -- it's unclear where the 0.10806 comes from
# though, although code comments in Amber say 'assuming M salt
# concentration, T=298.15, and epsext=78.5')
if unit.is_quantity(implicitSolventSaltConc):
implicitSolventSaltConc = implicitSolventSaltConc.value_in_unit(unit.moles/unit.liter)
# 0.73 to account for ion exclusions. The Amber Debye screening length
# units will be in 1/angstrom given these units. Not a pretty way to do
# this...
implicitSolventKappa = 0.73 * sqrt(0.10806 * implicitSolventSaltConc) / unit.angstrom
sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString, nonbondedCutoff=nonbondedCutoff, sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString, nonbondedCutoff=nonbondedCutoff,
nonbondedMethod=methodMap[nonbondedMethod], flexibleConstraints=False, gbmodel=implicitString, nonbondedMethod=methodMap[nonbondedMethod], flexibleConstraints=False, gbmodel=implicitString,
soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, implicitSolventKappa=implicitSolventKappa, soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, implicitSolventKappa=implicitSolventKappa,
......
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