Bouke P. van Eijck, Department of Crystal and Structural Chemistry, Bijvoet Center for Biomolecular Research, Utrecht University, Padualaan 8, 3584 CH Utrecht, The Netherlands. E-mail: or

For information on the standard UPACK package, see the main manual .
The ab initio option is implemented beginning with version 7.1, April 2006.


Code words
Ab initio charges
Ab initio energies and forces
Example: methyd


The possibility to calculate ab-initio charges and forces during the calculations has now been incorporated in UPACK. To this end the external programs MOLDEN [1] and GAMESS-UK [2] are used as subroutines. For copyright reasons we cannot include these programs in the UPACK package, each user has to make his own arrangements.

The primary function is to create ab-initio charges. GAMESS produces wave function output which is read by MOLDEN. This program then calculates charges by fitting to the electrostatic potential. Additions to this program enable the user to obtain distributed multipoles (written by Wijnand Mooij for use with the XTINKER package, not included) and to obtain charges at non-atomic sites (written by Bouke van Eijck for use in UPACK). UPACK does not have the facilities to employ distributed multipoles, but additional charge centers might (at least in principle) be able to attain the same goal.

A second, unfortunately rather computer-intensive, application is to incorporate forces from the gradient vector and the hessian matrix [3] of the free molecule(s) calculated by GAMESS. This option may use a different basis set (usually smaller to reduce the computer time) compared to the charge calculation.

In both cases the calling of GAMESS may be repeated if certain geometry parameters deviate so much from their starting values that the initial results become of doubtful reliability.


This feature uses the directories upack/abinitio/molden and upack/abinitio/gamess.

Obtain the UPACK-XTINKER version of MOLDEN from Gijs Schaftenaar. Create an executable, call it molden.e and place it in the directory upack/abinitio/molden. Read readme.1 and readme.2 in upack/abinitio/molden. In particular, go to the subdirectory testfit which contains an output from GAMESS-UK called fromgamess. To find multipoles, call:
molden-dma fromgamess
To find charges, call:
molden-esp fromgamess
Compare the resulting files dma.out and esp.out with dma.res and esp.res, respectively.

In upack/abinitio/gamess, place the GAMESS-UK executable which must be called gamess.e. It should be possible to use another ab-initio program, provided that the wave function can be read by MOLDEN and that the program can deliver first and second derivatives of the energy to the Cartesian atomic coordinates. For the latter application some adaptions in UPACK will be necessary, at least in the formats in the subroutine leesg.f which reads the ab-initio output file.

One example is the use of the DFTB method from theoretical chemistry, to replace GAMESS-UK. This program can be obtained from the SCM company, Amsterdam. It was tested in 2016 and at that time it was not sufficiently developed for satisfactory results.

Codeword dftb in prep, enters through bter1, btra1, bter2, btra2, bter3, btra3.

The subroutine dftb is defined in func3g and is called in prep and func3g.

List of files containing additional subroutines


As func3, with call of subroutine ygames


Calls GAMESS, calculates energy and forces


Removes translational and rotational parts and returns cleaned matrices.


Reads and analyzes GAMESS-output


The calculations are performed with the UPACK program prep. When this feature is active, a file may be present to define the internal coordinates. This file may be placed in the main directory upack/subst.
It contains the following data (for one group, atom numbering as in topology file):

naib: the number of bonds. Then for every bond 1... naib:

iaib, jaib: corresponding atom numbers

naih: the number of angles. Then for every angle 1... naib:

iaih, jaih, kaik: corresponding atom numbers

naid: the number of dihedrals. Then for every dihedral 1... naib:

iaid, jaid, kaid, laid:: corresponding atom numbers

Make sure that naid + jaidi + naid equals 3 times the number of atoms - 6 times the number of molecules. 

The relevant code words for ab initio work are:

bintra btra1 btra2 btra3 {none none none}
Basis set for intramolecular energy: initial calculation, energy minimization, and final calculation, respectively. None means a standard force field calculation.
The basis sets btra1 and btra2 must be the same except when btra1 = import (see below).

binter bter1 bter2 bter3 {none none none}
Same for intermolecular Coulomb interaction through ab-initio charges.

aintra aiofst {0}
Energy offset (atomic units) to bring energies to a manageable scale. Its exact value is unimportant; it should be estimated for all molecules in one group together.

maxshi saib saih said {0.1 5 5}
A new ab-initio calculation is done if a shift in bond length exceeds saib (A); same for bond angles (saih; degrees) or dihedral angles (said; degrees).

maxcyc maxcyc {15}
The maximum number of cycles of ab-initio calculations.

chadd ichad {0}
ichad = 0: charges from topology file (only if the code word binter is not set)
ichad = 1: charges from file or, but only for real atoms (default if binter is set; see below)
ichad = 2: charges from file or, also for virtual atoms (see below).

hess ihess {0}
If ihess = 1, the final ab initio energy calcution with basis bintra btra3 will also produce derivatives, and therefore take longer time. This is useful if the energy minimization has to be continued in a subsequent run. It can also be used (although ihess = 1 may hardly be worth the additional expense) if lattice vibrations are desired:
ihess = 0 or 1: the intramolecular hessian is taken from the last ab initio cycle
ihess = 2: it is taken from a previous calculation through file hessian

Ab initio charges

Upon calling from UPACK, MOLDEN creates a file, containing the charges. For each distinct molecule there will be a block of data in XYZ format:

nxyz: the number of atoms (real + virtual)

titel: comment. Then for every atom 1... nxyz:

atom name (X for virtual atoms), various data (not used) with last entry: charge

Apart from the virtual atoms, the atom order is the same as in the topology file. In later work these charges can be re-used without the time-consuming ab initio calculations. In the absence of, the program searches for a file called This feature can be used to compare several sets of charges without calling MOLDEN again and again.


We recall the possibility to introduce charges on "virtual atoms", i.e. on sites which are defined in terms of three real atoms a, b, c, as follows:

where p and q are dimensionless fractions and s is in A-1.
For each separate molecule (01, 02, 03, ...) these sites must be defined in the files mol01.def, mol02.def, ..... containing for each virtual atom:

a, b, p [,c, q [,s]]

where a, b, c are atom numbers as defined for that molecule (starting at number 1, which is different from the topology file for all molecules after the first one!). Absence of such a file means that for the corresponding molecule only the normal atomic sites are considered. So this feature can also be used to simply supersede the charges defined in the force field file. Note that the atom numbers start at number 1 for each separate molecule.

Ab initio energies and forces

Upon setting bintra, ab initio forces are calculated. Here the file is obligatory. The most time-consuming work is the calculation of the second derivatives. Therefore, it is useful to start with structures that are pre-optimized with an empirical force field that is designed to produce structures close to an ab initio energy minimum for a free molecule. Otherwise time may be wasted in adjusting trivial parameters, e.g. C-H distances taken from an experimental structure which are always too short.


Example: methyd

The probably non-existing substance methanol monohydrate, methyd, was chosen to enable tests which can be run in a manageable time. It has two different molecules in the topology file, and a situation with G = 2 in P21/c was chosen as a general example. 10000 structures were generated and energy-minimized in the OPLS force field, after which the 20 best structures were chosen for further study. Free energies in this force field are found in methyd/p21c/pack.33, methyd/p21c/pack.pr3 and methyd/p21c/allres.

Ab-initio energies for these 20 structures were calculated in the directory methyd/aidemo. To reproduce these results, copy the 20 SPF files from methyd/p21c into directory methyd/test and run there:

runppall.spfmethyd opls

Energies are given in allres (sort after E-tot in column 9), detailed results in files *.prep.pra, optimized structures in files *.prep.spf, charges in files * Finally, ab-initio derivatives are conserved in the files *.hessian.



1. G. Schaftenaar and J. H. Noordik, J. Comput.-Aided Mol. Des. 14 (2000) 123-134.
2. M. F. Guest, I. J. Bush, H. J. J. van Dam, P.Sherwood, J. M. H. Thomas, J. H. van Lenthe, R. W. A. Havenith and J. Kendrick, Mol.
Phys. 103 (2005), 719-747.
3. B. P. van Eijck, W. T. M. Mooij and J. Kroon, J. Comput. Chem. 22 (2001) 805-8154.

Back to beginning