THE UPACK PROGRAM PACKAGE

Version   1, 1994. Specially written for hexapyranoses
Version   2, 1996. General molecules, more space groups possible, polarization.

Version   3, 1998. Code-directed input, automatic threshold setting, force field through code numbers.

Version   4, 1999. More than one independent molecule possible.

Version   5, 2000. Treatment of lattice vibrations possible.

Version   6, 2004. Automatic reduction of high-symmetry space groups to lower symmetry.

Version   7, 2006. Charges on non-atomic sites possible.

Version   8, 2010. Complete update.
Version   9, 2011. Topology structure improved.
Version 10, 2012. Rigid structures possible, double precision throughout.

Contents

Definitions
Procedure
Coordinate systems
Energy calculations

Overview

Script files
Input files
Output files
Simple force field

Preliminaries
The search
The experimental structure
Example: ethanol
Example: benzene
Preparations for further work

Space group file
Codes file
Molecular connectivity file
Molecular topology file
Molecular geometry files
Force field file

Program pack12
Program pack3
Program dist
Program prep
Program crysmd

Introduction

The program suite UPACK (Utrecht Crystal Packer) is made available under the following conditions:

The programs are available free of charge for academic use only, excluding commercial purposes.
In publications where UPACK has been used acknowledgment must be made to the author.
No responsibility is taken for errors in the codes or documentation.

Definitions

The substance to be studied is characterized by a group of separate molecules. Such a group is the smallest possible repeating unit, for instance ABB in a dihydrate where A is a certain molecule and B is water. The constituent molecules are building blocks for the generation of hypothetical crystal structures. Generally several distinct molecular conformations are thinkable. It is possible to specify their geometries precisely, and to investigate a list of those that are deemed likely to occur. Alternatively, one may rely on automatic generation of possible conformations. This is especially useful for hydroxyl groups where the possible crystal conformations are hard to predict.

In the group a number of bonds may occur about which significant internal rotation is possible. We shall denote the number of rotatable hydroxyl groups by N_h and the number of other rotatable bonds (not connecting a hydrogen atom) by N_z. The symbol  will be used for the dihedral angles. Furthermore, let the group contain M_g molecules and A_g atoms. Of course, all these properties of the group are independent of the actual molecular conformation and crystallographic environment.

To generate a crystal, a number (G) of independent groups is placed in a certain space group which has S general equivalent positions. We shall exclude the possibility that a symmetry element of the group coincides with a symmetry element of the lattice: such a structure must be described in a space group with lower symmetry, disregarding that symmetry element. Thus the G groups always coincide with the asymmetric unit, and the total number of groups is S G. The number of really nonequivalent molecules in a unit cell will be denoted by Z''. In the example of the dihydrate ABB it can be seen that Z'' must be between 2 and 3G.

Procedure

Possible crystal structures are generated by placing molecules with varying positions and orientations into a unit cell of varying dimensions. This gives at most twelve (P-1) and at least eight (e.g. Pna2₁) parameters to be varied; this number can be considerably larger if the molecules are flexible or if more than one independent molecule must be placed. One of the cell parameters can be eliminated by using an approximate crystal density.

The calculations must be repeated for all probable space groups. Fortunately, relatively few space groups account for about 90% of the actually observed structures for organic molecules with Z'' = 1: P2₁/c, P-1, P2₁2₁2₁, P2₁,Pbca, C2/c, Pna2₁, Cc, Pca2₁, C2, P1, Pbcn, Pc. If sufficient computer power is available, one may extend the coverage by including P2₁/c, P1, P-1, P2₁2₁2₁, P2₁ with G = 2 [4]. If the substance is chiral, only P1, P2₁, P2₁2₁2₁, and C2 remain. For a group with internal symmetry the generated structures often belong to a space group with higher symmetry. The determination of the actual space group of these structures is not possible in UPACK; programs like PLATON [5] may be employed.

Coordinate systems

Except for input and output, where fractional coordinates may be used, the programs work with Cartesian coordinates in nanometer units. The crystallographic a-axis coincides with the Cartesian x-axis and the b-axis is positioned in the x-y plane. The crystallographic axes are projected on the Cartesian axes and the components a_x, b_y, c_z, b_x, c_x and c_y are used to characterize the unit cell.

The crystal is built from rigid molecules, defined in a local axes system with aid of three reference atoms for each molecule. The axes system is chosen in such a way that:

the plane through these three atoms is parallel to the x, y-plane,
the vector from atom 2 to atom 1 is parallel to the positive x-axis,
the y-coordinate of atom 3 is larger than the y-coordinate of atom 2,
the origin coincides with the average coordinates of the non-hydrogen atoms.

Energy minimization is the central technique of the method. The energy is always reported in kJ/mol for one group (although other units may be used in the force field files). The programs were developed with routines from the GROMOS molecular dynamics simulation program package [6] as basis, and we are grateful to Professor W. F. van Gunsteren for his permission to do so.

As discussed above, in the searching stage fast calculations are essential. A Lennard-Jones 12-6-1 potential is used here. Ewald summation may be used for the Coulomb and attractive r ^-6 dispersion terms, and for further work a Buckingham-type exponential repulsion is often preferable to the r ^-12 term. Atomic polarizabilities may also be included, although we have not found that refinement very helpful.

To allow fairly reliable calculations without an extremely large cutoff radius, the molecules can be split up into charge groups. The cutoff criterion is based on charge groups rather than on individual atoms. If the charge groups are chosen to have (almost) zero net charge, the electrostatic interaction decays at long distance with r ^-3 rather than r ^-1. To allow for small errors in the force-field setup, the total molecular charge is set to zero by adding equal amounts to each atomic charge.

The energy minimizations are performed by a number of steepest descents steps, followed by a more or less exhaustive conjugate gradients procedure. The use of pair lists for nonbonded interactions (as well as lists of vectors in reciprocal space for Ewald summations) is essential for stable convergence [7]. These lists are only updated upon a significant change of any geometric parameter.

Unavoidably, the same structures will be encountered repeatedly, probably with different choices for the unit cell parameters. Of course, it is very important to eliminate such equivalent structures by a clustering procedure. A general algorithm, based on comparison of radial distribution curves, is available.
 
 

Installation

The following directions apply to Linux systems. Little testing on other operating systems has been done, but all programs are in simple Fortran77 and should be fairly portable. Deviations from the standard are:

the use of include files.
the use of the do ... enddo syntax rather than labels (a more standard version can be produced on request).
the specification of some integers: two bytes for memory reduction, four bytes if necessary.
the use of ANSI-type color codes (in case of problems replace color.i: cp   no-color.i   color.i).

The program package is normally delivered in a tar archive: upack.tar.gz. Decompress the tar archive:
gunzip upack.tar.gz
Finally, unpack the resulting file in your home directory:
tar xvfk upack.tar
This should create a subdirectory upack directly under your home directory, with subdirectories progs, scripts, tops, data, upackdoc and abinitio. Note that the option k prevents overwriting existing files with the same names!

The subdirectory progs contains:

the programs and subprograms (all starting with the letter y)
the include file params.i, which contains most array dimensions (to be changed if required by your needs)
the other include files (all starting with the letter c and having the extension .i) containing the common blocks
various script files for compiling and linking, assuming the Fortran compiler gfortran

To compile and link the programs, go to subdirectory progs and call
compall
This should create the executables (*.e). An individual program or subroutine may be recompiled after a modification, e.g.:
f yprep
linkprep

The subdirectory tops contains the force field files (ending with .ff), the codes files (ending with .cod), and the space group file all.spg.

The subdirectory upackdoc contains the actual version of this manual.

The subdirectory abinitio contains a framework to combine UPACK with ab-initio calculations. Not all necessary programs are included for copyright reasons; see the ab initio manual.

The subdirectory data is meant to contain your data files. It should contain a subdirectory subst for every substance to be investigated, and here the topology file(s), e.g. subst.tps, should be placed. The best policy is to create sub-subdirectories for every space group to be considered. If you wish to use experimental data for comparison, they should be placed in a sub-subdirectory exp.

Some subdirectories under data are already included for demonstration and testing purposes. The example for the random search is ethanol which has two independent molecules in space group Pc. It is an exceptionally lucky case where the experimental structure is found quickly. Another simple example is benzene where a complete structure generation is demonstrated. The last examples are in data/glucsa, which is meant to demonstrate the now hardly ever used grid search, and in data/methyd for the the ab initio option.

List of files containing programs

ycodes.f	Creates framework for topology files
ypack1.f	Generates structures with one independent molecule
ypack2.f	Generates hydroxyl torsional angles
ypack12.f	Generates structures with more than one independent molecule
ypack3.f	Minimizes the energy
yclus.f	Removes equivalent structures with one independent molecule
yzoek.f	Searches for one equivalent structure with one independent molecule
ydist.f	Combines the latter two functions for the general case
yprep.f	Minimizes the energy for one structure
ycrysmd.f	Performs molecular dynamics calculations for one structure.
ymktop.f	Creates topology file with force field values.

List of files containing subroutines or functions

ybonde.f	cbond, choek, cdiha: calculates intramolecular energies
ycart.f	cart: constructs coordinates for a free molecule
ych.f	ch: calculates eigenvalues and eigenvectors (from EISPACK)
ycongra.f	congra: performs conjugate gradients energy minimization
yconn.f	conn: finds connected atoms from coordinates for a free molecule
ydclus.f	dclus: compares possibly equivalent structures
ydens.f	dens: estimates initial density for pack12
ydiff.f	diff: compares possibly equivalent structures
ydihan.f	dihan (function): returns the value of a selected dihedral angle
ydlist.f	dlist: makes list of nonbonded distances for dclus
yexpsym.f	expsym: expands symmetry operations to the full group
yfiles.f	files, fdat: makes mdi, spf, cif or dat type of coordinate file
yfunc0.f	func0: calls energy routines for a free rigid molecule
yfunc1.f	func1: calls energy routines for a rigid molecule in the crystal
yfunc2.f	func2: calls energy routines for a rigid molecule with rotatable hydroxyl groups in the crystal
yfunc3(g).f	func3: calls energy routines for a flexible molecule in the crystal
ygeom.f	geom: calculates and compares geometric parameters
ygetsym.f	getsym: finds all relevant symmetry elements
yident.f	ident: rearranges and relabels coordinates to conform with molecular topology
yintre.f	intre: calculates intramolecular energy
yintrot.f	introt: rotates a selected hydroxyl group
ykwald.f	kwald: calculates reciprocal part of Ewald correction
ylatvib.f	latvib: calculates lattice vibrations
yleesc.f	inpt (function): reads and stores card-directed input
yleesmd.f	leesmd: reads data for molecular dynamics run
yleesr.f	leesr: reads general space group information
yleess.f	leess: reads coordinates and expands equivalent positions
yleest.f	leest: reads molecular topology information
yleesx.f	leesx: reads geometry information, places hydrogen atoms
ylist.f	list: calculates nonbonded energy and forces from pair list
ylistk.f	listk: calculates reciprocal part of Ewald correction from list
ylistp.f	listp: as list, with polarization
yltw.f	nword (function), realstr, intstr: internal read procedure
ynonbe.f	nonbe: calculates nonbonded energy and forces
ynonbep.f	nonbep: as nonbe, with polarization
yopth.f	opth, listh: calculates energy for various hydroxyl torsional angles
yrunmd.f	runmd, stopcm: performs molecular dynamics run
ysecder.f	secder: calculates second derivatives for lattice vibrations
yshake.f	shake: keeps bond lengths constant during MD run
ysort.f	sort: sorts data in ascending order
yspepos.f	spepos: tests for special positions and reduces the symmetry if necessary
ystedes.f	stedes: performs steepest descents energy minimization
ysubs.f	random (function), norm, prod: various utilities
ytest.f	test: numerical check on calculation of derivatives
ytrans.f	trans: transforms to standard axes system
ytransf.f	transf, matmat, invbox: transforms to another crystal setting
ytri.f	tri: transforms to oblique coordinates and back

Overview

1. For every substance the molecule(s) in the group must be defined through a topology file.  Originally this file contained code numbers for every force field contribution (like bond angles, atomic charges, and so on). This scheme was designed to treat a large number of related compounds with the same force field. Nowadays it is obvious that such an approach is bound to fail, and that a force field should be designed for each compound individually. For that purpose it is easier to include bonded force field terms and charges directly in the topology file, and this is now possible.

2. For every conformation to be studied separately, a coordinate file for the constituent free molecule(s) has to be prepared. Note that often one input conformation can suffice, since rotational freedom about single bonds can be taken care of automatically.

3. The structures are generated by a random search for one or more independent molecules by the program pack12. Each molecule is placed with random orientation, random position, and (possibly) random rotatable dihedral angles in a crystal with random cell axes and cell angles. The energy of the generated structures is provisionally minimized with either rigid molecules or full geometry relaxation. A Lennard-Jones type force field is recommended here to avoid infinitely negative energies which can occur in exponential repulsion when atoms are placed close together.

4. The results are clustered by the program dist. This program also determines the number of crystallographically independent molecules (Z'').

5. The minimization is continued with the program pack3. Usually a first cycle is done with the same force field as used in the search, after which another force field may be preferred. After each cycle the results are clustered again.

6. Collecting and clustering the results from the different space groups delivers a list which may be considered as final, but preferably used for further study in a subsequent stage with a more sophisticated force field. There is the possibility of including polarization, relaxing the space group symmetry by expanding to P1, and to perform a primitive molecular dynamics "shake-up". In this way the stability of the obtained crystal structures can be investigated. It is also possible to study the effect of lattice vibrations on thermodynamical quantities [8]. Preferably, using other program packages, atomic multipoles and/or ab-initio energies should be used. A word of warning: it cannot be taken for granted that there is a one-to-one correspondence between the crystal structures in different force fields.

If the experimental structure is known, one will want to know whether or not it was found by the packing procedure. The program dist can be used for this purpose too, if the experimental structure is first energy-minimized with the program prep.

Molecular dynamics calculations can be carried out individually for one structure with aid of the program crysmd.

Script files

Files for input and output have fixed names in the various programs, to be described in detail below. These file names should be linked to actual files present on disk. To this end a number of C-shell scripts is available which link the necessary files and execute the programs. In this manual a certain convention for directory structure and file names will be used; other users may prefer to create a different system.

File names are given in brown color. Instructions to be given for executing the script files are red. Many programs require code-directed input: on each input line a code word is followed by a number of numerical values which are identified by colored names. To aid in understanding the Fortran sources, these names generally (but not quite always) correspond to identifiers used in the programs. The symbol # indicates that a certain number must be inserted in that place.

A script file to execute a certain program typically looks like this:
runprog s/a/... [filenum] subst [coord] ffname
Here runprog is one of the execute scripts, to be discussed in detail later. Some or all of the following parameters must be specified; entries between square brackets [...] are optional. A short explanation:
s/a/... : select the appropriate topology file, defined by the extension in subst.tps, subst.tpa, ...
filenum: a number (#) referring to a pack-type file (vide infra).
subst: a chosen name for the substance under study, which is also the directory name.
coord: the name of a coordinate file.
ffname: the name of the force field file ffname.ff to be used.

Input files

Most input files must have a precisely defined structure, but the ones that specify the course of the calculations are code-directed: parameters have default values unless specified by a line ("card") starting with one of a well-defined set of code words. An overview is given here, while a detailed description is given in a later section. Some explanation is also given as comments at the beginning of the subroutines.

The following files are more or less general, and have to be updated only incidentally:

all.spg: data on space group symmetry.
codname.cod: defines atom types and code numbers for various terms that contribute to the potential energy.
ffname.ff, various force field files to provide the actual numerical parameters corresponding to the code numbers.

subst.tps, subst.tpa: the topology file. The program codes may be helpful here; it requires the preparation of a file subst.con which gives the molecular connectivity as well as the atomic charges. The topology file may now also contain the bonded force field parameters and the charges.
cor.# (# = 001, 002, ...): atomic coordinates for every molecular conformation of the group to be packed. These files may be prepared by any relevant program; in UPACK the program prep can be used which will ask for the values of all dihedral angles.

subst.pa# (# = 1,2,3,...): input specifications for the various stages in the crystal structure generation. These files may be omitted if the default specifications are adequate.

coord: crystal cell data and atomic coordinates.
subst.pp: input specifications for prep.

Output files

The generated structures are collected in parameter files pack.#, the number # indicating the stage of the calculation in which the file was prepared:

# = 10: first stage (from pack12)
# = 20, 30, ...: subsequent stages (from pack3)

These files contain general information about the generated structures. Detailed atomic coordinates are delivered in corresponding files pack.19, pack.29, ..., respectively. These files generally contain equivalent solutions, which can be removed by the clustering program dist. This increases the file numbers by 3 (e.g., pack.30 is clustered to pack.33).

It is also possible to produce individual coordinate files of the hypothetical crystal structures in formats mdi, spf, cssr, cif or fdat (for further use in UPACK or other packages).

General output information is given on a file out.pri, which serves as a record of the simulation conditions. This file is usually linked to pack.pr#, distcl#, and so on. There is also output to the screen to signalize error conditions, and to obtain (optionally) more output for debugging purposes.

Simple force field

Since general force fields are of limited use, we may just as well use a very simple force field to construct a first version of the topology file, to be improved manually for each individual compound. It takes care of the most frequently occurring atom types:
 

HA	H in alcoholic OH
HZ	H in carboxylic OH
HN	H on nitrogen
HC	H on aliphatic C
HR	H on aromatic C
HW	H in water
HD	H on double bonded C
CP	C with single bonds (tetraedral)
CD	C with one double bond
CR	C in aromatic ring
NP	N with three bonds (pyramidal)
ND	N with a double bond
NR	N in aromatic ring
OD	O with double bond
OA	O in alcohol
OZ	O in carboxylic OH
OW	O in water
S	S with single bonds
SD	S with double bond
SR	S in aromatic ring
F
Cl
Br
I

Nonbonded parameters in the force field are from the OPLS Lennard-Jones potential (simple-lj.ff as in opls.ff) or from the Price-Williams Buckingham potential (simple-pw.ff as in will.ff); energy units are kcal/mol, distance units Angstrom. Torsional parameters are only set for aromatic systems; other ones have to be determined for each compound individually. Modern force fields depend on atomic charges from ab-initio quantumchemical calculations. For demonstration purposes only, the simple force field contains charge increments for end atoms (H, O) with some general usefulness. Of course, for real applications this is inadequate. Moreover, aromatic and conjugated systems are unique in principle: geometry parameters and charges have to be derived individually from theoretical chemistry.

In case the topology file contains also the information from bonded contribution and charges, now only nonbonded force field terms are needed. Then a "stripped" force field file is sufficient and can be used nearly universally. To this end there exist the files lj.ff and pw.ff.
 
 

Procedure

Preliminaries

1. Prepare the connectivity file(s) subst.con. This is helpful to obtain the molecular topology file(s) by running the program codes:
runcod subst.con [codname]
using a codes file codname.cod. The default codname is simple. Alternatives are opls and, for united atoms, unitat. If an atom type is missing, update the codes file and try again.

2. The output file out.tp can, after editing, be used as topology file. However, as discussed above, it is easier to have a topology file that includes bonded force field terms and charges for a specific molecule. For that purpose an auxiliary program mktop is now available:
runmk out.tp [forcefield]
using a force field file forcefield.ff. The default forcefield is simple-lj. This program produces the output file out.up which can be renamed to subst.tps (but see the next paragraph).
An additional new feature is that mktop can also use external charge input, e.g. from an ab-initio calculation. This information can be given in either .cssr format or .xyz format,  and the relevant file should be copied to charges.in. All previous charge information from force field or topology is discarded.

3. The output file out.up is a framework which usually needs editing. It is immediately necessary to inspect out.up for the warnings *****  which indicate the absence of a vital force field parameter. Replace the zero values by appropriate ones and check for further possible improvements. The simple force field produced in this way is good enough for testing purposes, but obviously inadequate for real work. So extensive further editing is necessary, especially for charges. Rename the edited file to subst.tp# (e.g. subst.tpsfor the simple force field, but replaced by subst.tpa, subst.tpb, ... in later improvements).

4. If necessary, update the space group file all.spg.

5. Decide which conformations are to be used for crystal packing, and prepare a coordinate file cor.# for each of them (# = 001, 002, 003, ...) for a free molecule. This can be done by your favorite molecule builder or by the UPACK program prep (for details see below):
runpp s subst none forcefield
The first call produces a file cart.out. Copy this file to cart.in and edit. It contains indices of dihedral angles that have to entered manually (degrees). To minimize the number of dihedrals to be assigned, the program omits hydrogen atoms in a first try; if no chain of four atoms can be found, it tries again with hydrogen atoms. If cart.in is present, the same execute script will do a regular calculation. Check the resulting geometry, especially the dihedral angles. Repeat the procedure for every desired conformation, and rename the output file(s) 0000.mdi or 0000.spf to cor.#, the coordinate file(s) to be used in the crystal structure generation.
 

The search

Decide which space groups are to be investigated. Avoid Buckingham-type force fields in this stage, as infinitely negative energies can occur when atoms are placed close together. Create subdirectories named after these space groups (with omission of slashes, e.g. subst/p21c). In every subdirectory a search must be carried out. Usually the default settings are adequate; otherwise input files subst.pa# must be created. Copy the coordinate files cor.# into each subdirectory. The random search consists of the following steps:

1. Generating and clustering structures:
runpa12 a subst forcefield
Output files are pack.10 (parameters) and pack.19 (coordinates); the parameter file must be clustered to pack.13:
rundist a 10 subst
To save computer time  a united-atom force field and/or a rigid-molecule search might be used. This is generally not recommended, and should only be done as a last resort.

2. Energy minimization: for instance, to continue with another force field:
runpa3 a13 subst opls
rundist a 20 subst
with final results on files pack.23 and pack.29.
 

The experimental structure

The program prep can also be used to minimize the energy of a given individual crystal packing. The crystal cell data and atomic coordinates are given in an input file coord. The space group information can be taken from this file or communicated through a codes file subst.pp.  Run the program:
runpp a subst coord forcefield
Output coordinate files are given in various formats.

An important application is the minimization of the energy of an experimental polymorph, if available. The geometry changes give an indication about the quality of the force field, as given by the parameter D. Also, the energy-minimized experimental structure should be present in the files from the search procedure. To this end the program dist can also be used:
rundist a filenum subst coord
where coord is the name of the file with the experimental coordinates (type spf). Since interatomic distances are strictly compared, this structure must be energy-minimized in the same force field as was used to create the hypothetical structures.
 

Preparations for further work

It is useful to prepare a set of coordinate files in the last run of the structure generation phase. These files must be collected for subsequent studies with other procedures. As these procedures are usually time-consuming, the number of coordinate files must be limited by selecting an energy range with respect to the global minimum encountered so far. A useful script is:
select.cssr  -12500 [directory]
Now all directories are scanned for CSSR files within a window of 12.5 kJ/mol (note that the range must be given in J/mol!). Hopefully the force field used is good enough to contain the experimental strtucture(s) within the chosen range. Adding the last keyword provides a directory into which the selected CSSR files are copied. The corresponding script for SPF files is:
select.spf  -12500 [directory]
In that directory further studies can be carried out. In UPACK the program prep can be used, where all present files with the following extensions can be studied in one run:
runppall.spf a subst forcefield
runppall.cssr a subst forcefield
runppall.cif a subst forcefield
 

Example: ethanol

Go to upack/data/etanol. Study the given connectivity file etanol.con.s, note the two charge groups and the definition of one dihedral angle. Make the topology file by running:
runcod etanol.con.s
runmk out.tp
Compare the resulting file out.up with etanol.tps.

To create a coordinate file many programs exist, for instance MOLDEN. In UPACK it can be done by running (for details see below):
rm cart.*
runpp s etanol none lj
cp cart.out cart.in
Edit the file cart.in and add the required values of the torsion angles (180, 60, 180, 300, 60, -60) at the end of the lines in that file. Then run again:
runpp s etanol none lj
Compare the resulting file 0000.mdi with cor.001.
Check for trivial errors:
runpp s etanol cor.001 lj
This confirms that the coordinate file can be optimized without problems.

For further demonstrations, go to the subdirectory etanol/exp. In the database we find the experimental structures ETANOL.CIF (Pc, Z"=2) and ETANOL01.CIF (a high-pressure form in P2₁/c, Z"=1). Perform energy minimizations:
runpp s etanol ETANOL.CIF lj
or:
runpp s etanol ETANOL01.CIF lj
This creates an energy-minimized structure in two 0001.* files. Observe that 0001.spf equals etanol.spf or etanol01.spf, respectively.
These files are the target structures for the structure prediction in the simple-lj force field.

For a demonstration of the random search method, go to subdirectory etanol/simple, and run:
runpa12 s etanol lj (the actual search)
rundist s 10 etanol (clustering)
runpa3 s 13 etanol lj (first energy minimization, Lennard-Jones)
rundist s 20 etanol (second clustering)
rundist s 23 etanol etanol.spf (search for the experimental structure)

This run is done with few trial structures (500 rather than the more usual 10000) to reduce the computer time for this demonstration. It is an exceptionally favourable situation that this setting should take less than one hour, even on an outdated PC, to find the ethanol structure (Z' = 2 in space group Pc).  Note that starting with another random seed is hardly liable to be so successful! Indeed, it is important to verify that, after clustering, low-energy structures are found more than once in the structure generation.
Compare the resulting files with the corresponding ones in etanol/demo. Do not worry about slight discrepancies, the results are sensitive to processor-dependent roundoff errors.
As seen from file distzk.pr2, the experimental structure was found and has Delta-E = 4 kJ/mol in this force field. Note its number and find from files cluslist.pr2, pack.23 or allstruc.pr2 that the rank of that structure is around 75.

Compare with the Buckingham potential:
runpa3 s 23 etanol pw (second energy minimization, Buckingham)
rundist s 30 etanol (further clustering)
There is no improvement at all. Note that searching by
rundist s 33 etanol etanol.spf
would not work: ETANOL.CIF should be optimized in the Buckingham force field, delivering 0001.spf, and then:
rundist s 33 etanol 0001.spf

It is interesting to repeat the entire procedure with the OPLS force field, which is supposed to work relatively well for this compound:
runcod etanol.con.a opls (note different atom labelling)
runmk out.tp opls
Here you shall note a discrepancy between out.up and etanol.tpa. The reason is that the dihedral angle C-C-O-H should not be taken from the carbohydrate version of OPLS [9] as used in UPACK, but rather from the standard parameter set. Continue the structure generation in directory opls:
runpa12 a etanol opls
rundist a 10 etanol
runpa3 a 13 etanol opls
rundist a 20 etanol
rundist a 23 etanol 0001.spf  (of course, after optimization of ETANOL.CIF in OPLS)
As seen from file distzka.pr2, the experimental structure was found and has Delta-E = 1 kJ/mol in this force field. The  ranking is now below 10, an enormous improvement over the simple force field. But, of course, this refers to only one space group. A complete search would give quite disappointing results.

If you wish to study details of, say, the first 10 structures, SPF-type files can be obtained by running:
runpa3 a 23 etanol opls
rundist a 30 etanol
The input file etanol.pa3, made for this run, also codes for calculation of free energies. It can be seen that there is no improvement.

The high-pressure form ETANOL01.CIF crystallizes in P2₁/c with Z"=1.  Interestingly, it should also have been found in the structure generation in Pc with Z"=2 because this is a subgroup of P2₁/c [4]. However, a much longer search is necessary to find this form.  Indeed, as stated earlier the success in finding the other ethanol structure so quickly is a stroke of luck.

Example: benzene

Go to upack/data/benz. The coordinate file cor.001 is already present. Make the all-atom topology file by running:
runcod benz.con.s
runmk out.tp
Compare the resulting file out.up with benz.tps.
A special force field for benzene  has been published by Williams and Starr [10]. Make the all-atom topology file by running:
runcod benz.con.a opls
runmk out.tp benz
Compare the resulting file out.up with benz.tpa. The corresponding force field is benz.ff.

For the experimental structure, go to the subdirectory benz/exp. In the database we find the experimental polymorphs I (BENZEN01.CIF, Pbca, Z=4) and III (BENZEN03.CIF, a high-pressure form in P2₁/c, Z=2). Note how UPACK treats the special positions in the energy minimizations:
runpp s benz BENZEN...  pw
or:
runpp a benz BENZEN... benz
Due to the sixfold symmetry of benzene there are six equivalent possibilities for identifying the atoms, of which you are asked to select one. In this case the choice is irrelevant.
The special positions make the structure prediction less straightforward, as subgroups with the correct vale of Z must be used. For benzene several space groups are possible, for instance P2₁2₁2₁ for benzene I and P2₁ for benzene III [4]. Fortunately, for this simple compound 50 trial structures turn out to be sufficient (rather than the more usual 10000! ). This demonstration shows how 13 space groups can be treated easily with the following scripts (note that no intramolecular nonbonded interactions are taken into account):

mkdirs (create 13 directories, each one containing cor.001)
rundirs1  s benz lj (the search, Lennard-Jones, followed by clustering)
rundirs2  s benz lj (first energy minimization, Lennard-Jones, followed by clustering)
rundirs3  s benz pw (second energy minimization, Buckingham, followed by clustering)
rundirs4  a benz benz (final energy minimization, Williams-Starr, followed by clustering)
rundirs5  a benz benz (produce SPF and CSSR files from the clustered list)
The intramolecular nonbonded energy has been omitted from all calculations.

To perform a clustering over all space groups, create subdirectory benz/clus and run there:
rundist s -23 benz
rundist s -33 benz
rundist a -43 benz
To search for the experimental structure(s), run there:
rundist s -33 benz 0001.spf
rundist a -43 benz 0001.spf
Before each run, the experimental structure 0001.spf must be created in subdirectory benz/exp with the corresponding force field and for both polymorphs separately. Note in allstruc.pr# that, due to the special positions, the experimental forms are found in several space groups.

Collect all resulting structures within 8 kJ/mol and cluster again: do in the main directory:
select.spf 8000 collect
In subdirectory benz/collect:
cat *.spf >> extra.1
rundist a 0 benz (clusters extra.1; no pack files, so number 0 is arbitrary here)
remove.1  (remove equivalent structures; or, after editing, move to another directory)

Additional experiments can be performed in directory benz/collect by adjusting benz.pp, for instance change to 30 kbar:
runppall.spf a benz benz
format.up (shows summary of all results in allres and dihedrals, if present, in alldih)

Sort the file allres according to column 9 (E-tot). Differences with the energies published in 1998 [11] arise because the benzene molecule was kept rigid in that work. Compare the rankings in both force fields. The Williams-Starr force field is excellent, and the experimental structure for form I comes out at rank 1 (rank 4  in simple-pw). For form III the rankings are 5 and 1, respectively.  Comparing the Williams-Starr rankings for form I and form III we find 1 and 5 at normal pressure, to be compared with 5 and 2 at 30 kbar, respectively. So the ranking interchanges, as expected.
 
 
 

Detailed instructions for the preliminary stage

Space group file

The space group file upack/tops/all.spg is read by the subroutine leesr in various programs. It must have the following entries for every space group:

naam
where naam is the space group name

nspgr  is  imp
where nspgr is the number of space group symmetry operations (E excluded), is is the space group number and imp is the importance number for discarding structures in the clustering program dist. Space groups may occur in different settings, with numbers usually increased by multiples of 1000.

For every operation 1... nspgr:
msp(x) tsp(x) msp(y) tsp(y) msp(z) tsp(z)
where an equivalent position x' (in fractional coordinates) is found from:
x' = msp(x).x + tsp(x). Each msp is either +1 or -1 and each tsp is either 0 or 0.5, which means that space group Fdd2 is not implemented.

iref (1..12) = 0/1: disables/enables optimization of parameters:

cell axes:	1: ax	2: by	3: cz
cell angles:	4: bx	5: cx	6: cy
molecular centers:	7: X	8: Y	9: Z
molecular orientations:	10:	11:	12: ;

itest(1...3) = 0/1 disables/enables skipping certain areas in the grid search:

itest(1) = 1: enables skipping of cells with a_x > b_y
itest(2) = 1: enables skipping of cells with b_y > c_z
itest(3) = 1: enables skipping of cells with a_x > c_z

The file upack/tops/codname.cod is needed for the program codes. It defines symbols for various types of atoms and code numbers for terms in the force field. It has to be prepared only once for a certain class of molecules; if desired, several files ending with the extension .cod may be present. It must contain the following blocks of data, separated by one blank line (=), where A1, A2, ... indicate the fixed symbols for certain atom types and X stands for any atom type:

regel:  one line of arbitrary text
=
A1 amass mpol:  atom symbol, atomic mass (a. u.), code for polarizability (800-899).
=
A1 A2 code  for bond stretching (100-199)
=
A1 A2 A3 code  for angle bending (200-299)
 X   A2  X   code also allowed
=
A1 code  for improper dihedral deformation (300-399). This covers all sequences A1  X   X  X .
=
A1 A2 A3 A4 code  for dihedral angle torsion (400-499)
 X   A2 A3  X   code  also allowed
=
Optional: any number of sets of two lines:

A1 A2 A3... (defining one charge group)
codes for corresponding charges (500-599).

Note: in case of repeated entries only the first one is used, except for dihedral angle torsions where all relevant contributions are added. The general "atom" type X  is only used in the absence of the exact atom sequence, to allow the possibility of different periodicities. (This was not so in earlier versions of UPACK.)
 

Molecular connectivity file

The molecular connectivity file subst.con correspoonds to a certain codes file codname.cod. It must contain the following data:
regel: arbitrary text
nrp: A_g, the number of atoms in the group of molecules.
For every sequence number 1...nrp:
panm code mcg:

panm : atom name (freely chosen, but starting with C, H, O...; e.g.C17)
code : atom type code, corresponding with A1 , A2, ... in codname.cod  (e.g.CR)
mcg : atomic charge code as in the force field to be used. A negative mcg denotes the end of a charge group. If only entries +1 or -1 are given, the charge code number is set to 599, corresponding to zero charge.

Molecular topology file

For every substance subst a topology file subst.tp# must be prepared. This can be done manually, but it is more easy to run the program codes:
runcod subst.con codname

This program needs two input files:

in.con (unit 18), linked to subst.con, giving the molecular connectivity (vide infra).
in.cod (unit 19), linked to upack/tops/codname.cod, containing the force field code numbers (vide supra)

out.tp (unit 20), a framework for the required topology file.

Parameters may be missing in the file codname.cod. This situation is indicated by the warning ***** on the screen and in out.tp.
For united-atom force fields it is necessary to fix the chirality in groups like (CH)R₁R₂R₃ where unwanted inversion might occur. This is done by means of improper dihedral angles like C-R₁-R₂-R₃ which should take values of about -35 or +35, and so for every asymmetric atom the choice between corresponding codes must be made.
In aromatic rings the proper dihedrals can have the same function as improper ones. In saturated rings it may be useful to add improper dihedral angles with weak force constants to prevent unwanted ring inversions; for instance, improper dihedrals for the ring atoms (of about +60 or -60) can be used to fix one of the two possible chair forms in six-membered rings.
The program tries to find atoms C1, C2, C3 to define a reference axes system for the first molecule and equivalent ones. Numbers of other reference atoms may be provided. For different molecules (usually crystal water) the first three atoms are taken.
The torsional codes may be followed by a multiplicity, mul. Normally this entry is not used, but some force fields are designed to use only one of the mul possible combinations X-A-B-Y for all atoms X and Y around two given atoms A and B. This can be implemented by setting np negative for that entry in the force field file. UPACK will still calculate all torsions defined in the topology file, but with np divided by mul.

Examples are to be found in the directories data/glucsa, data/benz and data/etanol. It is instructive to compare the files out.up with the given topology files. As seen for ethanol, the all-atom carbohydrate OPLS force field (vide infra) is not adequate, and manual correction was necessary.
 

Molecular geometry files

For the crystal structure prediction molecular geometry files are needed for every conformation to be studied. They are named cor.001, cor.002, ... Allowed are the  spf, mdi, and cssr formats. The program prep may be helpful for the construction of these files (see also under Preliminaries).

The mdi format is compatible with the GROMOS program package. It uses fixed format:

regel (a80): one line of arbitrary text
nfound (i5): the number of following atoms. For every atom 1... nfound:
name xr (1...3) (10x,a5,5x,3f8.5): atom name, x, y, z (Cartesian nanometer coordinates)
box (1...6) (free format): a, b, c (nm), ,,(degrees). This entry is omitted for a free molecule.

TITL regel: one line of arbitrary text
CELL box (1...6): a, b, c (Angstrom), , , (degrees). This entry is omitted for a free molecule.
SPGR spname: space group name (also omitted for a free molecule).
ATOM name xr (1...3): atom name, x, y, z (one line for each atom). These lines have fractional non-orthogonal coordinates if CELL is present, and Cartesian Angstrom coordinates if CELL is absent.
END (optional if there is only one coordinate set present)

The cif format should be able to read the standard crystallographic information file. Please report errors.

The cssr format is useful for interfacing with other program packages, for example TINKER or CERIUS. It is in fixed format:

box (1...3) (38x, 3f8.3): a, b, c (Angstrom)
box (4..6) ispg spname (21x,3f8.3,9x,i3,a20): , , (degrees), space group number, space group name
nfound ispf titel (2i4,1x,a60): number of atoms; ispf= 0: fractional coordinates; ispf= 1: orthogonal Angstrom coordinates; first title
regel (a80): second title. Here used for structure number and energy.
For every atom 1... nfound:
(i4,1x,a4,2x,3(f9.5,1x),8i4,1x,f7.3) atom number, atom name, x, y, z (fractional coordinates), connectivities (max. 8), and optionally charge.

Force field file

As discussed above, the force field is defined by the file upack/tops/ffname.ff. In the ethanol example we encountered simple.ff and opls.ff. As an example for the Buckingham potential, the Williams and Starr [10] force field benz.ff was used for benzene. More details of these and other force fields with fixed charges, now mostly outdated as charges should preferably be calculated for each substance individually, can be found separately elsewhere.

The force field file upack/tops/ffname.ff is read by the subroutine leest in various programs. To allow an easy change of force fields, the relevant information is transmitted through code numbers as defined in the file codname.cod (vide supra). Upon reading a force field file these code numbers are translated into actual numerical values. The force field files must have the following structure:

regel: arbitrary text
ibuck = 0 (or -1, see below): use Lennard-Jones potential
           = 1: use Buckingham potential
cutoff alphac alpha6 hmaxc hmax6 tole tolk where

cutoff: radius for nonbonded interaction (on a charge group basis)
alphac (nm^-1): Ewald-Coulomb damping by alphac * r
alpha6: same for dispersion (r ^-6) terms
hmaxc: Coulomb cutoff radius in reciprocal space (nm^-1)
hmax6: same for dispersion (r ^-6) terms
tole: omit an energy term from the pair list if it contributes less than tole
tolk: same for the list of reciprocal vectors
If alphac= 0 or alpha6 = 0 the corresponding Ewald-summation is skipped.

npp = 0: do not calculate polarization
npp = 1: calculate intermolecular polarization, neglecting the contribution of the induced dipoles to the electric field.
npp > 1: maximum number of cycles in iterative calculation of the complete polarization energy. Then add (optionally):
poltol schaal gamma nnn mmm nte nttt
if you wish to change the default values 0.01 0.1662 0.85 3 32 1 1. See comments in subroutine nonbep for explanation of these symbols. This option is not recommended since it consumes extremely much computing time and because it may lead to catastrophically large polarization energies, for which unconvincing ad-hoc solutions are introduced [12].

sc14c sc14lj multiply Coulomb 1...4 interactions and Lennard-Jones 1...4 interactions, respectively
iddc [epsil]: regulates use of a relative dielectric constant  (default value for epsil is 1.0)
relative dielectric constant  = epsil       if iddc = 0
relative dielectric constant  = epsil * r  if iddc = 1 (note: r in Angstroms!)

Now the force field data are read, in arbitrary order and possibly mixed with blank lines. The programs use the units nanometer for length and kJ/mol for energy. However, in the force field file different units may be used, which are converted by the factors fact and facr. Angles are in degrees, but angle force constants  and are in energy units rad^-2. The following contributions to the force field are implemented:

code (200-299) ct t0: , ₀ in the harmonic angle bending potential:

code (300-399) cq q0:  , ₀ in the harmonic improper dihedral angle potential:

These terms are necessary to define the chirality around united-atom CH groups and to prevent the deformation of rings and planar groups.

code (400-499) cp np [cp np ...]: V_n,n, ... in the cosine torsional potential:

Some force fields are designed to use only one of the possible combinations X-A-B-Y for all atoms X and Y around two given atomsA and B. This can be implemented by setting np negative for that entry in the force field file. UPACK will still calculate all torsions defined in the topology file, but with np divided by the multiplicity (see also the topology file).

code (500-599) cg: q_i (elementary charge units) in the Coulomb interaction:

where  is the relative dielectric constant, and q_i is multiplied internally by the factor 11.787 needed to obtain the energy in kJ/mol. The 1...4 interactions are multiplied by sc14c.

The van der Waals parameters for nonbonded interactions are treated differently since they are based strictly on atom types rather than on individual atoms. The cards start with two atom types A1 A2, followed by numerical data. Entries for the combination of two equal atoms must be present; if other interaction parameters are absent (or are given the value -1), they are found from a combination rule. For the combination of atom types i and j this rule is generally the geometric average:

but sometimes the arithmetic average is used:

Intramolecular 1...2 and 1...3 interactions are always excluded. For the Lennard-Jones interaction it is possible to specify separate parameters for 1...4 interactions. Missing (or negative) parameters for the combination of two equal atoms are taken from the corresponding general parameters. For missing mixed 1...4 parameters the rule is:

Subsequently, all 1...4 interactions are multiplied by | sc14lj|.

The following input formats are used:

These are the parameters for the Lennard-Jones interaction:

The geometric combination rule applies.

If ibuck = -1: A1 A2 E R [E(1...4) R(1...4)]

As above, but now C₁₂ andC₆ are found indirectly:

The geometric combination rule applies for E and the arithmetic rule for R. This precludes the use of Ewald summation for the attractive terms.

These are the parameters for the Buckingham interaction:

There are no special 1...4 interaction terms in this case. The geometric combination rule applies for A and C, and the arithmetic rule for B.

: _i (volume units), the polarizability of atom i in:

Charges on non-atomic sites

It is possible to introduce charges on "virtual atoms", i.e. on sites which are defined in terms of three real atoms a, b, c, as follows:

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

If this file is present, there must also be the file chadd.xyz, containing the charges. For each distinct molecule there must 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

Note: although this feature looked promising, it turned out not to give better results than the usual atomic point charges.
 
 

Detailed instructions for the search stage

All programs have code-directed input. After each code word one or more numerical entries (separated by one or more spaces) may be required or optional. If a certain code word is not specified, the parameters keep their default values which are given between parentheses {...}. The same goes for unspecified optional parameters (indicated by [...]). Code words that are unknown, or preceded by one or more spaces, are disregarded.

A number of code words occurs in several programs. Some are discussed here, and will be mentioned in the individual input descriptions without repeated explanation.

cutoff cutoff [alphac alpha6 hmaxc hmax6 tole tolk]
Overrides the cutoff value given in the force field file, and possibly also the settings for Ewald summation .

epsil epsil
Overrides the dielectric constant value given in the force field file.

inner inner
A surface charge on polar crystals causes a shape-dependent energy term. Normally it is assumed that this term is compensated by external charges; this corresponds to omitting the h=0 term in Ewald summation [13]. If the same situation is to be simulated without Ewald summation, a term 2p² / (3 V) must be added. In UPACK this is done if inner> 0; if alphac > 0 inner is set to zero. In lattice vibrational calculations, the contribution of this term to the second derivatives is omitted unless inner = 2.

pol npp [poltol schaal gamma nnn mmm nte nttt]
Include polarization, possibly overriding the default values given in the force field file.

sd# ncy#sd tol#sd stp#sd nre#sd shm#sd
Specifies steepest descents energy minimization in various stages (#) of a program: the minimization is stopped after ncy#sd cycles or earlier if the energy decreases less than tol#sd. The initial step size is determined by stp#sd. A new pair list is made after nre#sd cycles or if any parameter changes more than shm#sd. The steepest descents minimization is usually followed by:

cg# ncy#cg tol#cg stp#cg dep#cg shm#cg
Specifies conjugate gradients energy minimization in various stages (#) of a program: the minimization is stopped after ncy#cg cycles, or earlier if the root mean square value of the energy derivatives is less than tol#cg and a recalculation with a new pair list does not change the energy more than dep#cg. The initial step size is determined by stp#cg. A restart with a new pair list is made if any parameter changes more than shm#cg.

files ifils ifilm ifild ifilc
Produce output files if these parameters are larger than 0:
ifils = 1: a separate file (format spf) is made for each structure.
ifils = 2: the results are also written to one continuous file (all.spf).
ifils = 3: the results are only written to one continuous file (all.spf).
ifilm = 1: format mdi.
ifild = 1: format cif.
ifild = 2: format fdat.
ifilc = 1: format cssr.
ifilc = 2: format cssr, the charges are written out too.
 

Program pack12

This program generates structures with any number of groups of molecules. It performs the same action as the combined grid search programs pack1 and pack2, but in quite a different way. Rigid molecules are placed with random positions and random orientations in a cell with randomly chosen crystallographic parameters. Dihedral angles may be set to random values too. The resulting structures are optimized with respect to energy, assuming either rigid or fully non-rigid molecules. An arbitrary number of conformations can be studied consecutively.
It is possible to make use of partial knowledge by limiting the search range for the cell axes and cell angles. It is also possible to take fixed starting values for known positions and orientations of molecules, as well as for dihedral angles.

Execute script: runpa12 s/a/... subst ffname

Input files:

in.tp (unit 20), molecular topology, linked to upack/data/subst/subst.tp#
in.spg (unit 11), space group, linked to upack/tops/all.spg
in.ff (unit 12), force field, linked to upack/tops/ffname.ff
in.pa (unit 13), code-directed input, linked to subst.pa12
cor.001, cor.002, cor.003, ... (unit 14), coordinates for various conformations.
mol#.def (unit 48, optional):  a file defining additial charge sites on virtual atoms.
chadd.xyz (unit 48, optional):  a file containing charges to superside the values from the topology file, and possibly also defining the additial charges on virtual atoms.

Output files:

out.pri (unit 16), general results, linked to pack.pr1
pack.10 (unit 10), structure parameters
pack.19 (unit 19), coordinates

The number of distinct groups of molecules (G) is called ngr. This number must certainly not be larger than 4 to avoid excessively long and usually fruitless computations. Each energy calculation is discontinued immediately if any repulsion term exceeds the value ebig. This will occur in the majority of cases. If not, ncysd cycles of steepest descents and ncycg cycles of conjugate gradients energy minimization follow.

The coordinate file(s) must contain the atoms in exactly the same order as the topology file. If ngr is larger than one, the coordinates can either be specified for each group individually or for a number of groups. In the latter case ngr must be a multiple of that number, and the given coordinates are copied as many times as necessary.

If one decides to stop the calculation before the originally desired number of structures is found, a "clean" stop can be forced by creating a file called STOP.

Input file

The input data discussed above can be manipulated by cards in file in.pa.

ngr ngr{1}
The number of groups of molecules, G.

rand nrand{5000}
The search continues until nrand random structures are found (for every conformation cor.001, cor.002, cor.003, ...).

spgr naamsp
The space group name (upper or lower case, slashes may be omitted). This card must always be present, except when the name of the directory is the same as that of the space group.

emax emax1 [emax2] {40 10000}
Structures with energy higher than emax1 (with respect to the lowest energy found at that time) or higher than the absolute value emax2 are not included in the output file.

axes amin amax bmin bmax cmin cmax {0.4 2.5 0.4 2.5 0.4 2.5}
The cell axes are scanned: a_x = amin...amax, b_y = bmin...bmax, c_z =cmin... cmax. This option disables autmax.

autmax autmax {3.7}
For large cells the axes are increased to sqrt(autmax) * (cell volume). This option is disabled if autmax is negative or if axes is set.

angles bxmin bxmax cxmin cxmax cymin cymax {0 1 0 1 0 1}
The cell angles, represented by the projections of the cell axes on the Cartesian axes, are also varied:
b_x=(bxmin...bxmax) a_x, c_x=(cxmin...cxmax) a_x, c_y=(cymincymax) b_y.

grvfix xxfix yyfix zzfix (1...M_gG_i ) {1000...1000}
Fixed starting centers of gravity for each molecule (but random for each value > 999)

eulfix phifix thefix psifix (1...M_gG_i ) {1000...1000}
Fixed starting Euler angles for each molecule (but random for each value > 999)

hydfix hydfix((1...N_z+N_h),1...G_i )) {1000...1000}
The starting value of each dihedral angle is set to the corresponding value of hydfix, except:
if hydfix >  999, random values are assigned to dihedral angles
if hydfix < -999, the dihedral angles from the input geometry are retained

celopt icopt {1}
Cell axes and angles are kept fixed if icopt = 0.

sd ncysd tolsd stpsd nresd shmsd {100   1  0.001 10 1}
cg ncycg tolcg stpcg depcg shmcg {1000 1    0.01 1   0.05}

ebig ebig{2000}

dexp dexp {-100} [volcor {1.10}]
dexp > 0: The experimental density dexp (in gram cm-3) is used to find a_x.
dexp = 0: No such restriction on a_x. However, this is usually very inefficient since much time is wasted on unproductive starting structures.
dexp < 0: The density is estimated from -dexp random structures (for each conformation) before starting the calculations. Its value is found from the smallest volume encountered.
volcor: The smallest volume encountered is multiplied by volcor to facilitate finding new structures.

inner inner {1}

seed iseed {281197}
The integer iseed initializes the random number sequence.

append iseed
If this keyword is present, the nrand random structures are appended to the results of a previous search (pack.10 and pack.19).  It is necessary to redefine seed iseed here to prevent just repeating the previous calculation! If the previous calculation was stopped uncleanly, take care that pack.10 and pack.19 end with the same entry.

coul icoul {1}
If icoul = 0 no Coulomb interactions are calculated.

chadd ichad{0}
ichad = 0: charges from topology file
ichad = 1: charges from file chadd.xyz, but only for real atoms
ichad = 2: charges from file chadd.xyz, also for virtual atoms.

intra intra {1}
intra = 0/1: disable/enable the calculation of intramolecular nonbonded interactions

flex iflex {1}
iflex = 0: rigid molecules, only intermolecular nonbonded interactions. This option implies intra = 0.
iflex = 1: fully flexible molecules

pres pres {0}
The pressure in atmospheres.

print iprint {0}
Minimal output on the screen is given for iprint = 0. A larger value gives more output.

test dtest {0}
If dtest > 0, the analytical energy derivatives are compared with numerically calculated values using parameter shifts dtest (dtest = 0.00001 is usually a good value).

conn icn {0}
If the subroutine leesx fails, it calls the subroutine connect to try and repair the situation. The geometry is compared with the topology in order to assign the correct atom numbers to the coordinates. There may be more solutions; usually they are identical but the program must be guided to select the correct solution number icn. It will ask you to enter that number manually. However, as the script will lead to an error situation, it is necessary to set conn with the desired solution number icn.

cutoff cutoff [alphac alpha6 hmaxc hmax6 tole tolk]
Overrides the default values from the force field. However, the values from simple-lj.ff, opls.ff or unitat3.ff (0.8 3 3 2 2 0.02 0.02) are recommended: Ewald summation with a still fairly small cutoff radius and elimination of quite a lot of terms in the pair list.

epsil epsil
Overrides the dielectric constant given in the force field file.

stc stcin(1...4) {1 1 1 1}
Steepest descents weights for cell axes, cell axis projections, and atomic coordinates, respectively.

Output files

The file pack.10 contains first data concerning space group symmetry, torsional degrees of freedom and atomic coordinates of the free molecule. For explanation of the symbols see the space group and topology files. The following lines occur:

- ispgr icode
Space group number and sequence code. The sequence code is -1 to indicate that the file is generated by program pack12.
- ndihz ndihh: N_z, N_h (= 0). For every dihedral 1...N_z:
- iome jome kome lome omega
- nrp ireftr (1...3): A_g, followed by the three numbers of the reference atoms.
Then for every atom i = 1...nrp:
- i panm iac xr (1...3),
where panm is the atom name corresponding to the topology file, iac is the sequence number of the atom type code, and xr gives the Cartesian coordinates x_i, y_i, z_i in nm for the first geometry file (cor.001).
- cel  (1...12): two lines with starting values for the parameters
- iref (1...12): one line to control the optimization of 12 parameters

After these initial data, each line characterizes a certain structure:

code1: N₁, the sequence number of this structure in the random series (after the density test but before the energy test)
code2: N₂, the sequence number in this file.
This is kept as identification number in all subsequent calculations.
epot: the energy (kJ/mol)
par (1...): a_x, b_y, c_z, b_x, c_x, c_y (pm) unless constrained by symmetry (i.e, the corresponding iref must not be zero).
ntest: sequence number in the random search (before the density test)
en1   vol: the energy value at the initial calculation, and the molecular volume (both multiplied by 1000 to obtain an integer value)
tau (1...): the values of the dihedral angles
1 (the multiplicity, for later reference)

Clustering by dist yields the file pack.13 which has the same structure. Only code1 now gives the sequence number in the clustered file (code2 still refers to N₂) and two new items are added:

zflag: the number of independent molecules (Z'')
imult: the cumulative number of structures that have been clustered together for this entry

The file pack.19 contains Cartesian atomic coordinates for every structure:

- icod2 box (1...6): N₂, a_x, b_y, c_z, b_x, c_x, c_y
- x (1...3A_g): x₁, y₁, z₁, x₂, y₂,z₂, ...

Program pack3

This program performs the energy minimizations in subsequent stages of the crystal structure prediction. Either rigid models or fully flexible ones can be treated.

Execute script: runpa3 s/a/... filenum subst ffname

To illustrate the use of the file number, we take filenum = 13 (the first instance where pack3 is usually called). As an alternative example between {...} we take filenum= 42; this should make the general principle clear.

Input files:

in.tp (unit 20), molecular topology, linked to upack/data/subst/subst.tp#
in.spg (unit 11), space group, linked to upack/tops/all.spg
in.ff (unit 12), force field, linked to upack/tops/ffname.ff
in.pa (unit 13), code-directed input, linked to subst.pa2{5}
in.par (unit 21), structure parameters, linked to pack.13{42}
in.cor (unit 29), coordinates, linked to pack.19{49}
in.add (unit 14), optional set of additional coordinates, linked to file extra
mol#.def (unit 48, optional):  a file defining additial charge sites on virtual atoms.
chadd.xyz (unit 48, optional):  a file containing charges to superside the values from the topology file, and possibly also defining the additial charges on virtual atoms.

Output files:

out.pri (unit 16), general results, linked to pack.pr2{5}
out.par (unit 30), generated structure parameters, linked topack.20{50}
out.cor (unit 39), generated coordinates, linked to pack.29{59}
For every structure (identified with a five-digit number #) there is the possibility for coordinate output in four formats (unit 17), linked to #.spf, #.mdi, #.dat, #.cssr.
Additional possible output from molecular dynamics:
#.en (unit.15), energy snapshots
#.co (unit.17), coordinate snapshots
#.ve (unit.18), velocity snapshots
#.xvo (unit.19), final coordinates and velocities

If one decides to stop the calculation before the originally desired number of structures is found, a "clean" stop can be forced by creating a file called STOP.

The program automatically detects the number of molecular groups from the number of atoms on the coordinate file. It is also possible to read an (additional) hand-made file in.add with cell parameters and atom coordinates.

It is possible to expand the geometry to independent molecules in space group P1. But the energy still refers to one molecular group. The objective is to see whether the symmetry is destroyed when it is no longer enforced; if so, the structure is in fact unstable. In this situation a significant change in energy occurs, and the resulting symmetry-averaged coordinates do not make sense.

The same objective can be reached by performing a primitive MD-run before (or instead of) the energy minimization. Here too one has the choice whether or not to enforce the space group symmetry. It is obviously better not to do so; in fact, to make the simulation less primitive one should consider more than one unit cell. That option is available in the program crysmd.

It is possible to calculate the lattice vibrations and the associated thermodynamical quantities. Here the occurrence of imaginary frequencies is an indication that the crystal structure is not a stable one.

Subsequent runs

To continue the energy minimization, subsequent runs can be done after clustering. For instance:
rundist a 20 subst
creates pack.23, and
runpa3 a 23 subst ffname
uses the parameter file pack.23 , together with input file subst.pa3 and coordinate file pack.29, to produce results in pack.pr3, pack.30 and pack.39. This procedure can be repeated at will up to pack.pr9.

Input file

The options mentioned above can be manipulated by cards in file in.pa:

expand expand{0}
If this code word is present, the coordinates are expanded to P1. Random position disturbances (with maximum expand) can be added to move the system slightly away from a possible saddle point; a reasonable value would be 0.0005.

maxinp maxinp{5000}
The maximum number of input structures to be read from in.par. To select certain structures from this file, set maxinp < 0 and then give on the same line any desired quantity of sequence numbers N₂. If no sequence numbers are given, the program disregards the information from in.par and operates on all entries from in.cor.

emax emax1 [emax2] {30  10000}
The highest allowed energy: emax1 is relative to the first entry on the input file in.par, emax2 is the absolute value.

celopt icopt {1}
Cell axes and angles are kept fixed if icopt  = 0.

sd ncysd tolsd stpsd nresd shmsd {100    0.01       0.0001   5   0.02}
cg ncycg tolcg stpcg depcg shmcg {5000  0.0005  0.005  0.3  0 .02}
These values ensure good convergence, but this will cost a lot of computer time. If the number of structures is large, an intermediate calculation with fewer cycles and a larger tolerance will be more efficient. However, equivalent structures may then not be recognized in the clustering procedure.

extra [ncycsd tolsd ncyccg tolcg]
If a file in.add is present to provide additional structures, new values for the minimization parameters can be given here. The format of this file may be mdi or spf; in the latter case the data for the different molecules must be separated by an END card. If the file title starts with a number that has not been used before, that number identifies the structure. Otherwise a new number is assigned, starting with 90001.

inner inner {2}

coul icoul {1}
If icoul = 0 no Coulomb interactions are calculated.

chadd ichad{0}
ichad = 0: charges from topology file
ichad = 1: charges from file chadd.xyz, but only for real atoms
ichad = 2: charges from file chadd.xyz, also for virtual atoms.

intra intra {1}
intra = 0/1: disable/enable the calculation of intramolecular nonbonded interactions

flex iflex {1}
iflex = 0: rigid molecules, only intermolecular nonbonded interactions. This option implies intra = 0.
iflex = 1: fully flexible molecules

pres pres{0}
The pressure in atmospheres.

iprint {0}
Minimal output on the screen is given for iprint= 0. A larger value gives more output.

test dtest {0}
If dtest > 0, the analytical energy derivatives are compared with numerically calculated values using parameter shifts dtest (dtest = 0.00001 is usually a good value).

cutoff cutoff [alphac alpha6 hmaxc hmax6 tole tolk]
The default values from most Lennard-Jones type force fields are chosen for speed rather than accuracy.  For instance, the values 0.8 3 3 2 2 0.02 0.02 give Ewald summation with a fairly small cutoff radius and elimination of quite a lot of terms in the pair list. This is a reasonable compromise between speed and precision. But note that some calculations may not converge, or even give absurd results. Redo these (using the maxinp -1 ... option to save time) with, for instance, cutoff 1.2 3 3 2 2 0 0. Especially in the last stages of the analysis it is better to have tole =tolk = 0, and to experiment with increasing the other parameters.

epsil epsil
Overrides the dielectric constant value given in the force field file.

pol npp [poltol schaal gamma nnn mmm nte nttt]

files ifils ifilm ifild ifilc {0 0 0 0}

stc stcin(1...3) {1 1 1}
Steepest descents weights for cell axes, cell axis projections, and atomic coordinates, respectively.

If molecular dynamics calculations are desired:

mdrun nstlim {0} [nrun {1}]
Do molecular dynamics. This is done by subroutine runmd, which has been taken from the GROMOS program package and adapted for use in crystal simulations. See ref. [14] for background. Time units are ps.
The card mdrun must give nstlim, the number of MD steps. Optionally the number of runs (nrun) may be specified. The following cards may also be used:

mdcut cutoff[alphac alpha6 hmaxc hmax6 tole tolk]
Only needed if different from the values used for the energy minimization

mdtsc ntt tempi tautp {1 0 0.10}
Controls temperature scaling:
ntt = 0: disable temperature scaling
tempi: if tempi > 0 the initial velocities are taken from a Maxwellian distribution with T = tempi
tautp: relaxation time for temperature scaling

mdtemp temp0 {298}
temp0: reference temperature for velocity scaling

mdpsc ntp comp taup {4 0.00076 5.0}
Controls pressure scaling:
ntp = 0: no pressure scaling
ntp = 1: isotropic pressure scaling
ntp = 2: orthorhombic pressure scaling
ntp = 3: monoclinic pressure scaling
ntp = 4: triclinic pressure scaling
comp: compressibility (Gromos units), essentially only a scaling factor
taup: relaxation time for pressure scaling

mdpres pres0
pres0: reference pressure (atm)
Only needed if different from the pressure used in the energy minimization

mdshak ntc {3} [tol {0.00001}]
ntc = 1: disable the use of bond length constraints by the shake procedure.
tol: the geometric tolerance for that procedure

mddt dt {0.002}
dt: the time step in picoseconds

mdnscm nscm {500}
Stop translation of the center of mass after every nscm steps

mdnsnb nsnb {10}
Make a new pair list after every nsnb steps

mdpri ntpr ntwx ntwv ntwe ntpw {0 0 0 0 2}
Write output to screen after every ntpr steps (the quantity of output being determined by iprint).
If ntwx > 0, write coordinates to #.co (unit 17) after every ntwx steps.
If ntwv > 0, write velocities to #.ve (unit 18) after every ntwv steps.
If ntwe > 0, write various results to #.en(unit 15) after every ntwe steps.
If ntpw = 0, the coordinate and velocity files are unformatted.

If lattice vibration calculations are desired:

latvib [nq ...] {0}
If this keyword is present, lattice vibration calculations are performed. It is advisable to do energy minimization first. The calculations are first done for q = 0, and subsequently for nq vectors in q-space. Possible modes of input:

latvib nq
Calculations are done for nq randomly chosen q-vectors.

latvib n1 n2 n3
Calculations are done for nq = n1 * n2 * n3 q-vectors, filling a regular grid in q-space.

latvib f1 f2 f3 nq
Calculations are done for the vectors defined by f1, f2, f3 divided in nq equal parts.

temp vtemp(1....) {298, undefined ...}
The thermodynamical quantities (and the vibrational amplitudes, if requested) are calculated at temperatures vtemp(1), vtemp(2), vtemp(3), ...

cutlv cutlv tollv {1.5,  0.0}
There is no Ewald summation possible, so one must be able to choose a larger cutoff (cutlv) than in the other calculations. Matrix elements smaller than tollv are disregarded. Note that here the use of charge groups with appreciable net charge may lead to inaccuracies.

Output files

The file pack.20 contains first data concerning space group symmetry, torsional degrees of freedom and atomic coordinates of the free molecule. This is quite similar to the output from pack12, but note that now icode is equal to the sequence number of the call of pack3. After these initial data, each line again characterizes a certain structure:

code1: N₁, the sequence number in this file
code2: N₂, the sequence number from pack12. This is kept as identification number in all subsequent calculations.
epot: the energy (kJ/mol)
par (1...): a_x, b_y, c_z, b_x, c_x,c_y, X, Y, Z, , , , unless constrained by symmetry (i.e, the corresponding iref must not be zero). Lengths in pm, angles in degrees multiplied by 10.
omega (1...N_h+N_z): torsional angles (degrees)
imult: the number of structures that have been clustered together for this entry in previous runs.

If more than one independent molecule is present, only the cell parameters are given (i.e., centers of gravity, Euler angles and torsional angles are not included).

The file pack.29 contains Cartesian atomic coordinates, in exactly the same format as pack.19.

The molecular dynamics output file #.en contains for every selected time step:

ener (1...46):
1...10: energy components: total, kinetic, potential; a_x, b_y, c_z, b_x, c_x,c_y, volume
11...15: energy components: bonds, angles, impropers, torsions, nonbonded
15...20: not yet assigned
21...30: nine pressure components, total pressure
31...40: nine virial components, total virial
41...46: a, b, c (nm), , , (degrees)
 
 

Program dist

This program compares possibly equivalent crystal structures on the basis of a comparison of interatomic distances. It can be used to eliminate or collect equivalent structures ("clustering mode") or to search for the presence of a certain structure (usually the experimental one; "search mode"). Fairly well-refined structures are required: if the geometrical differences become too large, the structural equivalency may be no longer recognized.

Execute script: rundist u/a [-]filenum subst [coord]

Input files:

in.tp (unit 20), molecular topology, linked to upack/data/subst/subst.tp#
in.spg (unit 11), space group, linked to upack/tops/all.spg.
in.dcl (unit 13; optional), code-directed input, linked to:
- clustering mode: subst.cl
- search mode: subst.zk
parameter file(s) (fort.201, fort.202, fort.203, ...): type pack.20.
coordinate file(s) (fort.101, fort.102, fort.103, ...): type pack.29.
For instance, if filenum = 41 then fort.201 is linked to pack.41 and fort.101 is linked to pack.49.
extra.1, extra.2 ,... (unit 14, optional): additional structures in continuous spf-type files.
in.exp (unit 15, optional): an spf-type file containing one reference structure for search mode, linked to the file coord.

If filenum = -41, all corresponding files in directories at the same level are linked. This means that the call should normally be done in a subdirectory of subst (e. g.  subst/clus).

Output files:

out.pri (unit 16), linked to distcl#.pr# or distzk#.pr#.
In clustering mode also:
- new parameter file(s) (fort.301, fort.302, fort.303, ...): type pack.20.
- a combined list of surviving structures: cluslist, linked to cluslist.pr#.
- a list of clusters of all structures: allstruc, linked to  allstruc.pr#.
- for the continuous spf-type file: file(s) remove.# indicating the removed structures and file(s) select.# containing the remained structures. A file remove.# can be directly called to remove the indicated structures, or edited to move them to somewhere else.

Two structures are considered equivalent if the radial distribution curves for each atom in the asymmetric unit are the same. This is even more specific than comparing powder diffraction diagrams, since not only the distances but also the atom types must correspond.This type is normally equal to the atom type, but it can be set more specific than just, say, "hydroxyl O": in glycerol the two end oxygen atoms are equivalent, but the middle one is different.

The structures to be compared are characterized by their space group symmetry, but there is no necessity that this space group is the same as required by crystallographic conventions. For instance, a P2₁/c structure could also have been generated in P2₁ with two quasi-independent molecules or even (with cell doubling) in P1 with four. So the numbers of molecules in the various structures to be compared may differ. The program first determines for every structure the number of crystallographically independent molecules (Z'') before starting the actual comparison of structures. There, of course, only structures with the same Z'' can be really equivalent. The determination of the actual space group of these structures is quite another matter, for which programs like PLATON [5] may be employed.

Structural information

The program tries first to read cell parameters and atomic coordinates from a set of pack-type coordinate files, one for each space group to be considered. They must be assigned to files fort.101, fort.102, ..., and the corresponding pack-type parameter files (fort.201, fort.202, ...) must be present to select the desired structures, and to carry information on space groups and atom types. These files have the standard UPACK format, as described under pack3. In the script file rundist the second parameter gives the number of the file to be considered (e.g. 30 for the file pack.30); if this number is given as -30 all pack.30 files present in directories at the same level (usually a directory for every space group) will also be linked.

After that, the program tries to find additional structures from file(s) extra.#, which contain an arbitrary number of spf-type entries. Each entry must contain at least the cards:

CELL box (1...6): a, b, c(Angstrom), , , (degrees)
and cards with fractional non-orthogonal coordinates:
ATOM panm x (1...3): atom name, x, y, z (one line for each atom). The atoms must be ordered as in the topology file, and it must be possible to form molecules from their coordinates without using symmetry operations.
The information for every structure is concluded by the card:
END

Obviously, at least either this file or one pair of pack-type files must be present. Optional information can be given by a title card containing one line of arbitrary text:
TITL regel
and the energy can be transmitted by a card:
ENER ener
This is used for selection of structures to be compared; in the absence of this card no such selection is done for this entry.

Symmetry information

In each pack-type file the space group number is given on the first line. In the additional file(s) extra.# the space group information is presented for each entry separately, either by a space group name, e.g.:
SPGR P21/m
or by detailing a set of symmetry operations separated by commas:
SYMM -x, y+1/2, -z
SYMM x, -y+1/2, z
where the redundant card
SYMM -x, -y, -z
is optional. If both types of cards are present, the SYMM information is used.

Lines not starting with any of the described code words are skipped. When only symm cards are used to introduce space group information, the file in.spg is not used. In that case the parameter filenum is a dummy argument, and all space groups can be defined. But if the space group of any entry is defined otherwise, the space group file in.spg is needed and only triclinic, monoclinic and orthorhombic systems can be handled.

Distance selection

It may be necessary or desirable to modify the standard information on atom types as given in the files discussed above for one or more of the following reasons:
- Normally only intermolecular distances are considered. So, if the list of atoms contains distinct molecules, each last atom of a molecule must be marked.
- Atoms may have the same type (vide supra).
- The coordinates of some atoms, for instance hydrogen atoms bound to carbon, carry hardly any significant information. Deleting such atoms can lead to a substantial increase in speed, and also eliminate some structure differences (desirable or not!).
The specification of this information is done through an atom card (vide infra) in the file in.dcl.

Mode selection.

In the standard mode of operation ("clustering") the program removes equivalent structures. Output, sorted to energy, is written to new pack files with the number increased by three: for instance, pack.31 is clustered to pack.34.

Alternatively, a file in.exp may be present which defines a reference structure (in the spf format described above). In this case, the program switches to the "search" mode where each structure of the list(s) is only compared with the reference structure. If in.exp is not present, the script looks for it in directory subst/exp.

Details of the calculation can be manipulated by cards in file in.dcl:

cutoff cutoff {0.5}
The cutoff distance (nm) for the radial distribution curves

margin margin {0}
The number of contacts per molecule that are allowed not to be found

sym isymm {1}
isymm = 0: the space group is changed to P1. This overrules the space group information given on the pack-files.
isymm = 1: space group information is used.

dener dener {1.0}
The comparison between two structures is skipped if their energies differ by more than dener. This does not apply to spf-type entries where no energies are available. Also it applies only to the clustering mode.

rvol rvol {0.01}
The comparison between two structures is skipped if the ratio between their molecular volumes deviates from 1 by more than rvol. This applies only to the clustering mode.

tol tol {0.005 in clustering mode, 0.02 in search mode}
Maximum allowable deviation (nm) between two distances that are being compared.

print iprint {0}
Regulates the quantity of output: larger iprint gives more output.

import import(1) import(2) (...) {1, 2, 3, 4, ...}
Determines importance of the lists as read by the program. There is a default list determined by entries in the space group file all.spg, discarding structures in order of increasing space group symmetry. Priorities for the most common space groups are, for instance:
1      2     3      4     5    6     7        8          9          10        11     12      13
P1, P-1, P2₁, Pc, C2, Cc,Pca2₁, Pna2₁,P2₁/c,P2₁2₁2₁, C2/c,Pbcn, Pbca
and giving higher priority to sets with lower value of G. This order of importance can be changed by using the import card.

select isel {1}
Rules details of the clustering:
isel = 0: discard the structure with lowest energy.
isel = 1: discard the structure from the list with the lowest import number (lowest number within one list).
isel = 2: discard the structure from the list with the highest import number (highest number within one list).

intra intri(1) intrj(1) intri(2) intrj(2) (...)
Consider also intramolecular contacts between atom pairs numbered intri and intrj. These pairs should be chosen to discriminate between several possible conformations of flexible molecules; it is a probably superfluous safeguard. In the absence of this card all intramolecular distances are omitted.

atom
This card, if present, must be the last one in the file in.dcl and must be followed by nrp cards (one for each atom) containing:
(-)k: a sequence number (1... nrp); a negative k denotes the last atom of a molecule.
panm: the atom name
type: an integer number to indicate the atom type. These numbers can be arbitrarily assigned between 0 and 99. If type = 0 the corresponding atom is skipped.
In the absence of an atom card, the program takes the types from the topology file in.tp. Lacking that, it searches for the pack-type file fort.201 and sets the type equal to the type found there. If that file is not present either, then all atoms are given a different type. Note that in the last two cases all atoms are considered to belong to one molecule, and then for cases with Z''>1 an atom card is necessary.
 

Program prep

This program calculates and minimizes the energy of one or more single (crystal) structures, thus optimizing their geometries. The input data are obtained from an spf, mdi, cif or cssr type coordinate file. It is possible to expand the geometry to independent molecules in space group P1 or to consider more than one independent molecule in the asymmetric unit. Lattice dynamical calculations are also possible.

The space group is determined by the coordinate file, where a space group name may be specified. If this information is absent, the space group name is taken from the code-directed input file in.pp. In case of conflict the coordinate file has priority. If still no space group name is available, space group P1 is assumed. The symmetry operators in the coordinate file are combined with those of the space group. If no cell data are given, a free molecule is assumed. The number of independent groups of molecules (G_i) is determined from the coordinate file. Neither the names nor the order of the atoms need to corrrespond with the topology, and arbitrary cell translations of any atom are allowed. If hydrogen atoms are missing in positions that can be guessed easily, the program is able to provide these missing coordinates.

Structures with molecules on special positions or with symmetries higher than orthorhombic are treated in a roundabout way. In the former case the symmetry operators which cause atoms to coincide with themselves or with others are eliminated; in the latter case the molecules which are created by operators that cannot be handled by the program are replaced by "independent" ones. Thus, in both cases the structures are treated in a space group of lower symmetry.

Execute script: runpp s/a...  subst coord ffname

Input files:

in.tp (unit 20), molecular topology, linked to upack/data/subst/subst.tp$
in.spg (unit 11), space group, linked to upack/tops/all.spg
in.ff (unit 12), force field, linked to upack/tops/ffname.ff
in.pp (unit 13), code-directed input, linked to subst.pp#
in.cor (unit 29), one or more sets of coordinates, linked to coord
mol#.def (unit 48, optional):  a file defining additial charge sites on virtual atoms.
chadd.xyz (unit 48, optional):  a file containing charges to superside the values from the topology file, and possibly also defining the additial charges on virtual atoms.

Output files:

out.pri (unit 16), general results, linked to prep.pr#
For every structure (identified with a four-digit number #) there is the possibility for coordinate output in four formats (unit 17), linked to #.spf,#.mdi, #.cif, #.cssr.
For the grid search procedure (cf. pack1) the following two files are produced:
out.exp (unit 18), aummary of experimental data, linked to prep.exp.
out.tgt (unit 19), same data, but after energy minimization. For the experimental structure this is a target function to be retrieved "exactly", linked to prep.tg#.

The input data can be manipulated by cards in file in.pp:

spgr naamsp {0}
If this card is absent or if naamsp = 0 the calculation is done for a free molecule. However, this information can be overridden by space group information on the coordinate file, which may be different for each subsequent structure.

inner inner{2}

intra intra{1}
intra = 0/1: disable/enable the calculation of intramolecular nonbonded interactions

flex iflex {1}
iflex = 0: rigid molecules, only intermolecular nonbonded interactions. This option implies intra = 0.
iflex = 1: fully flexible molecules
iflex = 2: rigid molecules with rotatable OH groups, only nonbonded interactions. The option naamsp = 0 is incompatible with iflex = 2.

coul icoul {1}
icoul = 0: no Coulomb interactions are calculated.
icoul = 1: normal calculation of Coulomb interactions.
icoul = 2:  the topology charges are replaced by the charges from the CSSR file.
icoul = 3:  the topology charges are replaced by the charges from the file charges.cssr.
icoul = 4:  the topology charges are replaced by the charges from the file charges.xyz.

chadd ichad {0}
ichad = 0: charges from topology file
ichad = 1: charges from file chadd.xyz, but only for real atoms
ichad = 2: charges from file chadd.xyz, also for virtual atoms.

pres pres {0}
The pressure in atmospheres.

celopt icopt {1}
Cell axes and angles are kept fixed if icopt  = 0.

expand expand {0} [nx ny nz] {1 1 1}
If this code word is present, the coordinates are expanded to P1. Random position disturbances (with maximum expand) can be added to facilitate the minimization leaving a saddle point; a reasonable value might be 0.0005. Optionally, the structure is further expanded to nx, ny, nz unit cells in the a, b, c directions, respectively.

sd ncysd tolsd stpsd nresd shmsd {100 0.01 0.0001 5 0.02}
cg ncycg tolcg stpcg depcg shmcg {5000 0.0005 0.005 0.3 0.02}

print iprint {0} [ipr1 ipr2 ipr3 ipr4 ipr5 {identical to first parameter if unspecified}]
Minimal output on the screen is given for iprint= 0. A larger value gives more output. If desired, this additional output can be specified for special subjects:
   ipr1: reading of the coordinate file
   ipr2: calculation of bonded energies
   ipr3: calculation of nonbonded energies and forces
   ipr4: energy minimzation
   ipr5: lattice vibrations

test dtest {0}
If dtest > 0, the analytical energy derivatives are compared with numerically calculated values using parameter shifts dtest (dtest = 0.00001 is usually a good value).

conn icn {0}
If the subroutine leesx fails, it calls the subroutine connect to try and repair the situation. The geometry is compared with the topology in order to assign the correct atom numbers to the coordinates.  There may be more solutions; usually they are identical but in any case the program must be guided to select the correct solution number icn. It will ask you to enter that number manually. If you know beforehand which one to choose,  you may set conn with the desired solution number icn.

reset
If this card is present, the C-H and O-H distances are reset to values from the topology file.

cutoff cutoff [alphac alpha6 hmaxc hmax6 tole tolk]

epsil epsil
Overrides the dielectric constant value given in the force field file.

pol npp [poltol schaal gamma nnn mmm nte nttt]

roth dihin (1...N_h)
Reset dihedrals for hydroxyl groups to dihin (note: dihedral i corresponds to (N_z+i))

Vn Vn (1...N_h) in kcal/mol {0.3...0.3}
ntal ntal (1...N_h) {3...3}
The latter two cards may be used if iflex= 2 (cf. pack2)

stc stcin (1...6) {1 1 1 1 100 1}
Steepest descents weights for cell axes, cell axis projections, centers of gravity, Euler angles, hydroxyl torsional angles, and atomic coordinates, respectively.

files ifils ifilm ifild ifilc {1 1 1 2}

molc irfat1 irfat2 {0 0}
Transforms atoms from irfat2 (molecule2) upwards to another equivalent position, such as to make the distance between these atoms and irfat1 (in molecule 1) a minimum.

For a free molecule only:

tors dihin (1...N_z+N_h) {0...0}  [Vimp {10000}]
Target values for dihedrals. These values will be kept (almost) constant during energy minimization by means of a temporary improper dihedral potential Vimp centered at dihin (1...N_z+N_h). If the chosen dihedrals cause an internal collision, Vimp will have to be drastically reduced.

scan iscan dihin(iscan) dscan escan
This card is only operative if tors is present. The dihedral iscan is scanned from the given value dihin(iscan) (which supersedes the corresponding entry from tors) in steps of dscan up to and including escan. The other dihedrals wil be optimized during energy minimization.

Variant for a free molecule with unknown coordinates

This is useful to prepare starting geometries for the search procedure.
Execute script: runpp s/a/...  subst none ffname

The first call produces a file cart.out. Copy this file to cart.in and edit. It contains indices of dihedral angles that have to entered manually (degrees). To minimize the number of dihedrals to be assigned, the program omits hydrogen atoms in a first try; if no chain of four atoms can be found, it tries again with hydrogen atoms. If cart.in is present, the same execute script will do a regular calculation. Check the resulting geometry, especially the dihedral angles. Repeat the procedure for every desired conformation, and rename the output file(s) 0000.mdi or 0000.spf to cor.#, the coordinate file(s) to be used in the crystal structure generation.

If lattice vibration calculations are desired:

latvib [nq ...] {0}
If this keyword is present, lattice vibration calculations are performed. To this end double precision is required; it is advisable to do energy minimization first. The calculations are first done for q = 0, and subsequently for nq vectors in q-space. Possible modes of input:

latvib nq
Calculations are done for nq randomly chosen q-vectors.

latvib n1 n2 n3
Calculations are done for nq=n1 * n2 * n3 q-vectors, filling a regular grid in q-space.

latvib f1 f2 f3 nq
Calculations are done for the vectors defined by f1,  f2, f3 divided in nq equal parts.

temp vtemp(1....) {298, undefined ...}
The thermodynamical quantities (and the vibrational amplitudes, if requested) are calculated at temperatures vtemp(1), vtemp(2), vtemp(3), ...

cutlv cutlv tollv {1.5,  0.0}
There is no Ewald summation possible, so one must be able to choose a larger cutoff (cutlv) than in the other calculations. Matrix elements smaller than tollv are disregarded. Note that here the use of charge groups with appreciable net charge may lead to inaccuracies.

celopt icopt dbox(1....) {1 0.01}
if icopt = 0, the cell axes and angles are kept fixed (but the cell contents are optimized following sd and cg)
if icopt = 1, the cell axes and angles are optimized as usual
if icopt = 2, the cell axes and angles are optimized followed by a search for the global minimum in the free energy. This is done by a grid search of 3 points along each axis, separated by dbox(1). For determination of the free energy lowering, this is usually sufficient. For thermal expansion more iterations are needed, and to this end more entries dbox(2...) can be given. In that case the search is iterated once for every entry, around the previous minimum using the corresponding new grid separation. Since the minimum depends on temperature, this is only done for the first temperature vtemp(1).

boxfix box(1...6)
This card starts the grid search for the global free energy minimum directly at the given box parameters a_x, b_y, c_z, b_x, c_x, c_y, again for the first temperature only.
 

Program crysmd

This program performs a molecular dynamics simulation for one structure. The coordinates of that structure are given in file in.xvi which should be in mdi format (the coordinates may be followed by the three velocities). After the simulation the output file out.xvo can be used as input file for a continuation run.

Execute script: runmd s/a/...  subst coord ffname

Input files:

in.tp (unit 20), molecular topology, linked to upack/data/subst/subst.tp#.
in.spg (unit 11), space group, linked to upack/tops/all.spg
in.ff (unit 12), force field, linked to upack/tops/ffname.ff
in.md (unit 13), code-directed input, linked to subst.md
in.xvi (unit 21), coordinates, linked to coord

Output files:

out.pri (unit 16), general results, linked to md.pr#
out.en (unit 15), energy snapshots
out.co (unit 17), coordinate snapshots
out.ve (unit 18), velocity snapshots
out.xvo (unit 19), final coordinates and velocities
out.ave (unit 14), average results in MDI format

The input specifications are given in file in.md, using the same code words as under the option mdrun in program pack3. At least mdrun nstlim must be specified. Additional code words are now:

spgr ispgr {0}
The number (not the name!) of the space group

ngr ngr {1}
The number of independent groups of molecules in the unit cell

expand [nx ny nz] {1 1 1}
If this card is present, the structure is expanded to P1. Optionally, the structure is further expanded to nx, ny, nz unit cells in the a, b, c directions, respectively.

mdskip ig {0}
The number of runs to be skipped before block averages of the remaining ones are taken.

print iprint {0}
 
 

References

Back to beginning