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.pkl" )

system.DefineNBModel ( NBModelCutOff.WithDefaults ( ) )

system.Summary ( )

# . Get a selection with the indices of the atoms whose coordinates were built.

hydrogens = [ index for ( index, atom ) in enumerate ( system.atoms ) if atom.atomicNumber == 1 ]

built = Selection.FromIterable ( hydrogens + [ system.sequence.AtomIndex ( "A:LYS.8:CB" ), system.sequence.AtomIndex ( "I:SER.17:OG" ) ] )

# . Minimize the coordinates of these atoms.

system.freeAtoms = built

ConjugateGradientMinimize_SystemGeometry ( system ,

maximumIterations = _Iterations ,

logFrequency = _Logfrequency ,

rmsGradientTolerance = _Tolerance )

system.freeAtoms = None

# . Get a selection corresponding to heavy atoms.

heavies = Selection.FromIterable ( [ index for ( index, atom ) in enumerate ( system.atoms ) if atom.atomicNumber > 1 ] )

# . 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 restrain heavy atoms.

tethers = None

if ( forceConstant is not None ):

reference = Clone ( system.coordinates3 )

tetherEnergyModel = RestraintEnergyModel.Harmonic ( 0.0, forceConstant )

tethers = RestraintModel ( )

tethers["Tethers"] = RestraintMultipleTether.WithOptions ( energyModel = tetherEnergyModel ,

reference = reference ,

selection = heavies )

system.DefineRestraintModel ( tethers )

# . Minimize.

ConjugateGradientMinimize_SystemGeometry ( system ,

maximumIterations = _Iterations ,

logFrequency = _Logfrequency ,

rmsGradientTolerance = _Tolerance )

system.Energy ( doGradients = True )

# . Save the system.

Pickle ( "step7.pkl" )

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 assigning them to the freeAtoms attribute of system. The optimization is then performed and afterwards all atoms are freed by re-setting the freeAtoms attribute of system to 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.