MBG wiki: 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

mv pdb_file_containing_a_truncated_octahedral_system truncated.pdb
rasmol truncated.pdb
set axes on
set unitcell on
select protein
cpk

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

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 -minmax flag to build a (trunctated) cube of waters around the (assumed centered) solute. You will need something like :

#
# 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

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

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

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

NOTE about pressure : 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).