The final step in the set-up process is to equilibrate the solvated system. There is no unique way to do this but a reasonable procedure follows that of the refinement of step 7. This involves performing a series of molecular dynamics simulations in which the heavy atoms of the protein are harmonically tethered. The tether force constants are large initially but are gradually decreased to zero. The program is:
# . Parameters.
NSTEPS = 1000
# . Get the system with no fixed atoms and an appropriate NB model.
system = Unpickle ( "../step8/step8_b.pkl" )
system.DefineFixedAtoms ( None )
system.DefineNBModel ( NBModelABFS ( ) )
system.Summary ( )
# . Get a selection for the heavy atoms in the A and I chains.
# . Protein atoms - A and I chains.
protein = ( AtomSelection.FromAtomPattern ( system, "A:*:*" ) | AtomSelection.FromAtomPattern ( system, "I:*:*" ) )
hydrogens = AtomSelection.FromAtomAttributes ( system, "atomicNumber", 1 )
heavies = Selection ( protein - hydrogens )
# . Define a random number generator in a given state.
rng = random.Random ( )
rng.seed ( 511717 )
# . Loop over dynamics.
for forceConstant in ( 1000.0, 500.0, 250.0, 100.0, 50.0, 10.0, 4.0, 1.0, 0.25, None ):
# . Harmonically constrain protein heavy atoms.
tethers = None
if ( forceConstant is not None ):
tetherEnergyModel = SoftConstraintEnergyModelHarmonic ( 0.0, forceConstant )
tethers = SoftConstraintContainer ( )
tethers["tethers"] = SoftConstraintMultipleTether ( heavies, reference, tetherEnergyModel )
system.DefineSoftConstraints ( tethers )
# . Dynamics.
LangevinDynamics_SystemGeometry ( system, \
collisionFrequency = 25.0, \
logFrequency = 10, \
rng = rng, \
steps = NSTEPS, \
temperature = 300.0, \
timeStep = 0.001 )
# . Save the system.
system.configuration.Clear ( )
Pickle ( "step9.pkl", system )
For the purposes of illustration, 10 ps of Langevin (NVT) dynamics are performed here. In actual applications, however, equilibration will normally need much longer simulations, especially when no tether constraints are being applied, and it may also be preferable to use alternative sampling algorithms, such as NPT dynamics or even Monte Carlo. Likewise, some studies will require modifications to the above procedure. Thus, for example, a hybrid potential investigation of an enzyme reaction would probably involve a partial pre-equilibration using a MM potential before finishing off with simulations with a QC/MM energy model.