Tutorial - Part 1

We start with a very minimal system, a linear chain of a simple bead-spring polymer with just one atom type and to be used with reduced units. The coordinates are provided in the chain.xyz file attached below. Although this file is extremely simple, it has two problems when loading it into VMD. Its atom type "1" will be interpreted by VMD as hydrogen and thus the automatic bond search will not work as its heuristics do not allow for bonds between hydrogen atoms. Furthermore the atom radius associated with hydrogen atoms (or other "normal" atoms) is to too large for the bond search to work correctly. We'll have to address these issues before telling VMD to guess the bonds and then generate additional topology information based on those guessed bonds. Step 1a: All atoms are of the same type

The first part of the tutorial focuses on using TopoTools for building simple topology files from plain coordinate data of simple models. Step 1: Building topologies for simple bead-spring polymers

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

The initial coordinates have some pretty high potential energy, so the MD input in.step1d has to be a bit more complicated to get via a few quenches (running some MD followed by minimization) to a somewhat relaxed structure.

In some cases the distance based automatic bond detection may fail like shown in the top picture on the right which was generated by using the script from step 1a on the spiral.xyz file. The script spiralchaintodata.tcl shows an example for how to overcome this kind of problem. This is a system where the particles are arranged in the shape of a spiral (cf. second image on the right), but towards the center the atoms from different turns are closer than the atoms that are bonded together. In this case we have to search for and alternate way to construct the list of bonds. For the system at hand we know that all atoms for a chain according to their order in the coordinate file. With that additional information we can assemble the list of bonds manually with the following VMD Tcl code:set blist {}

set natoms [molinfo top get numatoms]

for {set a 0; set b 1} {$b < $natoms} {incr a; incr b} {

lappend blist [list $a $b A-A]

}

topo setbondlist type $blist

Since there is only one type of atom, we also set the bond type to A-A directly. If some atom types would be changed, the bond types can always be updated with topo retypebonds. The rest of the script code is like before with one addition. In order to easy visualization we write out the .psf file, that can be read in at a later point to reset the bonding information to what we have computed with out script:animate write psf spiral.psf