Automated fragmentation or calculating fragment energies automatically

First of all, for this example you need ​mpqc, which is a very reliable ab-initio quantum chemistry code featuring among others a very fast Closed-Shell Hartree-Fock implementation. We will use MPQC to calculate the fragment energies of fragments of a molecular system, see also FragmentingMolecules. Here, we will not stop at purely creating fragments but do the whole thing, from scratch to end.

Second of all, we need a molecular system: See the file attached to this page with an alkane with 20 carbon atoms. Download it.

Last but not least, you need to have at least version 1.4.6 of MoleCuilder ready and we are good to go.

The whole procedure consists of the following steps and you will recognize them in the code snippet below:

• loading the molecule
• set the domain by using a boundary of 5 angstroem around the molecule
• create a bond graph (with correct bond degrees)
• select all atoms, i.e. marking them for the fragmentation to come
• create the molecular fragments up to order 3
• using the fragmentation automation procedure to create so-called jobs (it's all under the hood, you just do not specify output-types such that no files are written) and to call MPQC on each one of them sequentially
• analyse the resulting fragment energies to calculate the (approximate) total energy of the molecule

./molecuilder \
  --input alkane-20.xyz \
  --add-empty-boundary 5,5,5 \
  --subgraph-dissection \
  --select-all-atoms \
  --fragment-molecule BondFragment \
    --distance 2. \
    --order 3 \
  --fragment-automation \
    --fragment-executable mpqc \
  --analyse-fragment-results

After a brief moment of calculation you find a file called BondFragment_Energy.dat in the same folder that should have the following contents:

level   energy_total    energy_nuclear_repulsion        energy_electron_coulomb energy_electron_exchange        energy_correlation      energy_overlap  energy_kinetic  energy_hcore
1       -801.5418257    0       0       0       0       0       0       0
2       -779.3145265    0       0       0       0       0       0       0
3       -779.3905043    0       0       0       0       0       0       0

You see three lines beginning with the level (or fragmentation order) up to the specified order 3 and next to it the total Hartree-Fock energy approximation up to it.

Also there are more columns with certain energy contributions but these are only accessible with the additional JobMarket package (see section under [doxygen:install.html CompilationGuidelines]).

A few other files have been written which we briefly explain here (leaving out the above specified prefix BondFragment, i.e. you can change it to whatever you like):

• _Energy.dat - approximations of the total energy and certain energy contributions
• _Eigenvalues.dat - the list of approximated eigenvalues per level (only available with JobMarket)
• _Eigenhistogram.dat - a histogram of the eigenvalues per level (only available with JobMarket)
• _Forces.dat - the energy gradients for each atom per level (if gradient calculation is enabled by "gradient=yes;" (by default) in the MPQC parser-parameters)
• _IndexedEnergy.dat - similar to _Energy.dat only per fragment, you get twice as many columns, first the calculated energy per fragment and then the contribution to the total energy.
• Adjacency.dat - the adjacency list for each atom, i.e. the calculated bond graph
• KeySets.dat - the atomic indices of each fragment sorted by fragment (first just one id, then lines with two, then come lines with three ids, equal to the maximum order specified).
• _Times.dat - the times required for the calculation per level, level 3 gives the total calculation time (currently not available without JobMarket, as time is not parsed).

Modifications

We'll go through a few more details here for fine-tuning the above.

Modifying the solver

We have used the default values used for an MPQC output file. What if you want to use a different basis set? We just modify the default MPQC parameters:

...
  --subgraph dissection \
  --set-parser-parameters mpqc \
    --parser-parameters "basis=6-31G;" \
  --select-all-atoms \
...

Checking that what's going on is correct

How do we know that the above larger basis set is actually used when creating the MPQC output files, these "jobs"?

We increase the verbosity prior to fragment-automation and it'll show the output files:

...
  --order 3 \
  --verbose 3 \
  --fragment-automation \
...

If you do this, you will notice that we actually perform a Moeller-Plesset calculation of second order. (Change it by adding "...;theory=CLHF;" to the parser-parameters)

And you will also notice that the calculation takes longer because of the larger basis set and the approximate total energies per level are different.

Another way is by checking the files with

...
  --order 3 \
  --output-types mpqc

Note that we then cannot do the automatic fragment calculation as the files are no longer stored internally (referred to as "under the hood" above) but written to the hard drive. But the very same files will actually be produced by the automatic procedure.

Afterwards you have 57 files BondFragment...in, ready for inspection.