Input file Preparation -- Adding Water and Ions with tleap

We will now use tleap to add a water box to the molecule and counter ions to complete the model system.

> solvateBox mol WATBOX216 10
Solute vdw bounding box: 49.995 53.684 37.063
Total bounding box for atom centers: 69.995 73.684 57.063
Solvent unit box: 18.774 18.774 18.774
Total vdw box size: 73.439 76.851 60.149 angstroms.
Volume: 339472.593 A^3
Total mass 162173.797 amu, Density 0.793 g/cc
Added 8064 residues.

As shown above, the command to add water box is called "solvateBox".  There are several other ways to add water, but we'll just stick to the most straightforward way here.  WATERBOX216 is a pre-equilibrated box of TIP3P water.  The number 10 is the buffering distance between the edges of the box and the protein in angastroms.  Here you have to make some judgement calls about the size of buffer to use.  If you use too big a number, you end up with a big water box and waste a lot of unnecessary computing time on uninteresting water molecules.  However, if you use too small a water box, during the simulation, the molecule may undergo conformation changes and part of it may stick outside of the box.  If you want to simulate as close to experimental conditions and you know the concentration of your system, you can use that information to figure out the size of water box you need and set it explicitly, otherwise, I think 10 is a reasonable number to start.

The next thing we need to do is to add counter ions.  Before we issue the "addions" command, we need to figure out whether our system is positively charged or negatively charge.  If it is positively charge, we will want to add negatively charged Cl- ions to counter it and if it is negatively charged, then we'll add Na+ to counter it.  To calculate the charge of our system, we can use the command "charge" like this:

> charge mol
Total unperturbed charge: 2.00
Total perturbed charge: 2.00

So we see that our system is positively charged.  We'll add Cl- ions to balance out the charges. 

AMBER actually offers two algorithms to add ions.  The first approach implemented in "addions" is to simply draw a grid around the solute and places ions at grid points where the energies are lowest.  This approach will ignore water molecules in locating where to place the ion and if the chosen location overlaps a water molecule, the water is deleted and replaced with the ion.  If we use this algorithm, we ended up with Cl- ions next to our Mg2+ which is not what we want.  The second approach implemented by the command "addions2" does almost the same thing as "addions" except that it treats solvent molecules the same as solute.  We'll use "addions2" to ensure that the Cl- ions are some distance away from our molecule so its charges won't artificially distort our system.  The number "0" at the end of our command means that we want tleap to figure out the appropriate number of counter ions to neutralize the whole system.

As you can see, this calculation took 990 seconds on my Pentium IV PC, so if nothing seems to happen on your computer for a while, just exercise a little bit of patience before you reach for the reset button. 

> addIons2 mol Cl- 0
2 Cl- ions required to neutralize.
Adding 2 counter ions to "mol" using 1A grid
Grid extends from solute vdw + 2.47 to 8.47
Resolution: 1.00 Angstrom.
grid build: 4 sec
Calculating grid charges
charges: 990 sec
Placed Cl- in mol at (-18.42, 1.87, 29.98).
Placed Cl- in mol at (0.58, -17.13, 29.98).

Done adding ions.

Finally, we are ready to output our results and save the system as input files.

> saveAmberParm mol wt1mg.parm7 wt1mg.crd
Checking Unit.
Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
total 468 improper torsions applied
Building H-Bond parameters.
Marking per-residue atom chain types.
(Residues lacking connect0/connect1 -
these don't have chain types marked:

res total affected

WAT 8064
(no restraints)

You can now quit tleap.