Editing Molecules

Contents

To make a molecule, we first must make the atoms, then link them together, then do a molecular mechanics simulation to get the atoms to be in physically realistic positions. The result from this section of the tutorial will be a benzene molecule.

Making Atoms

The basic physical object we will use to make molecules is called a DesignAtom. Like the bouncing balls, it has physics defined for it, but instead of bouncing around the physics for a DesignAtom causes bonded DesignAtom's to move close together and unbonded DesignAtoms to repel one another until they are a certain distance apart. The equilibrium distance of two bonded DesignAtoms is close enough that Brenner's molecular modeling software will recognize that the two atoms are bonded, and the equilibrium distance of two unbonded DesignAtoms is far enough apart that the molecular modeling software will recognize that they are not bonded.

There are two keys that add atoms to the scene: comma and shift-comma. Shift-comma interacts with the user to determine what kind of atom to make. Comma makes whatever kind of atom shift-comma made last time, but without any user interaction. In either case, the new atom will be at the location of the mouse at the time the command was given, and the depth of the new atom will be equal to the average depth of the selection (if anything was selected) or the average depth of the scene (if there is anything in the scene) or the depth of the origin of the world coordinate system (if there is nothing at all). After the new atom is added, the new atom is selected and everything else is deselected.

Start this section of the tutorial with an empty scene. If there is anything in your scene, select it and delete it. Time should be stopped. To make sure time is stopped, you can press s. If the message is Stopping time, you are done. If instead the message is Starting time, then press s again to stop time.

Move the mouse to somewhere near the middle of the scene window and press shift-comma. You can then have the following interaction in the text window, where the user's inputs are in bold:

The configuration has these read/write fields:
Number  Field name      Field type  Field value
1       Atom to create  Atom        DesignAtom//C,SpaceFill
Give the number of a field to edit, or press 0 to exit: 1
This slot implements the protocol Atom.
The currently chosen factory is DesignAtom.
Do you want to change the factory and reset to the new factory's
default configuration? n
Editing the configuration...
The configuration has these read/write fields:
Number  Field name     Field type        Field value
1       Atomic Number  Positive Integer  6
2       Drawing style  Menu              Space Fill
Give the number of a field to edit, or press 0 to exit: 1
The old value of the slot is 6.
What do you want the new value to be? 1
The configuration has these read/write fields:
Number  Field name     Field type        Field value
1       Atomic Number  Positive Integer  1
2       Drawing style  Menu              Space Fill
Give the number of a field to edit, or press 0 to exit: 0
The configuration has these read/write fields:
Number  Field name      Field type  Field value
1       Atom to create  Atom        DesignAtom//H,SpaceFill
Give the number of a field to edit, or press 0 to exit: 0
This dialogue indicates that we're creating a DesignAtom and that we want it to have the atomic number 1, so it's a hydrogen. Just before the last 0 it reports the "Field value" as "DesignAtom//H,SpaceFill", which is confirmation that we configured it right. After the last 0 is entered, a small blinking white sphere appears in the scene graph. This is our new hydrogen atom.

Now let's create five more hydrogen atoms by positioning the mouse in various places and pressing comma. Arrange them in a largish hexagon, leaving us with a scene looking like this:

Now we'll make a ring of carbons in the middle. We could do this the same way we created the hydrogens, but we can do better. First we need to specify that we're creating carbons instead of hydrogens. Click the mouse on the black background; this deselects all atoms, and you will see that none of them are blinking any more. Then press "." (that is, period) and do the following dialogue in the text window to start creating carbons:
The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       Factory for replacement atoms  Atom        DesignAtom//H,SpaceFill
Give the number of a field to edit, or press 0 to exit: 1
This slot implements the protocol Atom.
The currently chosen factory is DesignAtom.
Do you want to change the factory and reset to the new factory's
default configuration? n
Editing the configuration...
The configuration has these read/write fields:
Number  Field name     Field type        Field value
1       Atomic Number  Positive Integer  1
2       Drawing style  Menu              Space Fill
Give the number of a field to edit, or press 0 to exit: 1
The old value of the slot is 1.
What do you want the new value to be? 6
The configuration has these read/write fields:
Number  Field name     Field type        Field value
1       Atomic Number  Positive Integer  6
2       Drawing style  Menu              Space Fill
Give the number of a field to edit, or press 0 to exit: 0
The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       Factory for replacement atoms  Atom        DesignAtom//C,SpaceFill
Give the number of a field to edit, or press 0 to exit: 0
Select something first.
In this case nothing was selected, so all the "." does is allow you to specify what type of atom future commands will create. The "." command normally modifies the selection, hence the complaint Select something first at the end of the dialogue. We will make more interesting use of the "." command later.

Now place the mouse in the center of the ring of hydrogens and press "6". This will create a hexagon of carbons, and the resulting display looks like this:

Making Bonds

To make a bond, use the l key ("ell", not "one"). If you select two atoms and press l, they'll become bonded together. Or, you can select one atom, hold the mouse over the other atom, and press l to bind them together. In this latter case the atom under the mouse will become selected, which makes it easy to bond atoms into chains or rings.

We will now use this to bind the hydrogen atoms to the corresponding carbons. Select a carbon, hold the mouse over a nearby hydrogen, and press "l". If the objects in the scene do not move, make sure time is started by pressing "s". You'll see the hydrogen and the carbon pull together:

Repeat for the other 5 hydrogens, giving this benzene molecule:

Saving and Restoring

We can make this benzene molecule more realistic by handing it to Brenner's molecular mechanics package, but this is a sufficiently drastic step that it is wise to save our work first. This is done by pressing control-s. The following dialogue can then happen in the text window:
The configuration has these read/write fields:
Number  Field name                             Field type  Field value
1       File name to write                     String      Unspecified.pdb
2       Write selection (true) or all (false)  Boolean     false
Give the number of a field to edit, or press 0 to exit: 1
The old value of the slot is Unspecified.pdb.
What do you want the new value to be? benzene.pdb
The configuration has these read/write fields:
Number  Field name                             Field type  Field value
1       File name to write                     String      benzene.pdb
2       Write selection (true) or all (false)  Boolean     false
Give the number of a field to edit, or press 0 to exit: 0
The configuration has these read/write fields:
Number  Field name                             Field type  Field value
1       File name to write                     String      benzene.pdb
2       Write selection (true) or all (false)  Boolean     false
Give the number of a field to edit, or press 0 to exit: 0
The need for the redundant 0 at the end is a minor bug.

Now you should have a file benzene.pdb. This is a normal PDB file. You can view it with Rasmol. If you have configured a pdb viewer (such as Rasmol) as a plugin for your browser, you can follow the benzene.pdb link and view the molecule; otherwise if you follow the link you'll simply see the text of the pdb file.

Now let's delete the benzene and try restoring. Click on the benzene with control-left-mouse to select it (this is the selection example we couldn't do earlier because we had no molecules), and press control-x to delete it. Type control-R to read a file, then answer the questions in the text window until you've read the benzene.pdb file back into the scene. The dialogue on the text window will look something like this:

The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       File name                      String      Unspecified file name
2       Select objects added to scene  Boolean     true
3       Where to put it                vector      (0.000000, 0.000000, 0.000000)
Give the number of a field to edit, or press 0 to exit: 1
The old value of the slot is Unspecified file name.
What do you want the new value to be? benzene.pdb
The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       File name                      String      benzene.pdb
2       Select objects added to scene  Boolean     true
3       Where to put it                vector      (0.000000, 0.000000, 0.000000)
Give the number of a field to edit, or press 0 to exit: 0
The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       File name                      String      benzene.pdb
2       Select objects added to scene  Boolean     true
3       Where to put it                vector      (0.000000, 0.000000, 0.000000)
Give the number of a field to edit, or press 0 to exit: 0
This causes the molecule to reappear in the scene, but we aren't done yet. If you hold the mouse over any atom in the molecule and press "/", you'll get output like this in the text window:
The cursor is pointing at object 13 which is a BoringAtom//C
with velocity (0.000000, 0.000000, 0.000000)
at position (2.657000, 19.898001, 0.000000)
with links to 12 14 7
The molecules we saved were DesignAtom's, but upon restoring they became BoringAtom's. This changes the behavior of the atoms; if you start time and drag around the atoms, you'll see that the BoringAtom's don't pull together the way the DesignAtom's did. We can fix this by converting them. Select the molecule and press "." (henceforth called "dot") to change the class of the atoms. Interact with the text window to indicate that you want the atoms to be DesignAtom's instead of BoringAtom's. This happens to be the default at this point, so to make the dialogue nontrivial let's convert the molecule to wireframe while we're at it:
The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       Factory for replacement atoms  Atom        DesignAtom//C,SpaceFill
Give the number of a field to edit, or press 0 to exit: 1
This slot implements the protocol Atom.
The currently chosen factory is DesignAtom.
Do you want to change the factory and reset to the new factory's
default configuration? n
Editing the configuration...
The configuration has these read/write fields:
Number  Field name     Field type        Field value
1       Atomic Number  Positive Integer  6
2       Drawing style  Menu              Space Fill
Give the number of a field to edit, or press 0 to exit: 2
The choices are:
0:  Wire Frame
1:  Space Fill
Which choice do you want? 0
The configuration has these read/write fields:
Number  Field name     Field type        Field value
1       Atomic Number  Positive Integer  6
2       Drawing style  Menu              Wire Frame
Give the number of a field to edit, or press 0 to exit: 0
The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       Factory for replacement atoms  Atom        DesignAtom//C,WireFrame
Give the number of a field to edit, or press 0 to exit: 0
You'll see this wireframe benzene molecule:
Now you may want to start time and drag the atoms around with the mouse to verify that they are behaving as DesignAtom's should.

Minimizing Energy

The molecular mechanics program is not integrated with the editing program, so to do molecular mechanics we have to set up a directory for running brennermd, save the scene in a special format, and run brennermd.

Making a Directory for brennermd

brennermd creates a lot of output files with fixed names, and the input file has a fixed name, so it's best to create a separate directory to store all these files. We'll assume here that the current directory for the Fungimol process has a subdirectory md that we will use when running brennermd. To create this directory, find a shell window and make sure it's directory is the same as the one for Fungimol, and give the command mkdir md.

Saving for Molecular Dynamics

To save the scene for brennermd, we must first change the types of the atoms in the scene. This gives us the option of specifying how the molecular dynamics package will treat them. There are, in principle, three choices for each atom: In practice, turning off the thermostat and strictly following the laws of physics introduces two problems. First, numerically integrating the laws of physics is inherently an approximation, so we can't exactly apply the laws of physics. The inaccuracies tend to accumulate and the symptom is that in long simulations, the temperature of the system increases without bound. Although this is an issue in general, the simulation in this tutorial is short enough that it is not a problem here.

The second problem is that if you've designed a system and you are simulating it the first time, it will typically have imperfectly placed atoms, which means a high potential energy which will quickly become a high kinetic energy and a too-hot system. If the components of your system can't take the heat, they'll disintegrate.

Therefore, by default each atom is controlled by a thermostat, and we will be using this default throughout this tutorial.

Before writing the input file to brennermd, we should indicate whether the atoms will be fixed, controlled by a thermostat, or dynamic by changing their factory again, like we did above just after restoring from the .pdb file. Select all of the atoms in the scene, and press dot. Change the factory for the atoms to BrennerAtom and accept the default motion flag:

The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       Factory for replacement atoms  Atom        DesignAtom//C,WireFrame
Give the number of a field to edit, or press 0 to exit: 1
This slot implements the protocol Atom.
The currently chosen factory is DesignAtom.
Do you want to change the factory and reset to the new factory's
default configuration? y
The available factories are:
Number  Name
0       BoringAtom
1       BrennerAtom
2       DesignAtom
What is the number of the factory you want? 1
Editing the configuration...
The configuration has these read/write fields:
Number  Field name                                               Field type        Field value
1       Atomic Number                                            Positive Integer  6
2       Drawing style                                            Menu              Space Fill
3       Motion Flag (1 = thermostat; 2 = fixed; other = moving)  Integer           1
Give the number of a field to edit, or press 0 to exit: 0
The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       Factory for replacement atoms  Atom        BrennerAtom//C,SpaceFill,THERMOSTAT
Give the number of a field to edit, or press 0 to exit: 0
Place the mouse over any atom and press "/" to confirm that this happened properly. The output will look like:
The cursor is pointing at object 15 which is a BrennerAtom//C,SpaceFill,THERMOSTAT
with velocity (0.000000, 0.000000, 0.000000)
at position (1.645125, 4.593740, -0.013181)
with links to 14 16 10
where the important part to check is the "BrennerAtom//...,THERMOSTAT".

Now we can save the file. If we're going to run brennermd from the md directory, we should save the file as md/coord.d. Type Control-S and save the scene to this file:

The configuration has these read/write fields:
Number  Field name                             Field type  Field value
1       File name to write                     String      Unspecified.pdb
2       Write selection (true) or all (false)  Boolean     false
Give the number of a field to edit, or press 0 to exit: 1
The old value of the slot is Unspecified.pdb.
What do you want the new value to be? md/coord.d
The configuration has these read/write fields:
Number  Field name                             Field type  Field value
1       File name to write                     String      md/coord.d
2       Write selection (true) or all (false)  Boolean     false
Give the number of a field to edit, or press 0 to exit: 0
The configuration has these read/write fields:
Number  Field name                             Field type      Field value
1       File name to write                     String          md/coord.d
2       Periodicity of space (Angstroms)       vector          (1000.000000, 1000.000000, 1000.000000)
3       Start time (Femtoseconds)              Float           0.000000
4       Time step (Femtoseconds)               Positive Float  0.500000
5       Write selection (true) or all (false)  Boolean         false
Give the number of a field to edit, or press 0 to exit: 0
Several things happened in passing here:

Running brennermd

Some basic parameters for the molecular dynamics run are defined in a file md/input.d. Typical text for this file looks like this:
20000 10 4 50  / # steps, # steps between data, thermostat, xmol writes 
0.6955 0.500 300.0  /Random # seed, neighbor list , temperature(K) 
 1                       / =1 REBO (C,H,Si,Ge), =2 tight-binding for C 
 6    12.0   51.2   2.28    / carbon Lennard-Jones parameters 
-1     1.0    8.6   2.81    / hydrogen 
-14    28.0   51.2   2.28   / silicon = carbon 
-10    20.0   47.0   2.72   / Neon 
-18    40.0  119.8   3.41   / Argon 
-36   131.0  164.0   3.83   / Krypton 
The parameters are:
Number of time steps, 20000 in the example above
How many time steps brennermd will do before it exits. The time per timestep was chosen to be 0.5 femtoseconds when we wrote coord.d, so this simulation will go for 10000 femtoseconds, or 10 picoseconds.
Steps between data, 10 in the example above
brennermd writes a file called output.d which includes lines for the total energy for the system and the total temperature at various times. Because this constant is 10, a line will be written to output.d every 10 timesteps, or 5 femtoseconds.
Thermostat, 4 in this example
How to control the temperature of the thermostatted atoms. Based on my reading of the code, the choices are:
Value Meaning
-1 Friction plus random noise. I have not tried this.
1 A Berndstein thermostat. We'll use this later in the tutorial.
2 Set all velocities to zero after each timestep. Things seem to barely move when I try this, so I do not use it.
3 An Evans-Hoover thermostat. I prefer using 1 over using 3; the Evans-Hoover thermostat seems to have too great a tendency to prevent things from moving.
4 No thermostat; just do the laws of physics. We'll use this very soon in the tutorial.
5 Probably the same as 1.
6 Energy minimization. One would expect a minimum-energy configuration to have a low temperature, but this does not appear to be the case reliably. Instead of energy minimization, I use thermostat 1 with a low target temperature, maybe 30 Kelvins.
8 This is used in the code but I never figured out what it means.
Beware that options 1, 3, 5, and 6 do not work when the temperature of the system is zero. The temperature of the system is normally zero when we write a fresh coord.d file, so after writing a fresh coord.d file you have to do a few timesteps with option 4 before turning on any real thermostat. See below for an example.
xmol writes, 50 in the example.
Number of time steps between writing the positions of the atoms in the system to xmol.d. This file is essentially a movie and Fungimol includes code that allows you to watch the movie. It is written incrementally during the molecular dynamics run, and Fungimol is robust when faced with a partially written file, so you can watch it during the run. Don't watch the movie continuously during the run unless you have multiple processors on your computer, since sharing one processor between movie watching and molecular dynamics will cause the molecular dynamics to take twice as long.
Random number seed, 0.6955 in the example
Used to determine the details of how the thermostat changes the velocities of the atoms. I have not tried changing it. I suppose it might be useful to run a simulation several times with different seeds to determine whether a design will behave randomly due to thermal noise.
neighbor list, 0.5 in the example
I don't know what this does.
Temperatore, 300.0 in the example
This is very useful. 300 Kelvins is a good temperature to run at because it's reasonably close to room temperature. If you want to do something like energy minimization but you want to be able to watch the movie, use a low temperature (maybe 3 Kelvins) and thermostat 1.
Everything else in the file
I haven't tried to understand these or run the program with different values for them.
Our immediate goal is to simulate the benzene molecule. First, as discussed above, we have to update coord.d so it has a nontrivial temperature. Revise input.d to run for 10 rounds with the no-op thermostat 4:
10 10 4 50  / # steps, # steps between data, thermostat, xmol writes 
0.6955 0.500 300.0  /Random # seed, neighbor list , temperature(K) 
 1                       / =1 REBO (C,H,Si,Ge), =2 tight-binding for C 
 6    12.0   51.2   2.28    / carbon Lennard-Jones parameters 
-1     1.0    8.6   2.81    / hydrogen 
-14    28.0   51.2   2.28   / silicon = carbon 
-10    20.0   47.0   2.72   / Neon 
-18    40.0  119.8   3.41   / Argon 
-36   131.0  164.0   3.83   / Krypton 
Run brennermd with the following shell commands:
$ cd md
$ rm -f output.d xmol.d
$ brennermd
 kend=  60
 TOTE:  -59.2041212
The command rm -f output.d xmol.d is present above because brennermd does not remove its output files before it begins work. If you do not remove them manually, then while the program is running, the incomplete output file generally consists of new data at the beginning followed by leftover old data at the end. This makes it awkward to monitor the progress of the molecular dynamics run.

It is interesting to look at the output.d file:

*CLASSICAL DYNAMICS SIMULATION OF  doc/tutorial/molecule/md/coord.d       

TOTAL STEPS=    100 DATA WRITTEN EVERY   10 STEPS

UNITS OF LENGTH, MASS, TIME AND ENERGY: 1 A  1 AMU  1 fs 103.6434 eV

LANGEVIN PARAMETERS:      300.000 k PSEED=  0.980083

        ENERGY(eV)     T      TSTEP(fs)  TIME(fs)

    10   -4.92309214    81.890    0.500    5.000

Hmm, it looks like the "total steps" count is wrong. Oh well. But more importantly, notice the plausible temperature: 82 Kelvins. This reasonable temperatured is an indication that our initial system was in a physically plausible state.

Our next step is to bring the system to room temperature. Change the input.d file to have thermostat 1 and 1000 timesteps and a temperature of 3 Kelvins, and run brennermd again:

[tim@infoscreen md]$ brennermd
 kend=  48
 kend=  24
 TOTE:  -51.4623215
 kend=  60
 TOTE:  -56.6558046
 TOTE:  -57.5855522
 TOTE:  -59.1336214
 ...
 TOTE:  -59.3366934
 TOTE:  -59.3364132
 TOTE:  -59.3374562
Note that the standard output from brennermd is generally pretty useless. We can't simply redirect the output to /dev/null because occasionally the program gets an error immediately after starting, prints the error message to standard output, and stops. I generally redirect the output to a temporary file and only look at the temporary file if the program exits suspiciously quickly.

Now we can do several things. First, let's look at the result and see if looks like a benzene molecule. Select all objects in the scene and delete them, then type Control-R in the Fungimol window and read md/coord.d:

The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       File name                      String      benzene.pdb
2       Select objects added to scene  Boolean     true
3       Where to put it                vector      (0.000000, 0.000000, 0.000000)
Give the number of a field to edit, or press 0 to exit: 1
The old value of the slot is benzene.pdb.
What do you want the new value to be? md/coord.d
The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       File name                      String      md/coord.d
2       Select objects added to scene  Boolean     true
3       Where to put it                vector      (0.000000, 0.000000, 0.000000)
Give the number of a field to edit, or press 0 to exit: 0
The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       File name                      String      md/coord.d
2       Select objects added to scene  Boolean     true
3       Where to put it                vector      (0.000000, 0.000000, 0.000000)
Give the number of a field to edit, or press 0 to exit: 0
Scene title is  md/coord.d                             
You'll see a nice regular benzene molecule:
If you position the cursor over any of these atoms and press "/", you'll see that the atoms are all BrennerAtom's, and that they have no bonds to each other. The lack of bonds is a problem because it prevents selecting molecules. We will soon fix the problem, but first convert all of the atoms to wireframe so we can see the fix happen. To keep unbonded atoms from being invisible, they are drawn as small spheres in wireframe mode:
Now you can infer the bonds by selecting all of the atoms and pressing "b". You'll see the wireframe benzene again:

Another thing to do is to check that the simulation brought the benzene to the desired temperature. Look at md/output.d; I see a temperature at the final timestep of 323 Kelvins, which is fine.

Another thing to do is watch the movie resulting from the molecular dynamics run. Select and delete everything in the scene, type Control-M to start the movie, and interact with the text window to indicate the movie you want. The default filename happens to be ideal:

The configuration has these read/write fields:
Number  Field name                     Field type  Field value
1       File name                      String      md/xmol.d
2       Repeat                         Boolean     true
3       Select objects added to scene  Boolean     true
4       Where to put it                vector      (0.000000, 0.000000, 0.000000)
Give the number of a field to edit, or press 0 to exit: 0
The movie looks like the benzene atom pictured above, except it wiggles around a bit. To stop the movie, select all of the atoms animated by the movie and delete them. This leaves the scene empty, which is a suitable state for starting the next section of the tutorial
Copyright 2000 Tim Freeman <tim@infoscreen.com>