3752
J. Phys. Chem. 1995, 99, 3752-3764
Crystal Packing without Symmetry Constraints. 1. Tests of a New Algorithm for Determining Crystal Structures by Energy Minimization K. D. Gibson and H.A. Scheraga" Baker Laboratory of Chemistry, Cornell University, Ithaca, New York 14853-1301 Received: July 18, 1994; In Final Form: November 16, I994@
We describe an efficient algorithm for determining crystal structures of rigid molecules by energy minimization, in which no constraint is imposed except for the existence of a fully variable lattice. The algorithm makes use of advances in secant methods for minimization, and in techniques for computing analytical derivatives of empirical potentials, to achieve both speed and accuracy of convergence. Tests of the algorithm with crystals of 25 organic molecules, in which the number of independent variables varied from 12 to 114 and energy minimization was started from an experimentally observed crystal structure, showed that (1) the final energies and lattice parameters were independent of the starting lattice (centered or primitive), (2) the final space group symmetry, determined by mathematical superposition of the molecules in the unit cell, was identical to the initial space group symmetry within close tolerances, and (3) in most cases, some or all symmetry was lost during energy minimization and regained before convergence. Varying the limits at which energetic contributions from nonbonded interactions were cut off, from values adopted frequently in the literature out to larger distances, had little effect on the final lattice parameters but could result in changes in total energy of as much as 10%. Tests with three different potentials, in which the pairwise interatomic contributions to the energy had the same mathematical form but different parameterization, indicated that, for full accuracy in reproducing crystal structures, changes might have to be made in the mathematical form as well as the parametrization of these potentials.
I. Introduction Crystal packing computations involving energy minimization have been used extensively for testing and refining empirical and semiempirical potential energy functions for organic Commonly, the energy is considered to be a sum of pairwise interatomic interactions and a mathematical expression derived from physical considerations, in as simple a form as possible, is assigned to the energy of each interatomic i n t e r a c t i ~ n . ' The ~ ~ ~adjustable ~ parameters in these expressions are chosen so that the lattice parameters and lattice energies of a set of organic crystals are reproduced sati~factorily!~~'-~~ The potential energy function is then ready for use in theoretical studies of these and similar compounds in the solid or liquid state, or in conformational studies of macromolecules composed of atoms of the same types. In almost all studies of this kind, the computational algorithm treats the molecules as rigid bodies and maintains the space group symmetry throughout the energy minimization. This frequently reduces the number of independent variables and the time required for the computation by a large factor. Two widely used programs that follow this general protocol are Busing's and the packing programs of Williams.25 Both programs allow the user to specify the space group symmetry, and both make some provision for internal motions within the molecules; however, in the great majority of applications the molecules have been regarded as rigid bodies, and the experimentally observed group symmetry has been imposed and maintained throughout the energy minimization. A further similarity between the two programs is in the mathematical methods used to search for local minima, which are not nearly as efficient as methods that have been developed during the past three decades.26 One point of difference is that Williams' programs use analytical derivatives of the energy, whereas in W I N derivatives are estimated numerically. @
Abstract published in Advance ACS Abstracts, February 15, 1995.
A different approach was taken by Lifson and Hagler and their collaborators. 1~12,27-29 These workers treated the molecules in the unit cell as independent rigid bodies and imposed no restrictions other than the existence of a variable lattice. Derivatives were calculated analytically, and an early, but effective, secant method30 was used for energy minimization. In a study that covered a number of crystals of organic acids and amides, inspection showed that in almost all cases, the final structure had the same space group symmetry as the starting structure; however, these studies did not address the question whether the symmetry was maintained throughout the energy minimization. Hagler and co-workers concluded that maintenance of symmetry is of little or no use as a criterion for assessing the accuracy of empirical potential^.^^,^^ Symmetry constraints have been relaxed in a few other studies, with results that were variously interpreted by their authors. Zielinski et al. minimized the energy of CBr4 with symmetry constraints and then removed all constraints except for the existence of a lattice.'* Although they considered the final structure to be triclinic, inspection of their Table I11 shows rather small deviations from the original C2/c symmetry. Dzyabchenko minimized the energy of benzene, using a unit cell with the number of molecules in the unit cell (Z)equal to 4 and no imposition of symmetry, from various starting structures, and found in almost all cases that the final structure had acquired a recognizable space group symmetry.31 In studies of crystal packing by molecular dynamics, symmetry has also been noted in the averaged equilibrium structure^.^^-^^ All these studies suggest that retention or adoption of space group symmetry is likely to occur in crystal packing computations as a consequence of the existence of a lattice, provided that the potential energy function is a reasonable one. Thus, maintenance of the observed space group symmetry, after energy minimization starting from an experimental crystal structure, should be considered a necessary criterion in the evaluation of any empirical potential energy. In this paper, we describe tests
0022-3654/95/2Q99-3752$09.QQIQ 0 1995 American Chemical Society
Crystal Packing without Symmetry Constraints
J. Phys. Chem., Vol. 99, No. 11, 1995 3753
of this statement that have been carried out with an efficient new algorithm for minimizing the energy of crystal structures composed of rigid molecules, with no requirement other than the existence of a lattice. This algorithm takes advantage of two recent advances in methodology to achieve speed and efficiency: (1) the dramatic advances in secant methods for minimization that have occurred over the past three decades,26 and (2) the observation that analytical gradients of the energy with respect to generalized coordinates can be computed at almost no cost if the Cartesian gradient of the energy is accumulated first.35 In the accompanying paper,36we use the algorithm to explore the effect of varying the starting structure on the space group symmetry of the final structure, using the crystal of benzene as a model. 11. Methods
A. Definitions and Conventions. We denote the spacefixed unit vectors along Cartesian axes by el, e2, e3, the basis vectors of the unit cell by 81, 8 2 , 8 3 , and the basis vectors of the reciprocal unit cell by bl, b2, b3. The parameters of the lattice are denoted by a, b, c, a, p, y. The basis vectors are chosen so that the direction of a1 coincides with el, the vector a2 lies in the (el, e,) plane, and the lattice vectors form a righthanded system; this choice of axes and vectors, in place of others that are often sed,^',^^ simplifies some of the algebra. The basis vectors are related to the Cartesian unit vectors by
L and U are related by
The positions and attitudes of the Z molecules in the unit cell are specified by translation vectors m t and Euler angles &, e,, q,,, (1 5 m IZ). The Euler angles are defined using a common convention, in which 4 is the angle between el and the line of nodes, 9 is the angle between e3 and the corresponding body-fixed axis k, and q is the angle between the line of nodes and the body-fixed axis i. During energy minimization, the origin and axes el, e2, e3 of the Cartesian coordinate system are held fixed; all translation vectors t, and Euler angles & Om, ly, and the six lattices parameters a, b, c, a, p, y , are allowed to vary independently. B. Energy Function. The energy is assumed to be a sum of pairwise interatomic interactions, in which each pairwise interaction is a function solely of the interatomic distance. Three types of interaction are considered (1) Coulombic, (2) nonbonded, and (3) hydrogen bond. In the work presented here, these terms take the respective forms:
a1 = L11e1
= L31e1 + &ZeZ L33e3 where Lij are components of the lower triangular matrix 3'
I" b
=
0
0
1
cos y b sin y 0 cos a - cos p cos y ) c c o s p c( sin y
(2)
In eq 2,
W = (1 - cos2a - cos'p
- cos2 y
+
2 cos a cos j3 cos y)lR (3) The matrix L is a lower Cholesky factoSg of the metric tensor G for the lattice:
LL~=G
(4)
The reciprocal basis vectors satisfy
b, =
In eq 9, qi and q, are partial charges on the atoms, multiplied by an appropriate factor so as to express the energy in kilocalories/mole (1 cal = 4.18635 J). In eq 10, EV is the energy (in kcal/mol) of the Lennard-Jones 6-12 potential at the energy minimum for the pair of atoms i and j , and ri is the interatomic distance at the energy minimum. The parameters E; and
(14)
the condition is satisfied. Let the lattice vector from the unit cell containing the origin to an image unit cell be
1 = nlal
where
+ n2a2+ n3a3
To derive a similar inequality for the remainder Rz from the sum in eq 20, let r, = maxij IrijJ, where i a n d j represent atoms in the unit cell. Then
(15)
and denote by n the vector ( n l , n2, n3)Tof integers. Then the distance d between corresponding points in the two cells is
d = (IT1)’” = (nTGn)’I2= (nTLLTn)’”
(16)
The summation over the lattice is carried out over all image cells that satisfy d Ie, i.e. for all cells that satisfy
In,l ln21 5 lnll 5
+
R2 5 wJ:yV
5 e/L3,
2 2112-L {(e2- ~,,n,)
Replacing summation outside the sphere of radius 4 by integration, setting c = 4 - r,, and using the fact that y r, 2yify 2 4
(17)
32 3 1 / 4 2
[{e2- ~i,n,2- (h2n2+ L32n3)2}1/2
erfc(ky)
erfc(ky) dy 2n1” +c exp(-pc2) k
- ~~~n~L3in3]/Lii
2. Coulombic Contributions. It is well known that the sum of terms of the form in eq 9 over a lattice converges very slowly, and it is necessary to transform the sum to improve the convergence. We use Ewald’s transformation in a form similar to that used by Catti?
The expressions in eqs 23 and 26 are both set equal to E and solved by Newton’s method to give values for h, and 4, for a given value of k. To determine a suitable value for k, use has been made of Catti’s estimates of the numbers n~ and n~ of lattice points or reciprocal lattice points within a The adjustable constant k is changed in increments of 0.025 from 0.55 to 0.15, eqs 23 and 26 are solved to give values of h, and 4, and these are substituted in Catti’s formulae to give estimates of n ~ ( 4 ) and nR(h,). The value of k for which ( 2 n ~ n ~is)smallest is chosen at the start of energy minimization and subsequently kept constant. The factor 2 arises because each term in eq 20 takes at least twice as long, in CPU time, to compute as does each term in eq 19. In the studies presented here, was determined from the formula
+
where and
rij
is the vector from atom i to a t o m j in the unit cell,
2k
F3(rij)= erfc(klrijl)/lrijl - -dij
n
In eq 19, V is the volume of the unit cell, and the sum is taken over all nonzero reciprocal lattice vectors of the form
h = m,b,
+ m,b2 + m3b3
(22)
with ml, m2, and m3 nonnegative. The constant term C in eq 18 represents the “self-energy” of the molecules in the unit cell, i.e. the sum of terms of the form in eq 9 taken over the atoms within the same molecule. The sums in e q s 19 and 20 include all terms within spheres of radii h, and 4, respectively, where h, and 4 are chosen so as to ensure that the remainder in each infinite series is less than a sufficiently small number E . If R1 is the remainder from the sum F1 in eq 19, Catti40 has shown that R, 5 Q $ e r f c E )
This value was sufficiently small to ensure that no detectable part of the electrostatic contribution to the energy was lost. If the unit cell has a nonzero dipole moment, it is necessary to add an extra term J to the electrostatic energy to allow for the fact that the Ewald series converges conditionally but not absolutely. We use the expression for an infinite sphere derived by Smith4]
where M is the dipole moment of the unit cell. For crystals of organic molecules, this term is sometimes a highly significant correction; for example, for the crystal of 1,3,5-trithiane, J accounts for approximately 50% of the total electrostatic energy (our unpublished observations). D. Gradient of the Energy. Partial derivatives of the energy with respect to the independent variables and the lattice parameters can be computed rapidly by following a procedure in which the Cartesian gradient is calculated first.35 Let ry) denote the vector of Cartesian coordinates of the ith atom in the mth molecule in the unit cell. We first accumulate the vector sums
J. Phys. Chem., Vol. 99, No. 11, 1995 3755
Crystal Packing without Symmetry Constraints
In eq 29, CP denotes a sum of functions of the forms in eqs 10, 11, 20, and 21, and Lis a nonzero lattice vector of the type shown in eq 15. The first sum is over all atoms in all molecules in the unit cell other than the mth molecule, and the second sum is over all molecules in the unit cell. An analogous set of vectors, arising from the part of the Ewald formula that involves the reciprocal lattice, is also computed:
The Cartesian gradient, V/!m)E,is the sum of cy’ and vj(m); the components of the gradient with respect to the rigid body variables are calculated from the Cartesian gradient as described elsewhere!2 To obtain the gradient with respect to the lattice parameters, we f i s t compute the three vector sums
( k = 1,2,3) (31) In eqs 34, nl, 112, and n3 are the (integer) coordinates of the lattice vector L(see eq 15). These sums are accumulated only over terms that involve the direct lattice, i.e. terms of the forms in eqs 10, 11, and 21. The contributions from these vector sums to the gradient with respect to the lattice parameters are
= -sl(2) cos y - sf) sin y
db
-dIdir a
sin a - s(3)c sin a(cos a - COS p COS y = S’(3)c siny W sin y
(3) $3
-I&
dY
c sin B(cos ,!?- cos a cos y )
W sin y
= s\”b sin y - si2’bcos y s(3) COS ,6 2
s(3)c[cos 3
y(cos2 a
+
- cos a cos y ) -
sin’ y cos2 - cos a cos ~ ( 1 cos2 y ) ]
+
w sin’ y (32)
Since sf), si’), and sb’) do not appear in eqs 32, they need not be computed. The contributions from terms involving the
reciprocal lattice can be expressed through three analogous vectors dk); details are given in Appendix A. It should be noted that the sums in eqs 29 and 31 assume nothing about the form of the functions CP(1ry) - $)‘ - 4). This procedure can be used to calculate the gradient of any pairwise interatomic potential depending only on interatomic distances. A final contribution to the gradient comes from the dipole correction term J in eq 28. This necessitates the addition of a term equal to (4d3V)qy’M to each cy’, and a term arising from the volume dependence of J to the gradient with respect to any lattice parameter. E. Starting Points for Energy Minimization. Since the original purpose of the present work was to test parametrizations of potential energy functions, observed crystal structures were used as starting points for energy minimization. Some effects of varying the starting points are discussed in the accompanying paper.36 To generate initial values for the rigid body parameters of the molecules in the asymmetric unit, the fractional coordinates of the atoms in these molecules, as obtained from the literature, were converted to Cartesian coordinates referred to the space-fixed axes el, e2, e3. Body-fixed coordinates were then defined by setting their origins at the centroids of the molecules and aligning the body-fixed axes with the principal axes of inertia (with unit weight allotted to each atom). Initial values for the translation vectors drn)were set equal to the Cartesian coordinates of the centroids; initial values for the Euler angles were determined from the orientations of the new bodyfixed axes. The remaining molecules in the unit cell were generated by applying the symmetry operations listed in the International Tables for Cry~tallography~~ for the experimentally observed space group to the fractional coordinates of the molecules in the asymmetric unit; these new sets of fractional coordinates were used to generate translation vectors and Euler angles as described in the preceding paragraph. To facilitate comparison of the relative movement of the molecules in the unit cell, the origin of coordinates was reset to the centroid of all the molecules in the unit cell before and after energy minimization. F. Determination of Symmetry Relations between Molecules. Rather than attempting to characterize the space group of the final structure, we choose to compare the symmetry operations that relate the molecules in the unit cell during or after energy minimization with those that relate them at the start. If these operations are unaltered, and if the lattice parameters whose values are fixed by the lattice class have not changed significantly, the space group type can be considered to have been conserved or restored during energy minimization. To determine the symmetry operations relating pairs of molecules, advantage was taken of the fact that the general space group operation and the operation of superposing two molecules are both affine transformations. Superposition of molecule p on to molecule m results in an orthogonal matrix S and a shift vectors referred to the origin of Cartesian coordinates; together, these determine an affine transformation T = {S, s}, which may or may not be a crystallographic symmetry operation. Associated with a crystallographic symmetry operation is the screw vector t, defined by t =s
+ s s + ... sk-’s
(33)
where k is the order of S. For T to be a crystallographic symmetry operation, the following criteria must be satisfied: (i) the trace of S must be an integer; (ii) if S is transformed from the Cartesian frame el, e2, e3 into the lattice frame 81, 8 2 , 8 3 , all matrix elements must become integers; (iii) if t is
Gibson and Scheraga
3756 J. Phys. Chem., Vol. 99, No. 11, 1995
TABLE 1: Crystal Structures compound acetic acid acetamide N-acety lglutamine N-acety lglycine adipic acid benzene benzene n-butane 1,4-diisopropylbenzene 1,4-di-tert-butylbenzene 4,6-dimethyl-l,3-dithiane ethane formic acid glutaramide glutarimide n-heptane malonamide N-methylurea naphthalene n-octane propionic acid sulfur ( s 6 ) sulfur ( S d 1,3,5-trithiane 2,4,6-trimethyl-1,3,5-trithiane p-xylene a
spacegroup Pna21 R3c P212121 P21k P21Ic
Pbca P21Ic P21Ic
P2Jn P2dn a/c P21h Pna21 a1c P21lc
P1 P21Ic P212121 P21la P1 P2llC
R3 Pnnm
Pmn21
cc
P2Jn
Z 4 18 4 4 2 4 2 2 2 2 8 2 4 4 4 2 8 4 2 1 4 3 2 2 4 2
a 13.255 11.513 13.811 4.589 10.01 7.39 5.417 4.110 5.584 6.3 12 20.08 4.226 10.241 6.19 10.226 4.15 13.07 8.4767 8.235 4.22 4.07 10.818 4.725 7.668 9.768 5.806
b
3.963 5.095 11.546 5.15 9.42 5.376 7.621 8.237 9.983 5.71 8.623 3.544 8.26 7.416 19.97 9.45 6.981 6.003 4.796 9.04 9.104 7.003 12.306 5.032
lattice parametersa C a 5.762 12.883 12.914 14.633 10.06 6.81 7.352 8.097 11.842 9.906 15.69 5.845 5.356 17.46 7.30 91.3 4.69 8.04 6.922 8.659 94.0 11.02 10.97 4.280 14.532 5.285 8.521 11.215
B
Y
138.29 136.75 110.0 118.603 102.11 96.02 118.4 130.83 102.28 74.3 73.0
85.1
ref 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
122.55 84.0 90.56
116.33 100.48
105.8
72 73 74 75 76 77 78 79
a, b, c in A; a, b, y in deg. Only the parameters not fixed by the space group are listed.
transformed into this same frame, its elements must all become integers if the lattice is primitive or if S is not the identity mat1ix.4~9~Transformation into the lattice frame changes S and t to S' and t':
S' = USLT t'=Ut
(34)
The three tests are applied sequentially to the transformation
T obtained by superposing a molecule in the unit cell on to a molecule in the asymmetric unit. If the trace of S is found to be an integer, the type and order of S can be determined from the values of its trace and determinant (ref 43, p 795). The screw vector t can then be calculated. If the transformed matrix S' and vector t' have integral components, further tests are applied to determine the position and direction of an axis of rotation, screw rotation, or rotoinversion, or the normal to a mirror plane. These tests are described in Appendix B. In the present work, each molecule that was not in the original asymmetric unit was superposed on to the molecule from which it was originally generated, and the tests described above were applied. To allow for roundoff errors and inexact convergence, all quantities that should be integers were allowed to deviate from exact integral values by f0.03. If a molecule was no longer related crystallographically to its original progenitor, a search was made for another molecule to which it might be related crystallographically. Only if no such molecule could be found was it considered to have become part of the asymmetric unit during energy minimization. G. Crystals and Potentials. The crystals used in the present work are listed in Table 1. The potential function used in most of the work was an early version of a potential that is currently under development in this laboratory. Partial charges were calculated separately for each molecule as d e ~ c r i b e d ? ~The ,~ parameters for nonbonded and hydrogen-bond interactions, which should not be regarded in any way as final, are listed in Table 2. For comparison with other potentials, parameters were
TABLE 2: Nonbonded and Hydrogen-Bond Parameters atom type 6 (kcallmol)" p (A,. H (-Csp3) 0.04164 2.81 H (-amide) 0.05 100 2.43 H (-Csp2) 0.01595 3.13 H (-hydroxyl) 0.37416 1.85 csp3 0.05059 4.02 0.10618 3.61 c (4) C (other sp2) 0.10797 3.912 N (amide) 0.06379 3.89 N (amine) 0.05214 4.06 0.21739 3.07 0 (carbonyl) 0 (-O-) 0.39220 3.315 S 3.66 0.80500 (kcallmol)b J (A)* 3.103 1.912 1.480 1.880 (O-)H* * e(+) 4.734 1.733 (O-)H.* O(
