# Step 8

The first step in solvating a system is to determine the size of the solvent box and the number of counterions, if any, that are to be added. This is achieved with the following simple program:

``# . Get the system.``
``system = Unpickle ( "../step7/step7.pkl" )``
``system.Summary ( )``

``# . Determine solvation parameters.``
``DetermineSolvationParameters ( system ) ``

The function `DetermineSolvationParameters` prints out information about possible (orthorhombic) solvation boxes for the system, including the numbers of counterions that it is necessary to add to have a solution of a particular molarity. It is usual to have a solvent layer of sufficient thickness such that proteins within adjacent boxes interact as little as possible. For PKA, we choose to have a solvent buffer of 10 Å minimum between the protein and the edge of the box (i.e. at least 20 Å between protein molecules) which means having a box with approximate dimensions 75 × 60 × 60 Å3.

The total charge of the PKA system is +3 so 3 monovalent anions would need to be added to have a simulation system that is neutral overall. However, we also would like to mimic the molarity of the intracellular environment within which the protein normally works. If this is taken to be 0.154 M, then approximately 24 monovalent cations and 27 monovalent anions are required for a box of the specified volume.

A solvent box of the appropriate size can be constructed using a program very similar to that of Example 27 in the pDynamo distribution. The program is straightforward and is:

``# . Parameters.``
``# . Density of water (kg m^-3).``
``_DENSITY = 1000.0``
``# . Box sizes.``
``_XBOX = 75.0``
``_YBOX = 60.0``
``_ZBOX = 60.0``
``# . Monte Carlo modifiers.``
``_MULTIPLIER1 =  1``
``_MULTIPLIER2 = 15``

``# . Define the solvent MM and NB models.``
``mmModel = MMModelOPLS ( "bookSmallExamples" )``
``nbModel = NBModelMonteCarlo ( )``

``# . Define the solvent molecule.``
``molecule = MOLFile_ToSystem ( "../data/water.mol" ) )``
``molecule.Summary ( )``

``# . Create a symmetryparameters instance with the correct dimensions.``
``symmetryParameters = SymmetryParameters ( )``
``symmetryParameters.SetCrystalParameters ( _XBOX, _YBOX, _ZBOX, 90.0, 90.0, 90.0 )``

``# . Create the basic solvent box.``
``solvent       = BuildSolventBox ( CrystalClassOrthorhombic ( ), symmetryParameters, molecule, _DENSITY )``
``solvent.label = "Water Box"``
``solvent.DefineMMModel ( mmModel )``
``solvent.DefineNBModel ( nbModel )``
``solvent.Summary ( )``
``solvent.Energy  ( )``

``# . Do Monte Carlo simulations to equilibrate the system.``
``nmolecules = len ( solvent.atoms ) / len ( molecule.atoms )``
``MonteCarlo_SystemGeometry ( solvent,                                     \``
``                            blocks          =                       100, \``
``                            log             =                      None, \``
``                            moves           = _MULTIPLIER1 * nmolecules, \``
``                            volumeFrequency =                         0  )``
``solvent.Energy  ( )``
``MonteCarlo_SystemGeometry ( solvent,                                     \``
``                            blocks          =                       100, \``
``                            moves           = _MULTIPLIER2 * nmolecules, \``
``                            volumeFrequency =                         0  )``
``solvent.Energy  ( ) ``

``# . Calculate and print the final density.``
``logFile.Paragraph ( "Solvent density = %.2f kg m^-3." % ( SystemDensity ( solvent ), ) )``

``# . Save the water box.``
``solvent.configuration.Clear ( )``
``Pickle ( "waterBox.pkl", solvent ) ``

The function `BuildSolventBox` works by placing solvent molecules evenly, but in random orientations, throughout a box of appropriate shape and dimension. This is followed up by two constant-volume Monte Carlo calculations, the first to obtain a configuration of reasonable energy (with few or no intermolecular overlaps) and the second to more fully equilibrate the system.

The third stage in the solvation process is to add counterions to the protein system. The full program is:

``# . Parameters.``
``# . Box sizes.``
``_XBOX = 75.0``
``_YBOX = 60.0``
``_ZBOX = 60.0``
``# . Number and type of ions to add.``
``_NNEGATIVE   = 27``
``_NPOSITIVE   = 24``
``_NEGATIVEION = "chloride"``
``_POSITIVEION = "potassium"``
``# . Reorient option.``
``QREORIENT = False``

``# . Define the solvent MM model.``
``mmModel = MMModelOPLS ( "bookSmallExamples" )``

``# . Get the system.``
``system = Unpickle ( "../step7/step7.pkl" )``
``system.Summary ( )``

``# . Reorient the system if necessary (see the results of GetSolvationInformation.py).``
``if QREORIENT:``
``    masses = system.atoms.GetItemAttributes ( "mass" )``
``    system.coordinates3.ToPrincipalAxes ( weights = masses )``

``# . Get the positive and negative ions.``
``if _NNEGATIVE > 0:``
``    anion = MOLFile_ToSystem ( "../data/" + _NEGATIVEION + ".mol" )``
``    anion.DefineMMModel ( mmModel )``
``    anion.Summary ( )``
``if _NPOSITIVE > 0:``
``    cation = MOLFile_ToSystem ( "../data/" + _POSITIVEION + ".mol" )``
``    cation.DefineMMModel ( mmModel )``
``    cation.Summary ( )``

``# . Add the counterions.``
``newSystem = AddCounterIons ( system, _NNEGATIVE, anion, _NPOSITIVE, cation, ( _XBOX, _YBOX, _ZBOX ) )``

``# . Save the combined system.``
``newSystem.configuration.Clear ( )``
``Pickle ( "step8_a.pkl", newSystem ) ``

After some initial parameter definitions, the program starts by reading in the protein system and then, optionally, reorienting the protein coordinates. This is not done here as the results of the program `GetSolvationInformation` suggest that this is not necessary to obtain a simulation box of reasonable size. After this, systems corresponding to a single anion (chloride) and a single cation (potassium) are defined and the function `AddCounterIons` is invoked which places the required number of anions and cations within a certain box around the protein. Finally, the new, merged protein/ion system is saved for subsequent use.

The function `AddCounterIons` places ions at random positions around the protein. It is, in principle, possible to use algorithms which use an energetic criterion for placement (i.e. choose positions with the lowest energy) but similar results can be obtained, if necessary, by carrying out Langevin molecular dynamics or Monte Carlo simulations on the combined protein/ion system. This would require using an implicit solvent model (e.g. an NB model with a dielectric appropriate to that of water – about 80), fixing the protein atoms, performing a simulation and then selecting energetically favorable ion configurations.

The final step in the solvation procedure is to merge the protein/ion system and the water box. This involves a procedure that is very similar to that of Example 28 in the pDynamo distribution. The program is:

``# . Define the NB model.``
``nbModel = NBModelABFS ( )``

``# . Get the solute system.``
``solute = Unpickle ( "step8_a.pkl" )``
``solute.Summary ( )``

``# . Get the solvent system.``
``solvent = Unpickle ( "waterBox.pkl" )``
``solvent.DefineNBModel ( nbModel )``
``solvent.Summary ( )``

``# . Create the solvated system.``
``solution       = SolvateSystemBySuperposition ( solute, solvent, orientSolute = False )``
``solution.label = "Solvated PKA Simulation System - Chains A and I Only"``
``solution.DefineNBModel ( nbModel )``
``solution.Summary ( )``

``# . Fix solute atoms (A and I chains only).``
``fixedAtoms = Selection ( AtomSelection.FromAtomPattern ( solution, "A:*:*" ) | AtomSelection.FromAtomPattern ( solution, "I:*:*" ) )``
``solution.DefineFixedAtoms ( fixedAtoms )``
``solution.Summary ( )``

``# . Optimization.``
``ConjugateGradientMinimize_SystemGeometry ( solution,                    \``
``                                           maximumIterations    =  200, \``
``                                           logFrequency         =    1, \``
``                                           rmsGradientTolerance = 10.0  )``

``# . Define a random number generator in a given state.``
``rng = random.Random ( )``
``rng.seed ( 517151 )``

``# . Do a Langevin dynamics calculation to refine the solvent and counterion coordinates.``
``LangevinDynamics_SystemGeometry ( solution,                   \``
``                                  collisionFrequency =  25.0, \``
``                                  logFrequency       =   100, \``
``                                  rng                =   rng, \``
``                                  steps              = 10000, \``
``                                  temperature        = 300.0, \``
``                                  timeStep           = 0.001  )``

``# . Save the system.``
``solution.configuration.Clear ( )``
``Pickle ( "step8_b.pkl", solution ) ``

As in Example 28, the program employs the function `SolvateSystemBySuperposition` to perform the solvation. This is followed by a refinement of the water and ion coordinates, keeping the protein atoms fixed, using both geometry optimization and molecular dynamics simulation.

The strategy used by `SolvateSystemBySuperposition` is to overlap the solute and solvent systems and then remove those solvent molecules which overlap with solute atoms. This is a straightforward method to implement but the solvated systems it produces depend crucially upon the initial configurations of the solute and solvent systems. This means, for example, that certain regions of the solute, such as internal cavities or nooks and crannies on its surface, may be insufficiently solvated. This problem can be alleviated by performing several cycles of solvation (by superposition) and, optionally, refinement using different solvent configurations. There also exist more specialized solvation algorithms, such as those that search for solvatable regions within or on the surface of a solute. These may be implement in future versions of pDynamo.