| Property Prediction |
Energy calculations take into account van der Waals interactions, Coulomb charges, hydrogen bonding, internal rotations, and hydrostatic pressure. The design of Crystal Packer permits greater flexibility in the use of these terms and in their parameterization.
The current version of Crystal Packer is completely independent of the Cerius2 Open Force Field. Crystal Packer maintains its own Dreiding forcefield parameters, which are available for editing through its control panels.
| For information about | See |
|---|---|
| Building crystals | The Cerius2 Builders book |
| Nonbond energy calculations | The Forcefield-Based Simulations book |
| Assigning charges | The Cerius2 Simulation Tools book |
Rigid units
The asymmetric unit of the crystal is divided into fragment-based rigid units. Nonbond (van der Waals, electrostatic, and H-bond) energies are calculated between the rigid units. During energy minimization, the rigid units can be translated and rotated, and the unit cell parameters varied.
Torsional subrotations within rigid units can be defined, thus introducing an internal degree of freedom to the unit and dividing the unit into rigid subunits. Nonbond energies are calculated between subunits, but with 1-2 and 1-3 valence exclusion applied.
minimization
You can perform a default minimization on the current model in three simple steps:
1. Open the Calculate control panel.
2. Click the Initialize Model button.
3. Click the Minimize Energy of Crystal button.
Defaults for all Crystal Packer options and controls have been carefully chosen so that reasonable results can generally be obtained. However, to take advantage of this module's potential and flexibility, you should take time to thoughtfully prepare the model, forcefield parameters, and energy expression and minimization options as described in this chapter.
General methodology
Model initialization
Before an energy calculation or minimization is performed for the first time on a crystal model, the model must be initialized. During initialization, Crystal Packer assigns rigid-unit groups and assesses the dimensionality of the model. For detailed information about model dimensionality, see "The van der Waals term" on page 118.
The rigid units are based on bonding, so that each molecule is a rigid unit. Rigid-unit groups can be edited after initialization.
1. Open the Calculate control panel by selecting the Calculate item on the CRYSTAL PACKER card.
4. Check the Load Up Hydrogen Bonds box to have hydrogen bonds automatically placed into the list box on the Hydrogen Bonds control panel when the model is initialized (see the online help for more control panel information).
5. Click the Rigid Unit Display button if you want each type of rigid unit to be displayed in a different color in the model window when the model is initialized. This option is useful only when there is more than one rigid unit type in the model.
1. Open the Energy Terms control panel by selecting the Energy Options/Energy terms item on the CRYSTAL PACKER card.
van der Waals interaction range
9. Using the Van der Waals control panel, set the interaction range
for the VDW energy calculation:
11. The Ewald Sum Constant, the Real Sum Limit, and the Reciprocal Sum Limit control the range and real/reciprocal space partitioning of the Ewald summation of electrostatic energy.
13. If hydrogen bonds have been searched for, displayed, and loaded during the model initialization step (see "Model initialization " on page 107), examine the H-bonds in the model and in the list to ascertain that they are what you want.
14. If hydrogen bonds were not assigned during the initialization step, assign them now using the Edit Hydrogen Bonding control panel.
15. If necessary, edit the assignment of functional parameters to the hydrogen bonds in the list box. Use the Parameters popup to specify whether the parameters are displayed as A and B coefficients or as well-radius and depth values. By default, all H- bonds are assigned the same parameter set.
16. If you do not want to calculate the van der Waals energy between H-bonded atoms, or if you want to modify the VDW parameters between H-bonded atoms, open the Acid Hydrogens control panel (see the online help) using the Acid H VDW... button at the bottom of the Hydrogen Bonds control panel.
18. All subrotations must be explicitly defined:
For more information about the underlying theory of this method, please see "Energy expression setup" on page 117
Minimizer constraints setup
A Crystal Packer minimization can contain several variable parameters:
, 
, 
).
Cell parameters must be constrained to maintain the lattice type. For example, in a cubic lattice-cell minimization, only the a dimension is allowed to vary. The b and c dimensions are tied to the a length, and the
, , and cell angles are held at 90°.Translational and rotational freedom may need to be constrained for various reasons. If the rigid unit contains only one atom, its rotational degrees of freedom should be fixed. If only one rigid unit is in the cell and no symmetries applied, the unit's translational degrees of freedom should be fixed.
1. Open the Minimization Constraints control panel (see the online help) from the Minimization Options/Constraints item on the CRYSTAL PACKER card.
3. To view and alter allowed rotations and translations for a rigid
unit:
a. Select the rigid unit in the model window.
4. Repeat Step 3 for each rigid unit in the model.
5. If subrotations have been defined (see page 111), they are listed
in the Variable Subrotations list box. Check those subrotations
that you want to vary during the minimization. Those not
checked are held fixed. Fixed subrotations are marked with an
F in the Subrotation Terms control panel (see the online help for
more control panel information).
Minimizer preferences setup
During minimization, cell parameters, molecular segment position, and orientation and subrotation torsion angles may all be varied in the search for energy minima.
Crystal Packer provides a choice of two minimization algorithms:
Modified Newton
One strategy is to construct a related positive-definite matrix through Cholesky factorization. This is the modified Newton method.
Another strategy is to use the steepest descents method, in which coordinates are displaced in a direction opposite to the gradient of the potential energy.
1. Crystal Packer calculates the energy, as well as the first and second
derivatives of the energy.
2. The derivatives are used to determine the best search direction
in the crystal-parameter space.
3. The minimum energy position in this search direction is determined.
Five user-specified maximum-increment values exist: cell edge, cell angle, molecule translation, molecule rotation, and subrotation torsion. Setting a maximum limit on movement allowed in one iteration promotes the investigation of local minima. This is because these limits restrict the minimizer to exploring only a small area of the crystal parameter space at a time.
Because atom-energy potentials have infinite range, but calculated energies are only between atom pairs within the cutoff distance, the crystal energy can vary discontinuously as the crystal space parameters change and atoms move in and out of the interacting radius. This sort of discontinuity can make minimization searches very difficult. Crystal Packer combats this problem by maintaining a nonbond list of interacting atom pairs that is updated only when necessary. This nonbond list, called the contact table, is used for the calculation of VDW, Coulomb, and H-bond energies.
Instead of updating the contact table at regular step intervals as some modules do, Crystal Packer updates its nonbond list when a given maximum interatomic displacement is reached. A running sum of atom displacements is maintained, and the nonbond list is refreshed when any interatom distance has changed by more than the threshold displacement since the last refresh.
The minimization ends when the maximum number of iterations has been performed or when the overall energy gradient falls to below the specified gradient tolerance value. (Energy gradients with respect to each degree of freedom are first calculated, then these components are combined in quadrature to produce the overall gradient.)
Associated with the gradient tolerance is the energy tolerance. The energy-tolerance check is performed after the gradient-tolerance condition has been met. The energy-tolerance check is necessary if the nonbond list is not updated after each iteration. If the nonbond list is out of date when the minimization converges, it could be that the convergence condition would not have been met if the nonbond list were current.
To set minimization preferences
1. Open the Minimizer Preferences control panel (see the online
help) from the Minimization Options/Preferences item on the
CRYSTAL PACKER card.
3. If necessary, edit the value in the Non-bond Update Displacement
Threshold entry box.
5. If necessary, edit the maximum number of minimization iterations.
Calculating energy and running the packing calculation
The energy calculation and energy minimization functions are the core of Crystal Packer. All procedures discussed until now have been done to prepare the model, the energy expression, and the minimizer for the energy and packing calculations.
Output from both calculations is sent to the text window. Additionally, for the minimization calculation, a plot of energy (kcal/cell) versus minimization step is sent to the graph window.
1. Open the Calculate control panel from the Calculate item on
the CRYSTAL PACKER card.
To perform a minimization of model energy
Before performing a minimization calculation on the model, the initialization procedure (see page 107) must be carried out. Also, unless you are using default settings for all options, you should first perform the energy terms setup procedure (see page 108), the minimization constraints procedure (see page 112), and the minimization preferences procedure (see page 115). To perform a minimization:
1. Open the Calculate control panel from the Calculate item on
the CRYSTAL PACKER card.
The energy expression set up by Crystal Packer can involve up to five terms:
Each energy term is discussed below.
VDW interactions apply to all atom pairs that are:
Polymer, surface, and
Nonetheless, some exceptions to the rule prohibiting interunit bonding must be made; that is, for polymer, surface, and 3D network models. Since polymers in crystals are covalently bonded throughout the crystal in the c-axis direction, Crystal Packer ignores interactions between a polymer segment and translational images of that segment in the c direction. Similarly, a 2D surface network is assumed to be bonded in the b-c plane, and interactions with images in cells (0, m, n) are excluded. A 3D network is bonded throughout the crystal, and no self-segment interactions are included. During the initialization step, Cerius2 automatically sets the rigid unit dimensionality according to the dimensionality of the model.
Crystal Packer provides an option to use convergence accelerating smearing (Pertsin and Kitaigorodsky 1987) to estimate the VDW energy of interactions between atom pairs further apart than the cutoff distance. This method works by assuming that, outside the truncation radius, atoms are uniformly smeared over space. This assumption leads to an easily calculated expression for estimation of long-range VDW interactions.
During the energy calculation, Cerius2 also checks for close contact atoms. For all atoms in different rigid units that are separated by less than the specified close contact distance, a warning message is sent to the text window. This warning has no effect on the calculation.
The VDW energy of atom pairs is calculated using the Lennard-Jones functional form (also called the 12-6 form):
Eq. 1
Where: 
Eq. 2
Where:
By default, Dij0 and Rij0 are calculated by combining values for the two elements involved in the interaction. The well depth is always combined geometrically, but the well radius can be either the arithmetic or geometric mean of the homo-atom radii:
Eq. 3
and
Eq. 4
or
(arithmetic mean)
Eq. 5
Where:
(geometric mean)
Crystal Packer also allows you to specify off-diagonal parameters for unique pairs of atoms. These pair potentials can be entered as well depth and radii (Dij0 and Rij0) or as R-12 and R-6 coefficients (A and B). Values defined for specific pairs take precedence over values determined by elemental combination.
The Coulomb term
Electrostatic interaction energy can be included in the Crystal Packer energy expression. Partial charges can be assigned or calculated using the Charges module (see the Cerius2 Simulation Tools book). The overall charge of the crystal must be neutral.
The electrostatic interaction is calculated between all charged atoms residing in different rigid units (or different subunits, if subrotations have been defined). To prevent summing over small, essentially random charges, you may want to set a minimum charge value, so that any smaller charges are excluded from the Coulomb sum.
The Ewald summation method (Ewald 1921, Karasawa and Goddard 1989) is used to calculate the Coulombic energy. The slowly converging real-space Coulomb sum is split into a quickly converging modified real-space sum and a summation in reciprocal-space.
The choice of these three parameters can be quite tricky. The larger the two limits, the more accurate the Coulombic energy, but the more expensive the calculation. Because the reciprocal-space sum is cheaper than the real-space sum, it is better to err on the side of a large sum constant, small real-space sum limit, and large reciprocal-space sum limit. However, the default values of 0.40, 10, and 0.5 for the Ewald sum constant, the real-space sum limit, and the reciprocal-space sum limit, respectively, are reasonable choices that the novice user may never need to alter.
| The Ewald summation technique converges only if the unit cell is electrically neutral. Therefore, Crystal Packer does not permit Coulombic energy calculations on cells with a non-zero net charge. |
The hydrogen bond term
Hydrogen-bond energy can be included in the Crystal Packer energy expression. Crystal Packer can be set to automatically find all H-bonds in a model according to specified bond search criteria, or you can assign H-bonds to the model yourself.
H-bond potential
Either one of two 12-10 functional forms can be used to calculated the H-bond energy:
A H D = Angle formed by the acceptor, hydrogen, and donor atoms
Parameters, expressed as R0 and D0, or as A and B, are supplied in the Hydrogen Bonds control panel; their source is the Dreiding II forcefield (Mayo et al. 1990). Some of these parameters are designed for use when Coulomb energies are included and others are for use when electrostatic interactions are ignored.
In some forcefields, H-bonds are parameterized assuming no explicit van der Waals interaction between the hydrogen atom and the acceptor atom. Crystal Packer lets you choose to exclude H-Acceptor van der Waals interactions from the van der Waals energy calculation.
The torsional energy term
During minimization, molecules in a crystal are given translational and rotational freedom. Additionally, internal degrees of freedom may be added by allowing torsional rotation within the molecules. Each rotatable torsion, called a subrotation, must be explicitly defined before the model is initialized.
Eq. 10
Where:
t = Torsion angle (as shown in figure above).0 = Reference torsion angle. t) = Torsion energy as a function of torsion angle.t, C1, C2, C3, and C6) from the Dreiding forcefield are supplied with Crystal Packer, but these are easily edited to suit your particular application (Mayo et al. 1990). When a torsion is initially defined, it is arbitrarily assigned function one of the Dreiding set, regardless of the chemical nature of the chosen subrotation. If function one is not suitable to you, assign a different parameter set from the list or enter a new set of your own.
The external pressure term
The hydrostatic energy term contribution to the energy (Ehyd) is:
Eq. 11
Where:
Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem., 94, 8897 (1990).
Clark, M.; Cramer, R. D.; Van Opdenbosch, N. J. Comp. Chem., 10, 982 (1989).
Ewald, P. P. Ann Phys (Leipzig), 64, 253 (1921).
Karasawa, N.; Goddard, W. A. J.Phys. Chem, 93, 7320 (1989).
Gill, P. E.; Murray W.; Wright, M. H. Practical Optimization,