My last experiments have analysed different conformations of ethane and butane where I've used MolecCuilder to calculate the ground-state energy of various configurations of these two molecules. FOr example, butane I've rotated around its central bond in equidistant steps. This allowed me to compare the resulting potential barriers against the literature.

I'm still on this train of thought: Directly connected to these potential barriers is the question of stability. Given a molecule and a specific configuration, how can I assess the (thermal) stability from ab-initio calculations? And the other question I had in mind: How can I extend this analysis to other (small, organic) molecules?

I will look at the second question first: My inspiration was an article by a former colleague of mine ​[Israels, Maaß, Hamaekers, 2017] on the octet rule in chemical space, Basically, they are generating -- what I call -- saturated molecules using the octet rule (i.e. fully occupied electron orbitals).

It was a low-hanging fruit to add this as some of the necessary ingredients have already been present in MoleCuilder.

1. Saturation of non-hydrogen atoms by placing hydrogens along pre-computed vertices of polyhedra (e.g., for carbon with the four important sp-orbitals, we look at the four corners of a tetrahedra).
2. Encoding knowledge of occupied and unoccupied orbitals per element (e.g., the specific angle in water is mostly because two corners of a tetrahedra are picked)

I then created a larger python script using the pyMoleCuilder library that does the following steps:

1. Start with a "seed" atom in the center of the domain.
2. Saturate the atom.
3. Pick randomly a number of newly added hydrogens (up to a maximum, e.g., 3).
4. Replace this set of picked hydrogens by a single new non-hydrogen atom of a chosen element using the mean position of the hydrogens.
5. Rescale the bond between old and new non-hydrogen atom to typical tabularized bond length (in MoleCuilder called "bond table").
6. Step to the new non-hydrogen.
7. When enough non-hydrogen atoms have been added, stop. Otherwise, go to step 2.
8. Calculate ground-state energy or optimize the created configuration.

In order to then generate "all" molecules of a given set of elements (I've picked {C,N,O}), I don't use random picking but the powerset of all hydrogens up till a given number and execute the procedure for every powerset combination.

One small challenge is the naming. I wanted a representation of the bond graph in the name. Informatics knows all about graphs and there's also a format for specifying them: graph6. There's even tools and websites to generate all (connected) graphs up to certain vertex count: See the ​[House of Graphs] and ​[nauty], on searching and generating graphs. I've added a small action to MoleCuilder that gives the currently selected molecule in graph6 representation and that uses BFS to give the element list starting from the least-connected non-hydrogen atom (ignoring hydrogens).

Let me show you the result for N=1.

graph6	elements	chemical formula	molecule name	energy [Ht]	largest force [Ht]	num neg EV	configuration
Ds_	C	CH4	​Methane	-40.181	5.06882e-10 2.80043e-10 6.06817e-10	18
Cs	N	NH3	​Ammonia	-56.166	2.02605e-05 2.88662e-05 4.92678e-05	16
Bo	O	OH2	​Water	-75.985	8.22233e-09 1.18472e-08 2.90999e-11	14

You find the graph6 string (of full graph including hydrogens) first, then the list of elements, the chemical formula and the molecule's name. Moreover, I give the ground state energy (see below for details), the largest remaining force component per dimension and number of negative eigenvectors. Finally, a small graphical representation of the optimized structure (automated snapshot taken with molecuildergui using a small python script). The most important bit about this: It's completely automated. I'm just giving the parameters, run the whole experiment and the output is the above table.

Energy calculations have been done in 6-31G basis (mostly for reasons of computational ease, not accuracy) and either Closed-Shell Hartree-Fock (CLHF) or Moeller-Plesset Perturbation Theory second order (MBPT). As you see, O and C are well-optimized, N is not yet. However, my goal here is not to give accurate ground-state configurations and energies but just show what's possible.

The tables for N from 1 up till 3 are given in the attachments, where you also find the necessary script files.

It's been quite a while since analysis_ethane, where I suggested that a look at butane might be interesting, too, as its structure of conformers is richer. As you can see by the following screenshots from MoleCuilder, there are four conformers.

syn/eclipsed methyl	gauche	120/eclipsed	anti

However, the more complex an experiment, the more prone one is to making errors as well. I have made quite many: overwritten output files, stumbling over the limited precision of PDB files (use tremolo format if you're optimizing to very high precision), confusion when using the Python interface because of no "named arguments" support so far (contained in upcoming v1.6.1), copy&paste errors when building on the script I used for ethane.

But one thing especially impeded me: The quality of the optimization.

If one plans on investigating butane's rotational barriers, then one should have a fully annealed "anti" conformer state at hand.

But to make a long story short, here it is:

As you see in the conformer screenshots, I have done the rotation around the middle CC bond in the same manner as with ethane before. The conformer names relate to the following angles: syn (0°), gauche (60°), 120 (120°), anti (180°). I used two initial states: an optimized "syn" and an optimized "anti" state - both optimized to 1e-5 a.u. absolute error in energy in roughly 2000 steps. I have used 6-31G basis set and CLHF theory for optimization and rotation. If you'd like to reproduce the figure, using the docker image, both optimized conformer states are attached to this blog post including the update script to rotate around the CC bond and compute the energy. You need MoleCuilder version 1.6.1 which is soon to be released.

You notice that the two minima left and right of the central peak are not the same when starting from the "syn" state. They both correspond to the symmetric state of the same conformer. This is a clear indication that one needs to start from the lowest energy state, which is "anti". On the other hand, there is an energy difference in the central peak when starting from either state.

Comparing the above with experimental data at ​introorganicchemistry, you'll notice that the central peak is strongly overestimated by theory. However, a recent article by ​[Yirong Mo, JOCNote, 2012] states that experimental data and theoretical predictions do not match. This is why there are articles such as ​[Murcko, Castejon, Wiberg, 1996] and ​[Allinger, Fermann, Allen, Schaefer, JCP, 1997] calculating the energy differences to very high precision, estimating the syn/anti difference at around 0.087 Ht, while experimentally it is just 0.006 Ht, in chemical units 3.78 kcal/mole.

Here, we notice that the "syn" state is closer to the true energy difference as expected from theory, namely 0.0115 Ht. Note that ​[Yirong Mo, JOCNote, 2012] obtains roughly 0.01 Ht using Moeller-Plesset 2nd order and 6-31G(d) (with diffuse basis functions) where he optimized in each rotational step.

Having not optimized each step itself, I am satisfied with the results.

Since the release of version v1.6.0, MoleCuilder is available as a debian package (installable at least under Ubuntu 14.04).

Since that it is pretty easy to build a docker image such that MoleCuilder can be run inside a docker container. A Dockerfile is provided in the same blog post.

To make it even easier, I created a public repository at ​Docker, where the image can be pulled from directly: frederikheber:molecuilder/latest

In other words and assuming you have a working docker installation at hand, all you need to do to try the latest version is the following.

docker pull frederikheber/molecuilder:latest
docker run -it frederikheber/molecuilder:latest

This will open a root shell in the container. Type

molecuilder --version

And you'll see the version installed.

For computations you need both molecuilder_server and at least one instance of molecuilder_poolworker running. The Poolworker is the number cruncher. It's a modified version of ​MPQC. The Server accepts computation jobs from the MoleCuilder program and pushes them on to an idle worker. We start either one in the background, logging to a file:

molecuilder_server --controllerport 10024 --workerport 10025 2&>server.log &
molecuilder_poolworker --server 127.0.0.1:10025 --listen 10026 --hostname 127.0.0.1 2&>worker.log &

The server listens to MoleCuilder on port 10024, the Worker enrolls with the Server on port 10025, and finally the Worker accepts jobs on 10026. Simply choose any ports available. You can also add more than one worker at a time using GNU's parallel, which is not installed on the container,

apt-get install parallel
parallel --gnu molecuilder_poolworker --server 127.0.0.1:10025 --hostname 127.0.0.1 --listen ::: 10027 10028 &

Check the server.log file to see that indeed 3 workers are in the queue now.

Well, that's all. Would you like to calculate the energy of a methane molecule?

molecuilder \
  --input methane.data  \
  --set-output tremolo \
  --set-tremolo-atomdata "id type x=3 F=3 neighbors=4" --reset 1 \
  --set-parser-parameters mpqc --parser-parameters "basis=STO-3G;theory=CLHF;" \
  --add-atom C --domain-position "10,10,10" \
  --select-atom-by-element C --saturate-atoms \
  --select-all-atoms --optimize-structure --steps 10 --deltat 1. --server-address 127.0.0.1 --server-port 10024

This should bring the largest remaining force components to about 1e-5 a.u..

Don't be overwhelmed by the lot of typing and let's explain step by step:

• First, we specify an input file (which is empty but will contain the methane on exit).
• Next, we set the output format, namely tremolo.
• And afterwards, we need to specify the tremolo format to contain the id of the atom, its element, three spatial and three gradient coordinates and at most four (covalent) neighbor ids.
• The ab-initio calculations will be done at Closed-Shell Hartree Fock (CLHF) with a simple STO-3G basis set. As we use MPQC as number cruncher, we have to give these information to the parser mpqc.
• Then, finally, we create the methane molecule by adding a carbon atom and saturating it.
• Last, all atoms are selected and the structural optimization is performed for 10 steps with an initial step length of 1 a.u.. Notice that we give the port 10024, which we told the Server to use for listening to molecuilder's jobs.

Wanna take it from here? Increase the optimization steps. Choose another basis set, e.g. 6-31G (any of MPQC's sets is fine), or Moeller-Plesset Perturbation theory 2nd order (MBPT2), ...

Some pro tips:

• have one executable running per container, i.e. start you server as

  docker run --name "MoleCuilderServer" -p 10024:10024 -p 10025:10025 frederikheber/molecuilder:latest \
  /usr/bin/molecuilder_server --signal 2 15 --controllerport 10024 --workerport 10025
  

• similarly, start you workers each in a single container, too, where we automatically extract the servers IP address

  SERVERIP=`docker inspect MolecuilderServer | grep IPAddress | tail -n 1 | awk -F":" '{print $2}' | tr -d \"\ ,`
  docker run --name "MoleCuilderWorker" -p 10026:10026 frederikheber/molecuilder:latest \
  mpirun -np 1 /usr/bin/molecuilder_poolworker --signal 2 15 --server ${SERVERIP}:10025 --listen 10026
  

• use GNU's parallel to start e.g. 8 workers at once, assume the above is contained in script start_worker.sh

  parallel --gnu --delay 1 --ctrlc ./start_worker.sh ::: `seq 1 8`
  

• you can also set up docker to start containers automatically at startup, see ​here. But this seems to be changing in the docker API.
• install molecuilder on you local machine and use molecuildergui for the methane optimization example above.

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.

In today's blog entry we look at the simple ethane molecule. More precisely, we interested in the rotational change between two of its conformations: staggered and the eclipsed.

We are interested in the energy barrier between these two conformations.

Let's first build the ethane molecule from scratch. We fire up molecuildergui and add two carbon atoms about 1.6 Angström apart from one another using add-atom. Then, we add hydrogens with saturate-atoms after having selected both carbons (either by clicking on them) or using select-all-atoms. The result looks quite reasonable already but the structure is not optimized. This we make up leeway by optimize-structure with 500 steps at a deltat of 2.5, all force components are then below 7e-7 a.u.. Note that beforehand we have used set-parser-parameters of mpqc to have a basis set of 6-31G and theory to MBPT2, moreover used random-perturbation with set-random-number-distribution as uniform_01 to randomly perturb all positions by at most 0.05 Angström. The latter shakes the configuration out of any possible unstable equilibrium. We obtain a final energy of -79.38819 Ht.

Now we step on to the actual analysis of the conformational change. We use a little python script (see C2H6_rotate-bond.py​) using the python module pyMolecuilder of MoleCuilder (version 1.6.1). The script rotates via rotate-around-bond one end of the ethane molecule by a prescribed angle, where we use select-atom-by-element to select both carbon atoms and thus specifying the bond around which to rotate. In total we go through the 360 degrees in steps of 5 degrees. After each rotation, the energy of the molecule in vacuum is calculated using MPQC by the action triplet fragment-molecule, fragment-automation, analyse-fragment-results. Before each rotation, we step on to the next discrete time step using step-world-time such that we obtain the complete trajectory of the rotation on output via output-as.

In the figure we see nicely the change in energy between the staggered and the eclipsed conformation, each occurring thrice over the range of 360 degrees. Note that we plot against the lowest energy of -79.38819027 Ht, hence negative difference actually means a higher energy. From the graph we read off the energy difference as roughly 0.005 Ht, equivalent to 0.136 eV. At room temperature of 300 Kelvin, we have a average thermal energy of 0.026 eV. Hence, the eclipsed state's probability is about $\exp(-0.136eV/(k_B * 300 K)) = 0.0052$ with respect to the staggered one, i.e. ethane will remain mostly in the staggered state.

For a nicer and more elaborate discussion of the underlying chemistry, see ​here. There, it is suggested to look also at butane. And one should replace one hydrogen atom on either end of the ethane molecule by a methyl group.

But we will leave this for the next blog entry ...

To allow for an easier installation of the debian package, MoleCuilder now integrates levmar, vmg, JobMarket, CodePatterns, LinearAlgebra, and MPQC as ThirdParty packages into its distributable tarball, i.e. these dependencies no longer need to be installed extra. The only extra dependence is Qt and Qt3D (the latter being unofficially provided under DownloadRedistribution).

This makes it easy to place MoleCuilder into e.g. a ​Docker container and have its molecuilder_poolworker sweat away on all those fragment calculations up in some cloud, see the attached Dockerfile​. Moreover, JobMarket is now fully integrated, meaning that a distinct molecuilder_server process takes care of a host of poolworkers that can munch away on the calculations in parallel. These calculations naturally are provided by molecuilder acting as the controller in this server-worker-controller implementation.

Furthermore, the fitting of empirical potential has been significantly enhanced by combining each potential fit function (e.g. dihedral angle bond) with a distinct bond model. This allows to precisely match the bond model with suitable fragments and therefore obtain more general fits of the empirical potential parametrizations. In the upcoming version, this will also work with partial charges and then we will provide some examples.

We added some more complex code for saturating atoms properly. It takes fully into account already present bonds and tries to match to internally stored SphericalPointDistributions as best as possible. This gives really good initial results, i.e. building up molecules just by adding a few carbon, nitrogens and oxygens, when subsequent structure optimization relaxes everything to its proper equilibrium positions.

Moreover, we fixed the fragmentation working on just a subset of atoms, i.e. leaving out saturation water.

In the GUI, the Fragment list can now be properly sorted (and still clicked to get the fragment as set of selected atoms).

FragmentationAutomation will actually fail when the Server is unreachable and not spin forever.

ParseTremoloPotentials is gone, as it is superseded by ParseParticleParameters which is also responsible for parsing partial charges. See also SaveParticleParameters.

Note that for versions in between this and v1.5.0 we refrained from producing debian packages as these did not really bring any good change for the user. Starting from this version, it is again worth to update ... especially with what's coming up!

1.5.3? Where is 1.5.1 and 1.5.2?!

To make things clear from the start, it took quite some effort to produce a clean interface between the GUI and the rest of molecuilder. Hence, v1.5.1 and v1.5.2 are sort of versions in between.

v1.5.3 basically brings mostly internal changes in form of the QtObservedInstanceBoard, namely all of the information about atoms, molecules, ... required for visual representation are copied and updated as needed. This is to make sure that things just live long enough for the sloppy GUI.

Right now, 1% of all test cases do not work in case of using the GUI as promised with v1.5.0.

Apart from that we have:

• The FitParticleCharge action is renamed to FitPartialCharge.
• long range forces working when using vmg.
• interdistance option when using Fragmentation action allows now to combine fragments up to a given order.
• electronic charges are additionally smeared for long-range electrostatic calculations which improves accuracy.

It is done!

From this version onward we provide debian packages (Ubuntu 12.04) for convenience to the user. Please note that MoleCuilder depends on several other packages, three of which had to be manually compiled and are made available at DownloadRedistribution. Especially, Qt3D was quite hard to get to work.

Furthermore, the GUI is now thread-safe!

Qt works with threads, MoleCuilder so far did not. Sure this had to clash. Now, the interface between MoleCuilder and the GUI done with Qt has been secured with mutexes and clean structures and the GUI is safe for use. All the more because it is contained in the debian package. Have fun fooling around with it.

Note that there still might be some issues. True stable execution will probably come up with v1.5.1 but early-adopters are welcome to try it out. For the moment 1 percent of all guicheck testsuite tests (see last versions about that) are still allowed to fail.

We need to switch from ​Codepatterns v1.2.9 to CodePatterns v1.3.0. Most of the GUI tests are failing because of race conditions in the code. CodePatterns v1.3.0 brings the Observer/Observable up to multi-threaded speed by having mutexes in place.

Hence, in this version nothing is really changed but we only switch underlying library versions on a version update.

Prepare for v1.5.0 because then the GUI tests will be working.

This is a smaller update that will some time later make a difference.

First of all, a new action has been added that allows saturating any present atom, i.e. creating methane is as simple as adding a carbon atom and saturating it.

Second, we finally have something like a testsuite for the GUI which should bring this interface up to a usable level real soon. The tests are mostly failing but this should soon be corrected.

Last, some GUI stuff has been fixed and enhanced such as the info boxes that now reveal a bit more about the atom or molecule the mouse pointer is hovering over.

First of all, you can now simply flip a switch and molecules are no longer shown via their atoms and bonds but as a tesselated surface, see the image below of an amylose molecule filled with a regular grid of water molecules. This is especially useful with very large molecules (think e.g. a Tobacco Mosaic Virus or simply some large biomolecule) because a lot less triangles have to be drawn on the screen. Molecules can still be highlighted and selected by moving over or clicking onto the tesselated surface.

Notice also in the picture that the list of molecules is fully working (and hopefully) perfectly. Molecules are aggregated by their formula and can be individually selected or their visibility changed as well as for all of the same formula at once. Changing a molecule by adding an atom immediately changes the list. This rewrite became necessary as empty molecules are now automatically removed.

Furthermore, there's a number of smaller changes that should make the GUI use more easy going:

• the 3D view can now be changed via the keyboard (up/down, left/right, page-up/page-down and zooming with +/-)
• New RemoveMolecule action for convenience (before you had to select-molecules-atoms and remove-atoms)
• Selections can be pushed and popped which is useful inside MakroActions.
• Bonds are restored when undoing a remove action
• Failing actions no longer clean the queue and thus the destroy the stored session file
• when checking the command-line options for a specific action,

  .../molecuilder --help fragment-molecule
  

  is enough, actionname is no longer needed.

In Germany there is a very famous and long-running TV show for kids called "Die Sendung mit der Maus" (the show with the mouse ... and not to forget: the elephant). I must admit that I still watch it frequently because you learn a lot, e. g. how a stuff is produced in factories, how walnut oil is made, and much, much more. I have even got a "Maus"-style tear-off calendar with a question and an answer on the flipside for every day of the year!

One of these questions (from 19th of June 2014) was: Why do noodles become soft during cooking?

My first guess was that cooking breaks up the starch chains to some extent and that’s why the noodles are softer afterwards. The tear-off page’s flipside however taught me a different answer: They consist of milled durum and come dry out of the package. Thus, they are hard but also last for a long time. They get soaked with water during cooking and that’s why they are soft.

I immediately had a ball-and-stick model in my head of starch chains with and without water and thought of how cool it would be to make a molecular dynamics simulation to check whether the answer could be validated this way. So this is me performing this virtual experiment. And the best thing of all: As there is no bond breaking involved, I can use the BOSSANOVA scheme, detailed in my thesis, and, of course, implemented in MoleCuilder (v1.4.9 required).

Actually, I had intended to give the experiment here in full but adding it in markup is impossible and as html still very hard. Attaching files is however easy. Hence, find the PDF to the experiment in the attachments, along with the amylose molecule if you want to try it out for yourself.

With this version we safeguard that the userguide contains at least one example for each and every Action.

The greatest change is however that ... how to put it ... convexification of tesselated envelopes is finally working. This means that volumes of envelopes and thus the space a molecule occupies can be calculated. This is because a convex envelope's volume can be easily calculated via an inner point and general tetrahedrons. If the volume changes can be traced when transforming a non-convex envelope into a convex one, then we also know the volume of the original non-convex envelope. The volume changes however consist of either removing points ("filling up craters") or adding general tetrahedrons, both volume changes are calculable.

As an example, we give two pictures to illustrate these envelopes, here for the amylose molecule, also known as potato starch.

There is a video among the files attached to this blog post where the progress is shown and the above description of removing points and adding tetrahedrons should become a little bit clearer.

This is a very important change because knowing the volume a molecule takes up is essential in calculating densities and thus knowing how many molecules to fill into a domain. See the upcoming amylose example.

Last but not least trajectories stored in the atom class has been changed. Now, we store them as maps, not as vectors. This means that steps in the trajectory do not necessarily have to be strictly consecutive. For example, when performing a structure optimization (with debugging on) the extra trajectory steps generated for a subset of atoms are no problem anymore.

Starting with version 1.4.8 we check for each version that valgrind runs through without glitches (except for some python sessions where valgrind shows too many python errors).

However, greater impact is that the GUI has been made more useful with a lot of smaller changes:

• No more clicking "ok" buttons when the dialog does not contain any queries.
• bonds are now properly updated in the 3D view.
• tooltips are shown for every query in the dialog, explaining its purpose.
• homology tab is working.
• file dialogs finally know something about the file types they are looking for.

With v1.4.7 a better hydrogen saturation procedure was implemented.

So far, we dealt with each bond to be saturated one after the other and could handle up to bond degree of three. The problem with the old code was that placement of saturation hydrogen was independent of that for another bond. This could result in hydrogens being placed on top of each other or at least very close. This occurred for example with the SLES molecule when the sulfur bonds were saturated due to its high valency.

Now, we calculate optimal places for saturation globally per molecule prior to the fragmentation procedure. After fragments have been determined and cut bonds need to be saturated, the required positions of the hydrogen are looked up from a table created by this prior calculation. This ensures that hydrogen placed by the saturation of two different bonds are optimally far away from each other.

Note that hydrogen saturation is a perturbation of the underlying Fock matrix (in the case a Hartree-Fock solver is used). Hence, resulting energies are affected and will now be slightly different.

There is a new make target extracheck that runs when mpqc (​http://www.mpqc.org/) is installed and can found in the path: It compares resulting values for some molecules against stored ones. There, relative differences (for bond order of 3) of up to 1e-5 are admissable!

Furthermore, the version now compiles on Ubuntu 14.04 systems. Some changes were necessary due to new boost and gcc versions.

Last but not least a simple Structural Optimization procedure based on the fragmentation has been implemented. This is so far just using calculated forces and some fixed stepwidth dampening scheme to look for the minimum. Enhancements are coming up but will take some more time.

With v1.4.6 the src/Actions/FragmentationAction/FragmentationAutomationAction.cpp​ can be used even without having [doxygen:install.html JobMarket]. Fragment Jobs (i.e. each fragment is stored as an MPQC file, MPQC is executed on it and the resulting energies and forces are extracted from its output) are launched serially.

Hint: You might be tempted to specify fragment-executable as mpirun -np 6 mpqc in order to have at least multiple threads working on one fragment at a time. However, this will fail cause the option's value is checked whether it exists as a file (and there is no file called "mpirun -np 6 mpqc"). One way out is to create a small script runmpqc.sh as follows

#!/bin/sh

mpirun -np 6 mpqc $@

Make it executable and use it as fragment-executable.

Another feature is that the console log is now shown inside the GUI (there is an extra widget called Log). You can click on an atom's name, e.g. C09, in the text and the atom will be (un)selected to highlight it. This comes in handy when something does not work as expected and you want to check what the code does while having the molecular system in front of you.

Last but not least: Each Action should give STATUS messages about success or failure, visible especially in the GUI (the bar at the bottom of the screen is the status bar).

In this blog we will report on changes and highlight some of them in upcoming versions of MoleCuilder.

Most likely, we will also feature small experiments or tell about interesting example usages.