Summary: CRYSTAL & all
Dear Netters,
Many thanks for your valuable replies. I tried to collect them all and
added several URL sites, so this summary is rather long. I tried to edit,
remove repititions, etc. but still...
-------------------- Original Question --------------------
Recent summary on CRYSTAL & EMBED rises the question: are there any
other programs available for periodical systems (I mean more advanced
treatment then just empirical force field)? Unfortunately, Crystal92 does
not optimize geometry, and Crystal95 is available only for mysterious EEC
Human Capital and Mobility Network members (as Crystal Homepage says).
AMPAC 2.2 had an option for 1D-periodicals, but it disappeared in later
versions instead of expansion into 2,3D-cases. Does anybody know something
else?
------------ My Comments and Relevant Information --------------
Besides Crystal95 and EMBED, Tourino group was distributing two DFT codes:
PW and LAPW(=WIEN)
Recently in Honolulu (Pacifichem'95) N.Tajima, H.Hirano, S.Tsuzuki &
K.Tanabe presented a poster entitled "Totally ab initio prediction of the
structures of CO2 molecular crystal". Authors used their own software.
PROMET by A.Gavezzotti is not on WWW, contact address could be found on URL
http://www.wkap.nl/natopco/NL50MEET.htm
Info on MCI's POLYMORPH:
http://www.msi.com/support/materials/applic/polymorphAnal/c2_PMAnal.html
Ab initio method for geometry optimization of translationally symmetrical
finite clusters:
http://www.chem.joensuu.fi/people/juha_muilu/Research/ab_initio_method.html
Info on ADF_BAND is available on URL
http://iodine.chem.vu.nl/adf.html
Info on CASTEP is available on URL
http://www.tcm.phy.cam.ac.uk/castep/
Info on CAMP is available on URL
http://hackberry.chem.niu.edu/ChemistrySoftware/MolecularDynamicsPrograms/CAMP
another (all-electron) implementation of Car-Pimentallo scheme is PAW:
http://www.zurich.ibm.com/~blo/
Seems that not only plane vawes but any quantum molecular dynamic program
(such as Argus
(http://hackberry.chem.niu.edu/ChemistrySoftware/SemiempiricalPrograms/argus)
or QMDCP (MOTECC91) is capable to optimize crystal structure from the first
principles...
---------------- New Questions ----------------
1. Can MOPAC93/MOPAC7 optimize 3D-periodicals?
2. What semiempirical/ab-initio MD software is available?
3. Does anybody know e-mail of H.Hirano, S.Tsuzuki or K.Tanabe?
------------------- Replies -------------------
From: Xiaoming Zheng <Xiaoming.Zheng { *at * } unsw.EDU.AU>
Why not use the Car-Parrinello approach and their alternatives?
You can look at Review of Modern Physics 64(1992)1045 and reference
therein. The software implementing this approach calls CAMP-Atami by
Japanese computational physics group (cost $25-$50).
From: Andreas Ehlers <ehlers { *at * } iodine.chem.vu.nl>
there is a good program called ADF_BAND that can calculate
periodic structures using the density funtional theory.
From: Ole Swang <oles { *at * } kjemi.uio.no>
The group of Baerends in Amsterdam maintains a program (ADF_BAND) for
calculation of 1d, 2d, and 3d band structures. It uses DFT and a
combination of Slater and numerical basis functions, and gives quite good
results. Unfortunately, it is rather resource-demanding.
Contact Bert de Velde (tevelde { *at * } chem.vu.nl) for further information.
From: Andrey Khavryutchenko <andrey { *at * } compchem.kiev.ua>
MOPAC 6 has 1D peridicity capability. It seems MOPAC 7 got 2D slabs and
3D solids, but I'm not sure.
From: Fransisco Carlos Lavarda <lavarda { *at * } qtp.ufl.edu>
Also for 1D periodic systems you have the Mopac Package, versions 6, 7, and 93.
From: Edoardo Apra' <og845 { *at * } adrian.pnl.gov>
CRYSTAL92 is a Hartree-Fock program for the study of crystalline systems;
to get informations about it, send an e-mail to
CRYSTA92 { *at * } itocsivm.csi.it
------------------------------------------------------------
From: Georg Ostertag <ostertag { *at * } alf.biochem.mpg.de>
I know a group at the university of Ulm. They are developing a
semiempirical method. As I know, they are now programming a facility to
calculate cristals or polymers. I think, they are also interested in
having contact with people using their method or program.
The people who are working on the program are Rainer Mallon (Dipl.
Physiker) and Manfred Doser (student)
E-Mail: Manfred.Doser { *at * } physik.uni-ulm.de
manfred { *at * } model2an.physik.uni-ulm.de
Rainer.Mallon { *at * } physik.uni-ulm.de
From: <v_uk { *at * } biosym.com>
In BIOSYM technologies we have a good ab-initio software for PBC system.
<see
http://www.sgi.com/Products/appsdirectory.dir/DeveloperIXBIOSYM_Technologies_Inc.html
>
DSOLID: DFT soft based on the numerical basis set for PBC systems
Plane_Wave is a first-principles (ab initio) quantum
mechanics soft package that performs accurate calculations
on solid-state systems.
Plane_Wave calculates variational self-consistent solutions
to the density functional theory (DFT) equations for insulating
and semiconducting systems that are subject to 3D periodic
boundary conditions. The method uses a plane-wave expansion for
the wavefunctions and a pseudopotential representation of the ions.
ESOCS is a first-principle (ab-initio) quantum mechanics packages that
performs accurate theoretical calculations on a wide range of 3D
solid-state systems including metals, semiconductors and surfaces. ESOCS
stands for Electronic Structure of Close-packed Solids. It is based on
spin density functional methods and the Atomic Sphere Approximation (ASA).
The distinguishing feature of the code is a spherical wave expansion of
the wave functions, which are centered at each atomic state of the
solids.Because of this construction, the projection of quantities like
densities of state onto atomic sites is particularly easily computed.
Since ESOCS does not use the frozen core approximation, changes in the
core spectra due to core relaxation can be calculated.
From: Lasse Hemmingsen <lasse { *at * } kvlvax.kvl.dk>
Concerning your question about first principle quantum chemical methods on
periodic systems (crystal-like that is), I suggest that you contact Prof.
Karlheinz Schwarz or Peter Blaha (Prof. Schwarz's email adress is
kschwarz { *at * } email.tuwien.ac.at). Their program, WIEN9?, is made for
calculations on periodic systems, and it has given excellent results for
example in calculating electric field gradients.
From: heinz { *at * } olymp.theochem.tuwien.ac.at
Reply to: Karlheinz Schwarz <kschwarz { *at * } email.tuwien.ac.at>
Our WIEN95 does allow a geometry optimization, since force are implemented
into the program. The optimization, however, is presently tested and thus
is not fully automated yet. Please note that we will organize a workshop
for the use of WIEN95 (10-13 April 1996, in Vienna).
For information on the full-potential linearized augmented plane wave
(LAPW) program WIEN95 see Worl Wide Web:
http://www.tuwien.ac.at/theochem/wien95/wien95.html
From: J Perlstein <jperl { *at * } stc.chem.rochester.edu>
If you are interested in predicting the crystal structure of a molecule
knowing only the atom-atom connectivity, there are at least four programs
presently available for this. One called POLYMORPH from BIOSYM/MSI based
on Karfunkels work another called PACK from Chemical Design, Inc. based on
my own work. Both programs use a Monte Carlo simulated annealing algorithm
to find the global and nearby local minima for a molecule in a particular
space group. The former program should be adequate for <up to three
symmetrically independent or chemically different> rigid molecules, the
latter allows for exocyclic conformational flexibility for
molecules containing up to 12 torsion bonds.
There is also a program available from D. Williams at the University of
Kentucky available from QCPE (and I believe also as a commercial package
from Williams) which uses energy minimization techniques and allows for
some flexibility as well as some software from Angelo Gavezzotti at the
University of Milan which uses a cluster approach to packing.
PACK can handle molecules with up to 130 atoms and 12 exocyclic torsion
bonds. It uses either MM2 for non-hydrogen bonded systems or CFF91 for
hydrogen bonded systems with one molecule in the packing unit.
PACK runs as an interface to the CHEMX/CHEMLIB molecular modeling program
on either SGI's or IBM RS/6000. The CHEMX/CHEMLIB program is available from
Chemical Design, Inc. in Mahwah,NJ. Contact Debbie Dunn at Chemical Design
for further details.(Tel: (201)529-2443).
There are several papers which discuss methods for predicting packing
geometries of molecules in solids:
1)Gavezzotti, A. JACS 113, 4622 (1991)
2)JR Holden, Z. Du, HL Ammon, J. Comput. Chem. 14, 422 (1993)
3)HR Karfunkel, B. Rohde, FJJ Leusen, RJ Gdanitz, GJ Rihs,
J.Comput.Chem 14, 1125 (1993)
4) J Perlstein, JACS 116, 11420, 1994
PACK is being used by a number of groups including
Jenekhe for predicting packing of polymer chains
(See Roberts et.al. Chem. Mater. 6, 658 (1994))
and the Whitten group who are interested in monolayer packing
(See Chen et. al. Phys. Chem. 98, 5138 (1994))
From: "D. E. Williams" <dew01 { *at * }
xray5.chem.louisville.edu>
Computer program mpa (molecular packing analysis) predicts the
structures of molecular clusters or crystals by energy minimization,
based on a variety of published force fields or a user-provided force
field. Molecular docking energy minimization and docking structure
prediction is simply treated as a molecular cluster calculation where
molecules in the cluster are not identical.
The intermolecular or nonbonded energy of the molecular assembly
is represented by a pairwise sum over atoms in different molecules.
The program accepts either (exp-6-1) or (n-6-1) nonbonded inter-
atomic potentials, referred to as Buckingham or Lennard-Jones
functions. A torsional potential is accepted for rotations about
selected internal bonds. The program can handle ionic species and
provision is made for net atomic charges or lone pair electron site
charges.
The structural variables considered by the program are the
rotations and translations of several molecules, and selected
internal rotations. Molecules may be related by space group symmetry
operations. For a crystal calculation lattice symmetry is added and
the six (or fewer) parameters of the unit cell may be selected as variables.
Crystal lattice energy calculations can be made in one, two, or three
dimensions with any space group symmetry. There can be more than one
independent molecule in the crystallographic asymmetric unit. Crystal
lattice sums are accurately evaluated using the accelerated convergence
method. Provision is made for summation over electrically neutral
subunits. The crystal unit cell is checked for correctness before
each contact table is generated. This involves tranformation of the
structure to the conventional primitive reduced cell and moving the
molecules as necessary to be inside the defined cell. First and second
derivatives of the molecular cluster or crystal lattice energy are
evaluated analytically (rather than numerically).
The program uses several energy minimization methods, designed
to find both local and global minima. Eigenvalues and eigenvectors
of the second derivative matrix (hessian) are calculated. If all
eigenvalues are positive (positive definite hessian) the Newton-Raphson
method (NR) is used to find shifts toward the energy minimum. If the
hessian is not positive definite off-ridge eigenvector minimization
(OREM) is used. In this method energy is minimized along eigenvector
directions corresponding to the most negative eigenvalues. In a few
special situations steepest descents (SD) minimization is used.
The program includes provisions for an annealing schedule. Annealing
may be useful to carry a structure from a subsidiary energy minimum to
its global minimum. In the annealing process the structure is given
random shifts or thermal kicks which simulate the effects of thermal
motion. The resulting higher energy structure is then relaxed, but
before the energy is completely minimized additional random shifts are
applied as specified by the annealing schedule. These thermal kicks
occur at specified effective temperatures corresponding to a threshold
energy decrement per minimization cycle. When the effective temperature
decreases to zero the structure energy is fully minimized.
Program mpg (molecular packing graphics) reads an output file
from mpa to graphically display predicted molecular cluster or crystal
structures. Mpg differs from the many currently available molecular
graphics display programs in that it is specifically designed to show
salient aspects of the packing structures of clusters and crystals in a
convenient and informative manner.
Limits on number of atoms, etc.
The default maximums are as follows. These limits can be changed
during program installation.
atoms in input list 1000
molecules in asymmetric unit 64
space group operators 24
molecules in contact table 900
contacts in contact table 100000
neutral groups 200
number of variables 192
potential sets 30
internal rotations 10
Selected publications of Donald E. Williams
D. E. Williams, "Molecular Packing Analysis", Acta Crystallographica
1972, A28, 629-635.
D. E. Williams, "Optimally Spaced Rotational Grid Points", Acta
Crystallographica 1973, A29, 408-414.
D. E. Williams and T. L. Starr, "Calculation of the Crystal Structures
of Hydrocarbons by Molecular Packing Analysis", Computers & Chemistry
1977, 1, 173-177.
D. E. Williams and D. J. Houpt, "Fluorine Nonbonded Potential Parameters
Derived from Crystalline Perfluorocarbons", Acta Crystallographica
1986, 286-295.
D. E. Williams, "Alanyl Dipeptide Potential-Derived Net Atomic Charges
and Bond Dipoles and their Variation with Molecular Conformation",
Biopolymers 1990, 29, 1367-1386.
D. E. Williams, "Science Citation Classic: Nonbonded Potential
Parameters Derived from Crystalline Hydrocarbons", Current Contents
1990, 30(11) 14.
D. E. Williams, "Net Atomic Charge and Multipole Models for the Ab
Initio Molecular Electric Potential", Reviews in Computational
Chemistry 1991, 2, 219-312.
D. E. Williams, "OREMWA Prediction of the Structure of Benzene Clusters:
Transition from Subsidiary to Global Energy Minima", Chemical Physics
Letters 1992, 192, 538-543.
D. E. Williams and T. R. Stouch, "Characterization of Force Fields for
Lipid Molecules: Applications to Crystal Structures", Journal of
Computational Chemistry 1993, 14, 1066-1076.
D. E. Williams, "Accelerated Convergence of R(-n) Lattice Sums",
International Tables for Crystallography 1993, Volume B, 374-382.
Y. L. Xiao and D. E. Williams, "Genetic Algorithms for Docking of
Actinomycin D and Deoxyguanosine Molecules with Comparison to the
Crystal Structure of Actinomycin D-Deoxyguanosine Complex", Journal
of Physical Chemistry 1994, 98, 7191-7200.
D. E. Williams, "Ab Initio Molecular Packing Analysis", Acta Crystal-
lographica 1996, Section A, in press.
QUOTATION (for academic use only)
mpa/mpg single computer license and installation . . . . . .$1,995.00
Installation requires telnet access (with temporary password) to
the computer directory where the software is to be installed
Test example calculations.
Several test files for small molecules are available to illustrate
various features of mpa. The first group of tests is for crystals
starting from the observed structure. The second group illustrates
crystal structure prediction, not starting from the observed structure.
The third group is for prediction of molecular clusters and molecular
docking. For each test there is a named input file (.mpa file extension)
and three named output files (.out, .summary, and .restart). The .out file
is the main output from the program, and the .summary and .restart files
are renamed versions of mpa.summary and mpa.restart produced by the test.
The computer time required for the examples will vary widely depending
on the speed and capacity of the computer. The folowing times were obtained
using a 100 MHz SGI Indy computer.
Test Description Time(seconds)
1 carbon_dioxide 4
2 benzene 14
3 benzene_neutral 26
4 benzene_reciprocal 1
5 pyrimidine 48
6 hexafluorobenzene 175
7 triphenylbenzene 176
8 benzene_large_shifts 38
9 benzene_abinit 2670
10 urea_abinit 1595
11 benzene cluster 198
12 benzene global minimum 221
13 actinomycin docking 318
Index to the test examples.
A. Crystal structure predictions starting from the observed structure.
These calculations relax the structure to an energy minimum; they also
test the accuracy of the force field selected.
1. Carbon dioxide
2. Benzene
3. Benzene: summation over neutral units
4. Benzene: reciprocal lattice sum
5. Pyrimidine: lone pair electron sites
6. Hexafluorobenzene: more than one molecule in the asymmetric unit
7. 1,3,5-triphenylbenzene: internal rotation
B. Crystal structure predictions not starting from the observed
structure. These calculations test the predictability of the crystal
structure. Subsidiary energy minima, as well as a global minima, may be
obtained.
8. Benzene, with large shifts
9. Benzene, ab initio prediction
10. Urea, ab initio prediction
C. Molecular cluster predictions.
11. Benzene four molecule cluster
12. Benzene cluster global minimum
13. ActinomycinD-deoxyguanosine molecular docking
1. Carbon dioxide (carbon_dioxide.mpa, carbon_dioxide.out,
carbon_dioxide.summary, and carbon_dioxide.restart)
This is a quick and easy crystal structure calculation to see
if the program is operational. The program gives a date and time stamp
to identify the run. After the title the unmodified input atomic
coordinates are listed followed by copies of the initial rotation matrix,
initial rotation center, and initial translation. The center of mass
is the new center of rotation but is not different from the initial in
this case.
The control parameters are then listed. There are three atoms in
this molecule (NA=3) and four molecules in the unit cell (NS=4). Atomic
coordinates are cartesian (NCEL=0) and no rotation or translation is
specified (NROT=0). Net atomic charges are specified (NCHRG=1) and 40
contact tables are to be generated (NTAB=40). Carbon-hydrogen bonds are
not to be foreshortened (XSHORT=0.0) and the initial energy is not known
(EZERO=0.0). No annealing is specified (NAGITATE=0) and full output is
requested (NPRINT=0).
The specified force field potential energy parameters are then
listed. Note that the user can supplement or customize force fields by
editing the file mpa.pot.
The program determines the number of molecules (rigid bodies, NTHE=1),
the number of neutral subgroups (NGRP=1), and the number of subrotatations
(rotations about bonds, NPSI=0). A check is made for net molecular electric
charge.
Carbon dioxide crystallizes in space group Pa3_, which has 24
symmetry operations. The reference molecule is placed at the origin pointing
along the diagonal threefold axis of the cell; the molecular symmetry
includes the crystallographic threefold and an inversion. When these
6 operations are deleted from the general operations (as shown in the
International Tables for Crystallography), only three screw axis operations
remain, plus the identity operation; the subgroup is equivalent to P212121
with a=b=c. The program notes that all symmetry matrices are diagonal
(NSDIAG=0). When this is the case atom pair index IA will always be
taken equal or greater than IB, and the energy and derivative values doubled
when IA is not equal to IB. If any symmetry matrix has nonzero off-diagonal
elements, NSDIAG is set to one, and all values of IA and IB are considered.
The input lattice constants are printed, along with the lattice
symmetry type and the cell volume. The center of mass of the molecule is
found and output as the new rotation center (unchanged in this case).
The first calculations conditions line requests energy minimization
to end when the cycle decrement is less than 0.00001 (ZEND); contacts are
to be considered up to 14A (SUMLIM). The accelerated convergence constants
are 0.125 for coulombic interaction (CONV1) and 0.15 for dispersion
interaction (CONV6). The reciprocal sum limit is 0.0 (RLIM) and the charge
scale factor is 1.0 (SCALQ).
The second calculations conditions line requests the use of cartesian
coordinates (ICAR=0, rather than principal axes of inertia), summation over
electrically neutral units (IREQ=1), no output of reciprocal lattice terms
(IROUT=0), and final output of contacts less than 0.5A plus the sum of van
der Waals radii (ILIST=1, ADJUST=0.5). The full hessian is to be calculated
and output (NBLOCK=0) and this is not a molecular cluster calculation
(NCLUSTER=0). If the hessian is not positive definite, one eigenvector
corresponding to the most negative eigenvalue is to be used in OREM refine-
ment (NEVNUM=1). Variable shifts are limited to 0.2A or rad.
The "a" lattice constant (only) is selected. No rotations or
translations are selected, since these are fixed by symmetry. Cell constants
"b" and "c" are treated as dependent parameters, set equal
to "a". by IDEP.
The main calculation begins by scanning cells from -4 to +4 from the
reference cell to generate a molecule list (LTSM=362 molecules) and a contact
table (table 1, ITABW=2172 contacts). The initial lattice energy is -27.4392
kJ/mol (EZERO), which is broken down into a coulombic component of -15.9697
(ECOUL), a dispersion component of -39.4472 (EDISP), and a repulsion component
of 27.9777 (EPACKR). The calculated values of the first (gradient) and second
(hessian) derivatives are listed. In this simple case the eigenvalue of the
hessian is just the second derivative itself. Since all eigenvalues are
positive, a Newton-Raphson energy minimization cycle is initiated. The energy
falls to -28.14341 by shifting the "a" lattice constant by 0.2000;
note that
the cubic symmetry is retained as specified by the dependent parameters, and
that the shift was limited by XLIMIT. In cycle 2 the same contact table is
used (to insure convergence, it is important that final refinement be based
on same contact list) and only the first derivatives are calculated.
Recalculation of the hessian is omitted because it is not expected to
change very fast and its calculation is very time-consuming. The lattice
constant decreases by 0.0625, showing the nonlinearity of the minimization
and the usefulness of shift limits. In cycle 3 there is a further small
adjustment of 0.0039 in the lattice constant.
It is now time to generate a new contact table (ITAB=2) with 320
molecules and 1920 contacts, after the new center of mass is computed (in
this case it is unchanged). Note that the initial energy (-28.2466) is
slightly different from the ending energy of the previous cycle (-28.2464).
The difference can be either positive or negative, and is a normal
consequence of a change in the contact table. Smaller summation limits and
larger variable shifts will increase the difference. In cycle 2 the
variable shift becomes smaller than ZEND and the calculation terminates.
The program lists the final energy values and the energy per molecule, and
the molecular rotation and translation (here zero).
A summary listing of short contacts is produced as requested by
ILIST and ADJUST. Contacts are sorted in order of potential type. A
summary of the lattice constant shifts, cell volume, and reduced cell
is given. Of course, the lattice constant shifts reflect the accuracy of
the force field. A list of the final atomic coordinates is given, both in
cartesian and fractions of the unit cell.
2. Benzene (benzene.mpa)
There are four molecules per cell in space group Pbca, with the
molecules on inversion centers. The published structure lists only half
of the atoms; the rest are generated by applying inversion symmetry. An
idealized symmetrical planar molecule with averaged distances was fitted
to the observed structure, and deletion of the inversion operation from
Pbca yields the operators of P212121. Notice that the C-H distance has
been set to a foreshortened 1.027A to be consistent with the WH86 force
field. Alternatively, the C-H distance in the input coordinates could
have been the normal 1.097A and XSHORT set to -0.070. Only rotations are
allowed, since the molecule must be constrained to the inversion center.
The angles of the orthorhombic cell are also not varied. The input
molecule is rotated by a previously determined matrix which gives the best
least squares fit to the observed orientation.
The program considers 79 molecules with 2767 nonbonded contacts
resulting in an initial energy of -51.6259. Since the eigenvalues of the
hessian are all positive Newton-Raphson refinement is selected. In
cycle 2 of the second contact table the energy drops less than
0.0002 and the energy minimization is finished. The final energy is
-52.5376. The angular shift was 2.909 deg degrees from the input
orientation. A list of short nonbonded contacts in the optimized
structure is given. The a, b, and c lattice constants changed by 0.0676,
-0.1134, and 0.3030A. The final cell and molecular volumes are
481.96 and 120.49 cubic A. The optimized atomic coordinates are listed
in both cartesian and fractional cell coordinates.
3. Benzene, summation over neutral units (benzene_neutral.mpa).
Using IREQ=1 and the same summation limit (10A), the program
selects 80 (instead of 79) molecules which overlap a 10A sphere drawn
around any reference molecule atom. There are now 6240 nonbonded
contacts, requiring more time for the calculation. The summation over
neutral units results in a more accurate initial energy of -51.5461.
The final energy is -52.4576, the angular shift is 2.965 degrees, and
the lattice constant shifts are 0.0697, -0.1139, and 0.3065. Setting
NPRINT=1 in the command line reduces the amount of output. Setting
NBLOCK=1 suppresses the printing of the hessian.
4. Benzene, reciprocal lattice sum (benzene_reciprocal.mpa).
The reciprocal lattice sum is evaluated with the data of
test 2. NTAB=0 causes the program to calculate the lattice energy
only. The individual terms in the coulombic and dispersion reciprocal
sums are listed (76 terms). The total reciprocal lattice sums are
0.0047 and -0.0036 to for coulombic and dispersion, respectively. This
type of calculation is useful for selection of optimum values of the
convergence constants conv1 and conv6. Generally, these constants should
be chosen so as to make reciprocal lattice sum contributions neglegible
relative to truncation error in the direct lattice sum. Note that mpa
does not include reciprocal sum contributions to first or second
derivatives; therefore conv1 and conv6 should be kept small.
5. Pyrimidine, with lone pair electron sites (pyrimidine.mpa)
The molecule is in a general position in an orthorhombic lattice.
Since the space group is polar along z this molecular translation is not
selected. The atom list includes lone pair electron sites X1 and X3 0.25 A
from the nitrogens. Only coulombic interaction is considered between the
lone pair sites. Because the charges are large, IREQ is set to 1 to cause
the program to sum over whole asymmetric units. This improves the accuracy
of the coulombic sum. The calculation quickly converges with a molecular
rotation of 1.719 deg and a molecular translation of 0.064 A. The a, b, and
c cell constants change by 0.0007, 0.2556, and -0.0471 A. The energy at the
input observed crystal structure is -56.5453 and the energy at the calculated
structure has been minimized to -57.3145.
6. Hexafluorobenzene: more than one molecule in the asymmetric unit
(hexafluorobenzene.mpa)
Hexafluorobenzene crystallizes in P21/n with six molecules per cell.
Since the space group has fourfold equivalent positions, one molecule is
in a general position and the second molecule is on an inversion center.
The mpa space group is P21 with three molecules in the input atom list
(asymmetric unit). Molecule 1 and molecule 2 are related by inversion
around one center, and molecule 3 sits on a second inversion center. To
retain the observed symmetry, IDEP numbers are assigned. ISEL selects
the monoclinic cell constants, the rotations and translations of molecule
1, and the rotations (only) of molecule 3. Molecule 2 is required by the
IDEP numbers to rotate the same way as molecule 1, and to translate in the
opposite way. This dependency retains the inversion relationship between
molecules 1 and 2. By not allowing translation, molecule 3 is fixed on
the other inversion center.
Note that in the published x-ray structure only half of the atomic
coordinates of molecule 3 are given; the others are easily generated.
Similarly, all of the coordinates of molecule 2 are generated from the
complete list of atomic coordinates for molecule 1. In the example file
an idealized molecule was fitted to the observed structure to give the
initial rotations and translations.
7. 1,3,5-Triphenylbenzene with intramolecular rotations (triphenbenz.mpa)
This molecule is used to illustrate the use of subrotations, or
internal rotations. The molecule is divided into four parts to enable
summation over electrically neutral groups. Each of the three peripheral
benzene rings is initially coplanar with the central ring. The nucleus
central ring has atoms assigned IPSI=0 including three bonded carbons and
three bonded hydrogens. Rotations are selected about the three
phenyl-phenyl bonds. IPSI values are assigned 1, 2, or 3 for atoms in the
peripheral rings. The specified subrotation potential has twofold
periodicity with a maximum conjugation energy of -40.00 kJ/mol for
coplanar phenyls. The observed cell and space group are used.
Note that in this case EZERO values include intramolecular strain
energy. Initially, the repulsion intramolecular energy ESTRNR is
208.2110, leading to EZERO=-24.7450. In table 4, ESTRNR decreases to
90.4128 and EZERO has decreased to -171.2733. Intramolecular coulombic
and dispersion energy is also included in the calculation (ESTRN1 and
ESTRN6). Including the shifts found from table 14, the final predicted
phenyl rotation angles are 318.74, 38.33, and 319.06, in good agreement
with the observed angles. Note one phenyl "propeller blade" is
reversed
in this structure.
This structure can be used to illustrate the characteristics of
summation over neutral groups to achieve greater accuracy in the coulombic
energy summation. With four neutral groups specified, the contact table
contains 38762 entries. If IREQ is set to 1 in the data file, then IGRP
numbers are ignored and the contact table has 52854 entries. Compute
time in the program is approximately proportional to the size of the
contact table. If only an atom-atom SUMLIM is used (NGRP=1 and IREQ=0)
there are only 8809 entries in the contact table. In this case the
calculation is much faster, at the expense of poorer convergence of the
coulombic sum because of summation over a non-neutral domain.
8. Benzene, with large shifts (benzene_large_shifts.mpa).
Instead of the observed RMINP matrix, the program is given
a unit matrix which gives a grossly incorrect molecular orientation
with an initial energy of 123.9. The observed orthorhombic lattice
constants and space group are initially specified. The cell constants
and rotations of the molecule are varied. The initial hessian has
four negative eigenvalues. OREM along eigenvector 1 reduces the energy
to 27.96017, along eigenvector 2 to 10.73014, along eigenvector 3 to
-9.53552, and along eigenvecto 4 to -12.89525.
In contact table 2 the energy calculated from a new contact
table is -13.1912. Three eigenvalues remain negative, and OREM decreases
the energy to -38.66789. Contact table 3 results in two negative
eigenvalues, and the energy decreases to -42.03437. Contact tables 4, 5
and 6 show single negative eigenvalues, and the energy decreases to
-52.05996. A supplementary SD cycle is done because the energy didn't
decrease much with the last table. With contact table 7 the calculation
enters a NR region and quickly refines to a minimum of -52.5121. Of
course, the rotation of the molecule is large (49.843 deg). The lattice
constants agree with the relaxed observed structure. The conventional
reduced primitive cell is given. In the .summary file the lines where
Newton-Raphson refinement was used are appended with NR.
9. Benzene, ab initio molecular packing analysis space group prediction
with Z=4 (benzene_abinit.mpa).
In this test no observed information is used regarding unit cell
constants, molecular position, or symmetry. To produce this test,
first several runs were made with NROT=3, in which random initial
angular orientation of all four molecules is specified. The initial
cell is 10x10x10A and initially the molecules are in face centered
positions. After several trials, the observed structure was predicted
correctly. The Lattman angles from this successful trial were then
inserted into the input file, with change of NROT to 5.
There are 27 variables in this calculation, enumerated as
follows: 6 cell constants, 3 rotations of the first molecule, and
6 rotations and translations of the subsequent molecules. The first
molecule is not translated in order to set the origin of the lattice.
The end of the .out file shows convergence to -209.5558, or
-52.3890 per molecule. The cell angles are all 90 deg, and the cell
edges agree with the observed values. To determine the space group,
look at the final cell coordinate list. Noting that all C and H
atoms are equivalent in benzene, molecule b is related to molecule a
by -x, y+1/2, -z+1/2, which is a screw axis along y. Similar consider-
ations lead to the identification of the observed molecular packing
space group P212121.
Now look at the .summary file. The initial energy for four
molecules is -53.6005. At first, OREM shifts are made in six different
directions. At table 5 the number of negative eigenvalues decreases to
five. There are reduced cell changes at tables 7, 8, 11, 15, and 18.
NR minimization kicks in at table 12, but is lost at table 14. NR
kicks in again at table 15, only to be lost again at table 16. Beginning
at table 17 a continuous NR path leads to the converged energy minimum.
This calculation is further described in reference Williams (1996).
10. Urea, ab initio molecular packing analysis with Z=2 (urea_abinit.mpa).
Since net atomic charges are large in this molecule, summation is
over a neutral domain (IREQ=1). The Biosym force field is selected
(IPOTR=4). As in the benzene example, a number of random trials were made
with NROT=3 to locate minima. The lowest minimum found agrees with the
observed structure with space group P4_21m. Note that this space group
is not observed with high frequency in the Cambridge Structural Database
of organic crystal structures.
The test file uses NROT=5 and inserts successful trial Lattman angles.
The starting model is body-centered in a 8x8x8A cell. There are 15 degrees
of freedom and the initial energy is -16.6118 for two molecules. The
summary file shows reduced cell changes at tables 4, 5, 9, and 14. The
program switches to NR refinement at table 11, but reverts to OREM in
table 13; at table 14 final NR refinement continues to an energy of
-174.8567 for two molecules. The predicted tetragonal structure is in
agreement with the observed structure; it can be displayed with mpg.
This calculation is further described in reference Williams (1996).
11. Benzene 4-cluster subsidiary minimum (b4cluster.mpa)
The starting model is tetrahedral, with the molecules on alternate
corners of an 8A cube. The lattice constants are set to 200 A with the
summation limit set to 20 A. A large shift limit of 1.0A is allowed.
Rotation and translation of molecules 2, 3, and 4 are selected. There
is only one lattice symmetry entry (LTSM=1) and 864 contacts (ITABW,
(12x12)x6). The initial energy is barely negative, -0.5543. Us of up
to 6 negative eigenvalues was specified, leading to EZERO=-1.5952 in
table 2. At table 14 NR refinement begins, temporarily interrupted at
tables 16 and 17. Final energy minimization at table 20 yields
EZERO=-50.8925. This tetramer structure is subsidiary minimum (1b) of
reference Williams (1992).
12. Benzene 4-cluster global minimum (b4anneal.mpa)
The test is the same as test 11, but with annealing added (OREMWA)
to assist in the search for the global energy minimum. The annealing
schedule is started at 200K and reduced by 20K decrements, with initial
maximum random shifts of 0.5A. Since annealing is a Monte Carlo
process, the results of each run will be different. File b4anneal.out
shows a run that successfully located the global minimum described
in reference Williams(1992).
Agitation begins immediately at 200K, leading to a starting
energy of -0.5585. Energy minimization proceeds along the 6 selected
eigenvectors, yielding -1.9849. Since the energy has not decreased by
at least eagitate (5.9861), another agitation takes place at 180K.
This process continues, lowering the annealing temperature, eventually
leading to a T=0K model with an initial energy of -54.2866. This model
quickly refines to the global energy minimum of -55.6281 reported
in reference Williams(1992).
13. ActinomycinD-deoxyguanosine molecular docking (actinomycin.mpa)
The starting model was taken from the crystal, but not including
the waters of solvation. There are 176 atoms in the actinomycinD
molecule and 32 atoms in each of the two deoxyguanosine molecules.
The lattice constants are set to 200x200x200 so that adjacent cells are
too far away to have any contacts included with SUMLIM=25. If SUMLIM
is large enough, all atoms of each of the three molecules will be
included in making up the contact table. In this case it can be
verified that the contact table contains (176)(32)(2) + (32)(32) =
12288 entries. Since this is not a lattice, the convergence constants
are set to zero. In addition to not varying the cell constants,
actinomycinD rotations and translations are not varied. Each
deoxyguanosine has six rotations and translations allowed, for a
total of twelve variables.
The starting energy is -55.2227. The calculation does not start
in a NR region, as there are four negative eigenvalues. In table 2
the number of negative eigenvalues decreases to one, and in table 3
a NR region is reached. Refinement proceeds to a minimum energy of
-129.4025. The two deoxyguanosine molecules have been rotated by
12.258 and 1.955 deg, and translated by 0.875 and 0.182 A. Of course,
this model represents an isolated gas phase cluster and is not expected
to agree exactly with the crystal configuration. In addition to
deviations normally experienced because of approximations in the force
field and neglect of thermal effects, waters of solvation are not
included in this model.
--------------------------- End of Replies --------------------------
Once again, thank you for replies.
Artem
__ _________
/ \ / _ _ \ Artem Masunov - amasunov { *at * } shiva.hunter.cuny.edu
/ \\ \\ \\ \ Chemistry Department, Hunter College
/ /\ \\ \\ \\ \ City University of New York
/ ____ \\ \\ \\ \ 695 Park Avenue, New York, NY 10021
/__/\__/\__\\__\\__\\__\ Tel: (212) 725-0317, Fax: (212) 772-5332
\__\/ \/__//__//__//__/