Tutorials‎ > ‎PDB Files‎ > ‎

Step 7

The aim of this step is to refine or to relax the structure of the current (vacuum) simulation system. This need not be very precise, because the system is to be solvated in a subsequent step, but should be sufficient to relieve unfavourable non-bonding interactions and other high-energy structural features.

There is no unique prescription for performing the refinement but it can involve a mixture of geometry optimizations and molecular dynamics simulations and should be done in a controlled fashion so that there are no abrupt changes in the system's structure. In the case of PKA, we decide to use only geometry optimization, but molecular dynamics simulations are easily included in the procedure if this is deemed necessary. The program is:

# . Parameters.
_Iterations = 1000
_Logfrequency = 100
_Tolerance = 1.0

# . Get the system with an appropriate NB model.
system = Unpickle ( "../step6/step6.pkl" )
system.DefineNBModel ( NBModelABFS ( ) )
system.Summary ( )

# . Get a selection with the indices of the atoms whose coordinates were built.
built = AtomSelection.FromAtomAttributes ( system, "atomicNumber", 1 )
built.Add ( system.sequence.AtomIndex ( "A:LYS.8:CB" ) )
built.Add ( system.sequence.AtomIndex ( "I:SER.17:OG" ) )
built = Selection ( built )

# . Minimize the coordinates of these atoms.
system.DefineFixedAtoms ( built.Complement ( upperBound = len ( system.atoms ) ) )
ConjugateGradientMinimize_SystemGeometry ( system, \
maximumIterations = _Iterations, \
logFrequency = _Logfrequency, \
rmsGradientTolerance = _Tolerance )
system.DefineFixedAtoms ( None )

# . Get a selection corresponding to heavy atoms.
hydrogens =
Selection ( AtomSelection.FromAtomAttributes ( system, "atomicNumber", 1 ) )
heavies = hydrogens.Complement ( upperBound = len ( system.atoms ) )

# . Loop over minimizations.
for forceConstant in ( 1000.0, 500.0, 250.0, 100.0, 50.0, 10.0, 4.0, 1.0, 0.25, None ):

# . Harmonically constrain heavy atoms.
tethers = None
if ( forceConstant is not None ):
reference = Clone ( system.coordinates3 )
tetherEnergyModel = SoftConstraintEnergyModelHarmonic ( 0.0, forceConstant )
tethers = SoftConstraintContainer ( )
tethers["tethers"] = SoftConstraintMultipleTether ( heavies, reference, tetherEnergyModel )
system.DefineSoftConstraints ( tethers )

# . Minimize.
ConjugateGradientMinimize_SystemGeometry ( system, \
maximumIterations = _Iterations, \
logFrequency = _Logfrequency, \
rmsGradientTolerance = _Tolerance )
system.Energy ( doGradients = True )

# . Save the system.
system.configuration.Clear ( )
Pickle ( "step7.pkl", system )

The program has two sections. In the first part, only the coordinates of the atoms that were built in step 6 are optimized. This is achieved by generating a selection, built, that contains the indices of the built atoms, and passing the complement of this selection, which corresponds to the indices of the non-built atoms, to the method DefineFixedAtoms of the class System. The optimization is then performed and afterwards all atoms are freed by invoking the DefineFixedAtoms method with the argument None.

In the second section, the coordinates of all atoms are optimized simultaneously. This is done by carrying out a sequence of optimizations in which the non-hydrogen atoms are harmonically-tethered about the positions that they had at the start of the current optimization. The force constant used for the tethering starts off with a high value but is gradually reduced until, in the final iteration, no tethering is applied and the optimization is unconstrained.

As a final point, it is worth remarking that the RMS gradient tolerance termination criterion for the optimizations is set to the reasonably high value of 1.0. This is more than adequate for the level of refinement that is required here.