Tutorial - Part 2

This section discusses the methaneaatodata.tcl script that builds an initial topology for a bulk methane simulation from the methane.pdb coordinate file. The process of building the topology mostly follows the strategy in the previous scripts. The first difference is in loading the coordinates. The atom names in the coordinate file allow VMD to correctly identify the elements and are given in angstrom, so we can trust the automatically computed bonds and don't have to adjust atomic radii:mol new methane.pdb autobonds yes waitfor all

Next we create three atom selections, one for the hydrogens, one for the carbon and one for all atoms and set parameters relevant to the force field:set selh [atomselect top {name H}]

$selh set type HC

$selh set mass 1.00800

$selh set charge 0.060

set selc [atomselect top {name C}]

$selc set type CT

$selc set mass 12.01100

$selc set charge -0.240 ; # = 4x 0.060

set sel [atomselect top all]

The carbon is a "tetrahedral" carbon and thus we set its type to CT. Hydrogens at tetrahedral aliphatic carbons have a charge of +0.060 and thus the carbon with four hydrogens would have a charge of -0.24. We still need to retype the bonds and guess the angles. OPLS-aa also defines dihedral interactions, however, the topo guessdihedrals command will (expectedly) come up empty for methane. The next step is to place the methane molecule into the center of the box and determine suitable box dimensions. The VMD Tcl script code to do this is (note the use of -withradii for measure minmax, which makes sure there is sufficient space between methane molecules; also note that we enlarge the box around the molecule by an additional, empirical 10% using the "vecscale" command to account for thermal expansion):

set center [measure center $sel weight none]

$sel moveby [vecscale -1.0 $center]

set minmax [measure minmax $sel -withradii]

set box [vecscale 1.1 [vecsub [lindex $minmax 1] [lindex $minmax 0]]]

pbc set $box

Now we need to create a 4x4x4 lattice of the single methane molecule to construct a bulk system with 64 molecules. This could be done with the replicate command in LAMMPS, but TopoTools has a special high-level utility method for it as well.

TopoTools::replicatemol top 4 4 4

This will create a new molecule containing the requested replications of the original single molecule data set. Now we can write out the result to a LAMMPS data file, and also write out the coordinate and topology information as .pdb/.psf file combo, so we can use the more space efficient and faster to read .dcd output in the MD code.

topo writelammpsdata data.step2a full

animate write pdb 64xmethane.pdb

animate write psf 64xmethane.psf

This input still needs to be equilibrated to be run at the desired temperature and density or pressure conditions. The attached input in.step2a sets up some initial equilibration at a temperature of 100K with constant volume. Methane should be liquid at this temperature and if the guess for the molecular volume via measure minmax was good, then the pressure should be fluctuating around 0.0 atm. The second input in.step2a-2 replaces 'fix nve' plus 'fix langevin' with 'fix nvt'. This removes the damping and randomization of velocities of the langevin thermostat, which helps in achieving equipartitioning of the kinetic energy quickly, with a proper NVT ensemble time integration. This results in a more realistic dynamics of the atoms and would need to be continued until the system is properly equilibrated.

Step 2b: Ethane

Setting up an initial configuration for ethane from the file ethane.pdb is very similar to the previous script. The commands in the script ethaneaatodata.tcl remain mostly the same, since we still have only two atom types (both carbons are "terminal"), however, the charge on the carbons is different (-0.18 = -1.0*3*0.06), since they are now bound to only three hydrogens instead of four.set selc [atomselect top {name C}]

$selc set type CT

$selc set mass 12.01100

$selc set charge -0.180 ; # = 3x 0.060