Ten ways to translate an atom

Usually, this blog involves a lot of self-praise. This entry is none the better: I usually find programs and their interfaces attractive when they offer many ways to achieve a goal. In C++ (also many other programming languages) for instance there are quite a number of ways to implement a loop because it contains many different programming paradigms: using for, using while, with recursive functions, or using std::list to just name a few.

Therefore, I would like to enumerate a few ways using MoleCuilder to translate an atom. So, the goal of all the following commands is to take an xyz file containing the coordinates of three atoms of a single water molecule centered in a domain of edge length 10 angström and translate the oxygen atom by one angström in the direction of the X axis.

You might want to consider this a riddle to tackle by yourself. Therefore, as a first I give you the list of commands that perform the actual moving of the oxygen (although more may be employed to recover from additional side-effects):

1. translate-atom
2. translate-molecules (variant a)
3. translate-molecules (variant b)
4. random-perturbation
5. fill-volume
6. fill-surface
7. remove-atom and add-atom
8. mirror-atoms
9. scale-box
10. stretch-bond and change-bond-angle

So, let's start with the actual solutions.

1. translate-atom

The first one is easy: select the atom and translate it.

molecuilder --input water_translateatom.xyz -o xyz \
    --load water.xyz \
    --select-atom-by-element O \
    --translate-atoms --position "1,0,0"

2. translate-molecules (and translate all back)

If we translate the whole molecule, then we do translate the oxygen but also its two hydrogens. Therefore, we have to translate the latter back to their old positions.

molecuilder --input water_translatemolecule1_stage1.xyz -o xyz \
    --load water.xyz \
    --select-all-molecules \
    --translate-molecules --position "1,0,0"
molecuilder --input water_translatemolecule1.xyz -o xyz \
    --load water_translatemolecule1_stage1.xyz \
    --select-atom-by-element H \
    --translate-atoms --position "-1,0,0"

As we need the position option twice (once for translate-atoms and once for translate-molecules), we have to split this action into two calls to molecuilder with an intermediate input file. This is only needed when using the CLI and not with any of the other interfaces.

3. translate-molecules (but have oxygen atom as single molecule by cutting bonds)

However, it is up to us to decide what a molecule looks like. As the xyz format does not contain bond information, the loaded water molecule is treated as a single molecule initially but only so because all atoms in a file go into the same molecule regardless of their bonds. Only update-molecule checks the bond structure and splits them apart. Hence, after calling it, we have three single molecules and may easily translate the oxygen "molecule". What if bonding information was contained? Simply remove it by destroy-adjacency.

molecuilder --input water_translatemolecule2.xyz -o xyz \
    --load water.xyz \
    --update-molecules \
    --select-molecules-by-formula O \
    --translate-molecules --position "1,0,0"

4. random-perturbation (with specific random number generator)

This is tricky as you have to remember that there are many random number distributions to choose from. Hence, simply pick one giving integers uniformly in the range [0,2], i.e. from the discrete set {0,1,2}, and pick the right seed such that we pick {1,0,0} for the "random" translation of the oxygen. (Naturally, we use a little script that runs through the seeds and stops when it has found an appropriate one.)

molecuilder --input water_randomperturbation.xyz -o xyz \
    --load water.xyz \
    --set-random-number-distribution "uniform_int" --random-number-distribution-parameters "min=0;max=2;" \
    --set-random-number-engine "mt19937" --random-number-engine-parameters "seed=120;" \
    --select-atom-by-element O \
    --random-perturbation 1.

5. fill-volume

Using fill-volume is very hacky as you simply have to guess what point inside the cylinder is taken. Note that a single filling point is not implemented so far. Hence, we use 12, get 4 valid filling points and use one of them. The three other copies of the water molecule are removed. Naturally, as the whole molecule was moved we have to shift back the hydrogens (or use update-molecule).

molecuilder --input water_fillvolume.xyz -o xyz \
    --load water.xyz \
    --create-shape --shape-name cyl --shape-type cylinder --stretch "1,1,1" --translation "5.28,6,5.628" \
    --select-shape-by-name cyl \
    --select-all-molecules \
    --fill-volume --count 12 \
    --clear-molecule-selection \
    --select-molecule-by-order 1 2 3 \
    --remove-molecule \
    --select-atom-by-element H \
    --translate-atoms --position "-1,0,0"

6. fill-surface

Using fill-surface is similarly hacky but a bit simpler as we may specify a single filling point on the surface of a sphere.

molecuilder --input water_fillsurface.xyz -o xyz \
    --load water.xyz \
    --create-shape --shape-name ball --shape-type sphere --stretch "1,1,1" --translation "6,6,5" \
    --select-shape-by-name ball \
    --select-all-molecules \
    --fill-surface --Alignment-Axis "0,1,0" --count 1 \
    --select-atom-by-element H \
    --translate-atoms --position "-1,0,0"

7. remove-atom and add-atom

As atoms are indistinguishable, we may also simply remove the oxygen at the wrong spot and add a new one at the right spot.

molecuilder --input water_removeaddatom.xyz -o xyz \
    --load water.xyz \
    --select-atom-by-element O \
    --remove-atom \
    --add-atom O --domain-position "6,5,4.628"

8. mirror-atoms

Mirroring atoms works only on the selected ones. Hence, we select the oxygen and choose a suitable mirror plane. That's all.

molecuilder --input water_mirroratoms.xyz -o xyz \
    --load water.xyz \
    --select-atom-by-element O \
    --mirror-atoms "1,0,0" --plane-offset 5.5

9. scale-box

Scaling the box scales both the simulation domain and all atomic positions such that their relative positions inside the box remain the same. Doing this, we have to translate each hydrogen back in a very specific manner. Note that we use translate-atom and mirror-atoms as we cannot use the same action twice on the command-line (this does not hold for any of the other interfaces!). Hence, we circumvent this by using mirror-atoms for the second translation.

molecuilder --input water_scalebox.xyz -o xyz \
    --load water.xyz \
    --change-box "10,0,10,0,0,10" \
    --scale-box "1.2,1,1" \
    --select-atom-by-id 0 \
    --translate-atoms --position "-1.144,0,0" \
    --invert-atoms \
    --unselect-atom-by-element O \
    --mirror-atoms "1,0,0" --plane-offset 4.708

10. stretch-bond and change-bond-angle

We may simply use stretch-bond and change-bond-angle to set both OH bond angles and then the HOH bond angle to value of the molecule after the O atom has been shifted. Note that this only works from v1.6.1 onward as before the bond graph is not treated when stretching bonds or changing angles (instead a cut plane is used whose normal vector is the bond vector). Note that we obtained the values for bond lengths and bond angle using the new geometry objects that are also introduced with v1.6.1..

molecuilder --input water_stretchbond_stage1.data -o tremolo \
    --set-tremolo-atomdata "Id type x=3 neighbors=4" \
    --load water.xyz -I \
    --select-atom-by-id 2 0 --stretch-bond 0.624310820025
molecuilder --input water_stretchbond.pdb -o pdb \
    --load water_stretchbond_stage1.data \
    --select-atom-by-id 2 1 --stretch-bond 1.80824887668 \
    --select-all-atoms --change-bond-angle 45.3788358886 \
    --select-all-molecules --rotate-around-self 70.008595817 --axis "0,1,0" \
    --add-empty-boundary "4.28,5,4.628"
molecuilder --input water_stretchbond.pdb -o xyz

This is a little tricky because of the limited output precision of the configuration files: PDB offers only three digits fixed point precision. Tremolo offers more but std::iostream precision is usually 6 digits at most. This is not enough for the subsequent rotations which need to be exact for 6 digits as well. Hence, we need the third step where we strip down precision to 3 digits by outputting in PDB format and then converting back toXYZ.

And that's all I could come up with so far. Have some more ideas? Mail me and I'll add them to the blog post.