MBG wiki: Editing Performing a molecular dynamics simulation in a truncated octahedral cell

A truncated octahedron is one of the most effective choices for approximately spherical molecules. The trouble is that a truncated octahedron is a somewhat difficult shape to visualise. So, before getting on with the calculus, get a copy of this [[pdb file containing a truncated octahedral system]]. Then do

<code>
mv pdb_file_containing_a_truncated_octahedral_system truncated.pdb
rasmol truncated.pdb
set axes on
set unitcell on
select protein
cpk
</code>

Study the system using the following figure (stolen from Schlick's book) as a guide :

[[image: Truncated octahedral PBC]]

If you are still not sure you can imagine its shape, have a look at 
http://mathworld.wolfram.com/TruncatedOctahedron.html or ask NMG for a copy of his paper-made models.

With so much of an introduction, preparing a system using a truncated octahedral PBC is relatively straightforward. The differences from the standard procedure outlined in previous pages are :

* When using /vmd/'s solvate script you will have to change both the script call *and* the way the solvent box is being build. Instead of building the water box by adding a rim of waters around the solute, you will have to use the <b>-minmax</b> flag to build a (trunctated) cube of waters around the (assumed centered) solute. You will need something like :
<code>
#
# Make water box
#
package require vexpr
package require toctsolvate
toctsolvate psfgen.psf psfgen.pdb -o hydrated -minmax {{-16 -16 -16} {16 16 16}}  -b 1.80 
</code>
In this example, the program will build a truncated octahedron starting from a cube of dimensions 32x32x32 Angstrom^3. It is your responsibility to determine  the size of cube that is sufficiently large to allow for the required separation between neighboring images of the solute.

* In the /NAMD/ scripts you will have to define the unit cell translations for the truncated octahedron (the /cellbasis/ keywords). The matrix you need is

||/d/  || 0   ||   0||
||0  || /d/   ||   0||
||/d//2|| /d//2 || /d//2||

where /d/ is the length of the edge of the cube you've specified in the /vmd/ script. For the example shown above you would have to specify

<code>
cellBasisVector1        32.00    0.00    0.00   
cellBasisVector2         0.00   32.00    0.00   
cellBasisVector3        16.00   16.00   16.00   
cellOrigin               0.00    0.00    0.00   
</code>

Can you work-out why ?

* Do not forget to change the PME grid sizes (taking into account that not all three unit cell dimensions are equal).

-----------

<b>NOTE about pressure :</b> When /vmd/'s solvate script prepares the hydrated system, it throws away all waters containing atoms whose coordinates are outside the specified box. The result is that the density of water molecules near the surface of the unit cell is low, leading to significant negative pressures at the beginning of the simulation. These are fixed (through volume reduction) when the barostat is turned on. You can probably avoid this by specifying in the /vmd/ script a somewhat larger cubic cell dimensions that those intended (for example, you can ask /vmd/ to prepare a cell that is 33x33x33 A, and then tell /NAMD/ that the cell edge is 32A).

Summary:

This change is a minor edit.

Username:

Replace this text with a file.