THE UPACK - ab INITIO OPTION
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: vaneyck@chem.uu.nl
or vaneijck@xs4all.nl.
For
information on the standard UPACK package, see the main
manual .
The ab initio option is implemented beginning with version 7.1, April 2006.
Contents
Introduction
Installation
Overview
Code words
Ab initio charges
Ab initio energies and forces
Example: methyd
References
Introduction
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
|
yfunc3g.f |
As func3,
with call of subroutine ygames |
|
ygames.f |
Calls
GAMESS, calculates energy and forces |
|
yintcar.f |
Removes
translational and rotational parts and returns cleaned matrices. |
|
yleesg.f |
Reads and
analyzes GAMESS-output |
Overview
The
calculations are performed with the UPACK program prep.
When this feature is active, a file subst.int
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 molden.xyz or chadd.xyz, but only for real atoms (default if binter is set; see below)
ichad = 2: charges from file molden.xyz or chadd.xyz, 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
Upon
calling from UPACK, MOLDEN creates a file molden.xyz, 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 molden.xyz, the program searches
for a file called chadd.xyz. 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.
Upon
setting bintra, ab initio forces are calculated. Here the file
subst.int 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.
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.spf a methyd opls
format.up
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 *.prep.xyz. 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.