/usr/share/pyshared/MMTK/Solvation.py is in python-mmtk 2.7.9-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 | # This module contains code for solvation.
#
# Written by Konrad Hinsen
#
"""
Solvation of solute molecules
"""
__docformat__ = 'restructuredtext'
from MMTK import ChemicalObjects, Units, Universe
from MMTK.MolecularSurface import surfaceAndVolume
from MMTK.Minimization import SteepestDescentMinimizer
from MMTK.Dynamics import VelocityVerletIntegrator, VelocityScaler, \
TranslationRemover, RotationRemover
from MMTK.Trajectory import Trajectory, TrajectoryOutput, SnapshotGenerator, \
StandardLogOutput
import copy
#
# Calculate the number of solvent molecules
#
def numberOfSolventMolecules(universe, solvent, density):
"""
:param universe: a finite universe
:type universe: :class:`~MMTK.Universe.Universe`
:param solvent: a molecule, or the name of a molecule in the database
:type solvent: :class:`~MMTK.ChemicalObjects.Molecule` or str
:param density: the density of the solvent (amu/nm**3)
:type density: float
:returns: the number of solvent molecules that must be added to the
universe, in addition to whatever it already contains,
to obtain the given solvent density.
:rtype: int
"""
if isinstance(solvent, basestring):
solvent = ChemicalObjects.Molecule(solvent)
cell_volume = universe.cellVolume()
if cell_volume is None:
raise TypeError("universe volume is undefined")
solute_volume = 0.
for o in universe._objects:
solute_volume = solute_volume + surfaceAndVolume(o)[1]
return int(round(density*(cell_volume-solute_volume)/solvent.mass()))
#
# Add solvent to a universe containing solute molecules.
#
def addSolvent(universe, solvent, density, scale_factor=4.):
"""
Scales up the universe and adds as many solvent molecules
as are necessary to obtain the specified solvent density,
taking account of the solute molecules that are already present
in the universe. The molecules are placed at random positions
in the scaled-up universe, but without overlaps between
any two molecules.
:param universe: a finite universe
:type universe: :class:`~MMTK.Universe.Universe`
:param solvent: a molecule, or the name of a molecule in the database
:type solvent: :class:`~MMTK.ChemicalObjects.Molecule` or str
:param density: the density of the solvent (amu/nm**3)
:type density: float
:param scale_factor: the factor by which the initial universe is
expanded before adding the solvent molecules
:type scale_factor: float
"""
# Calculate number of solvent molecules and universe size
if isinstance(solvent, basestring):
solvent = ChemicalObjects.Molecule(solvent)
cell_volume = universe.cellVolume()
if cell_volume is None:
raise TypeError("universe volume is undefined")
solute = copy.copy(universe._objects)
solute_volume = 0.
excluded_regions = []
for o in solute:
solute_volume = solute_volume + surfaceAndVolume(o)[1]
excluded_regions.append(o.boundingSphere())
n_solvent = int(round(density*(cell_volume-solute_volume)/solvent.mass()))
solvent_volume = n_solvent*solvent.mass()/density
cell_volume = solvent_volume + solute_volume
universe.translateBy(-solute.position())
universe.scaleSize((cell_volume/universe.cellVolume())**(1./3.))
# Scale up the universe and add solvent molecules at random positions
universe.scaleSize(scale_factor)
universe.scale_factor = scale_factor
for i in range(n_solvent):
m = copy.copy(solvent)
m.translateTo(universe.randomPoint())
while True:
s = m.boundingSphere()
collision = False
for region in excluded_regions:
if (s.center-region.center).length() < s.radius+region.radius:
collision = True
break
if not collision:
break
m.translateTo(universe.randomPoint())
universe.addObject(m)
excluded_regions.append(s)
#
# Shrink the universe to its final size
#
def shrinkUniverse(universe, temperature=300.*Units.K, trajectory=None,
scale_factor=0.95):
"""
Shrinks the universe, which must have been scaled up by
:class:`~MMTK.Solvation.addSolvent`, back to its original size.
The compression is performed in small steps, in between which
some energy minimization and molecular dynamics steps are executed.
The molecular dynamics is run at the given temperature, and
an optional trajectory can be specified in which intermediate
configurations are stored.
:param universe: a finite universe
:type universe: :class:`~MMTK.Universe.Universe`
:param temperature: the temperature at which the Molecular Dynamics
steps are run
:type temperature: float
:param trajectory: the trajectory in which the progress of the
shrinking procedure is stored, or a filename
:type trajectory: :class:`~MMTK.Trajectory.Trajectory` or str
:param scale_factor: the factor by which the universe is scaled
at each reduction step
:type scale_factor: float
"""
# Set velocities and initialize trajectory output
universe.initializeVelocitiesToTemperature(temperature)
if trajectory is not None:
if isinstance(trajectory, basestring):
trajectory = Trajectory(universe, trajectory, "w",
"solvation protocol")
close_trajectory = True
else:
close_trajectory = False
actions = [TrajectoryOutput(trajectory, ["configuration"], 0, None, 1)]
snapshot = SnapshotGenerator(universe, actions=actions)
snapshot()
# Do some minimization and equilibration
minimizer = SteepestDescentMinimizer(universe, step_size = 0.05*Units.Ang)
actions = [VelocityScaler(temperature, 0.01*temperature, 0, None, 1),
TranslationRemover(0, None, 20)]
integrator = VelocityVerletIntegrator(universe, delta_t = 0.5*Units.fs,
actions = actions)
for i in range(5):
minimizer(steps = 40)
integrator(steps = 200)
# Scale down the system in small steps
i = 0
while universe.scale_factor > 1.:
if trajectory is not None and i % 1 == 0:
snapshot()
i = i + 1
step_factor = max(scale_factor, 1./universe.scale_factor)
for object in universe:
object.translateTo(step_factor*object.position())
universe.scaleSize(step_factor)
universe.scale_factor = universe.scale_factor*step_factor
for i in range(3):
minimizer(steps = 10)
integrator(steps = 50)
del universe.scale_factor
if trajectory is not None:
snapshot()
if close_trajectory:
trajectory.close()
|