Tutorial - Part 1

This paragraph discusses the simplexyztodata.tcl script. The first block of script code loads the TopoTools and pbctools plugins and makes sure that they are of a sufficiently high version. This test makes it look a bit more complicated, but effectively this is the same as these two lines:package require topotools

package require pbctools

The next step is to load the molecular coordinate data. We use the autobonds no flag to inhibit any automatic bond searches:mol new chain.xyz autobonds no waitfor all

Now we need to assign data to the atoms that is required for writing a proper LAMMPS data file, but not available in the .xyz file. This is done by creating an atom selection and then assigning properties to all atoms in that selection:

set sel [atomselect top all]

$sel set radius 0.85

$sel set name A

$sel set type A

$sel set mass 1.0

This specific system is to be used with reduced units (units lj in LAMMPS input) and thus we need to set a reasonable radius and mass value. This now is sufficient to compute the bonds from atom distances. The algorithm implemented by VMD is simple: if the distance between two atoms is shorter than 0.6 times the sum of their radii, those atoms are considered to have a bond. In the example coordinate file, that atoms are 1.0 units apart and we have set the radius to 0.85. Thus 0.6 * (0.85+0.85) = 1.02 which is larger than 1.0, so there will be a bond between each atom and its direct neighbor(s):

mol bondsrecalc top

This step does not assign any bond type (or bond order) information, but rather define it as "unknown". In many classical MD packages, the bond type is derived from the atom types that form the bond. TopoTools contains a script that helps with assigning bond types according to this convention. In this specific system, all atom types are the same, so there should be only one bond type:

topo retypebonds

Now all bonds have been given the type "A-A" based on being between bonds of two atoms with type "A". The force field used for this example also requires the definition of an angle potential between neighboring bonds. This is done via:

topo guessangles

This step defines an angle between each triple of atoms that consists of two bonds sharing one atom. The angle typing is done in automatically in this step and in a similar fashion as in the bond typing step, i.e. all angles should have an "A-A-A" symbolic type.

As a final step we ask VMD to re-analyze the molecular data since we have changed bonds and angles and we need to update any data that VMD derives from that information, e.g. which atoms belong to a single molecule.

mol reanalyze top

Finally we need to set the simulation cell dimensions and write out the final state as a data file for the desired "angle" atom style:

pbc set {100.0 100.0 100.0 90.0 90.0 90.0}

topo writelammpsdata data.step1a angle

The last script command quit terminates the VMD session, since this script is meant to be run in batch mode via:

vmd -dispdev text -e simplexyztodata.tcl

The attached script file in.step1a is an example for how to use the resulting data file in a LAMMPS calculation. This input will do a very short MD run with random velocities to randomized the coordinates and then request a minimization. The restart written at the end of of the calculation can be converted with the restart2data tool to a data file, data.restart-1a, and this can be read into VMD with:

topo readlammpsdata data.restart-1a angle

Step 1b: Some of the atoms have a different type

Now we want to look at a model polymer that has different atom types at the beginning and end of the chains and we discuss the twotypesxyztodata.tcl script. This is similar to the previous one, so only the differences are explained. We use the same coordinate file and most of the initial steps are the same, but we need to define a second selection to change the properties of the two atoms at the beginning and end of the 11-atom chain (which are listed in VMD as index 0 and 10, respectively):set sel2 [atomselect top {index 0 10}]

$sel2 set radius 0.95

$sel2 set name B

$sel2 set type B

$sel2 set mass 1.5

The bond detection also works the same, however, after retyping we should get two bond types, "A-A" and "A-B" instead of one. A list of bond types can be retrieved through:

topo bondtypenames

This will also indicate the order in which those symbolic bond type names will be mapped to numbers starting from 1 in the data file. This output of the script indicates that bond type 1 will be set for bonds between one atom of type A and one atom of type B, which bond type 2 will be assigned for bonds between two atoms of type A. Similarly, we'll get two angletypes, "A-A-A" and "A-A-B", respectively.

The remainder of the script is equivalent to the previous one, however, to run the calculation in LAMMPS some more parameters will need to be defined, so a new input script, in.step1b, is required.

Step 1c: Generate multiple molecules from one coordinate file

In this example we will read in the same data file as before, but want to create a topology that has two molecules instead of one, with one of them being translated relative to the other. The corresponding VMD Tcl script, twochainsfromone.tcl, starts like the script in step 1a, but after reading in the coordinate data and assigning its properties, we need to define a new "molecule" in VMD that is able to hold all coordinate information, since VMD does not allow to change the number of atoms in a data set after it has been created. Please note that in "VMD-speak" a "molecule" refers to the whole data set, regardless of how many actual molecules in contains. The animate command reserves storage for coordinate data, which is not automatically allocated unless coordinate data is read from a file.set natoms1 [molinfo top get numatoms]

set natoms2 [expr {$natoms1 * 2}]

mol new atoms $natoms2

animate dup top

Now we collect all information from the old molecule that we want to transfer to the new molecule in a list.

set proplist {name type radius resname mass x y z}

set atomdata1 [$sel get $proplist]

Then we change the original molecule's coordinates by translating it 4 distance units in x-direction:

$sel moveby {4.0 0.0 0.0}

and retrieve the same information to a new list:

set atomdata2 [$sel get $proplist]

Finally we create a new selection for all new atoms, join the two lists of atom properties and assign them to the atoms in the new molecule:

set sel2 [atomselect top all]

set atomdata [concat $atomdata1 $atomdata2]

$sel2 set $proplist $atomdata

From here on the script is equivalent to the previous ones and the MD input in.step1c is essentiall the same as in.step1a, except for the changes in file names, e.g. for the read_data command.

Step 1d: Automatic bond detection does not work