III. Methods for Simulating Large Systems next up previous
Next: IV. Statistical Mechanics as Up: No Title Previous: II. Molecular Mechanical Potentials

III. Methods for Simulating Large Systems

Size of the system limited by storage on computer and speed of programs.
A. Nonbond cutoffs
-The force/energy loop is almost always the rate-limiting step O( tex2html_wrap_inline2128 ).
-The nonbonding terms are the most time-consuming.
To save computational time, apply a spherical cutoff for nonbonding terms:
V(r)=0 for tex2html_wrap_inline2132
Advantages
-a cutoff can reduce the number of neighbors explicitly considered by a factor of approximately tex2html_wrap_inline2134 (where R is radius of entire system)
- tex2html_wrap_inline2138 should be large enough so that the spherical cutoff should be a small perturbation -for LJ atoms, a typical cutoff radius is tex2html_wrap_inline2140 , which is just 1.6% of the well depth -both the van der Waals and electrostatic interactions are significant up to 15 tex2html_wrap_inline2142 or more
Disadvantages
-calculated properties are for truncated potential and may be different from properties of non-truncated potential tex2html_wrap_inline1966 must test for each new system
-sudden cutoff of nonbonded interactions at tex2html_wrap_inline2138 tex2html_wrap_inline1966 discontinuities in energy and forces
1. Shifted and shifted-force potentials
Energy is discontinuous at tex2html_wrap_inline2150 (i.e. when a pair of molecules crosses this boundary the total energy will not be conserved)
Avoid this by shifting the potential function V(r) by an amount tex2html_wrap_inline2152

displaymath2120

- does not affect the forces so does not affect equations of motion
-does affect the total energy, and its contribution varies depending on the particular configuration
Force is discontinuous at tex2html_wrap_inline2150 (i.e. the force is zero for tex2html_wrap_inline2132 ), which can cause numerical instability
Avoid this by using shifted-force potential to make derivative zero at cutoff distance

displaymath2121


2. Switching functions
Multiply V(r) by a switching function tex2html_wrap_inline2158
-smoothly turns off interactions over a range of distances tex2html_wrap_inline2160 to tex2html_wrap_inline2162
- tex2html_wrap_inline2158 decreases from 1 to 0 between tex2html_wrap_inline2160 and tex2html_wrap_inline2162
Example: Brooks et al., J. Comp Chem. 4, 187-217 (1983)

displaymath2122

This brings the potential V and the force tex2html_wrap_inline2170 down to zero in a smooth way between tex2html_wrap_inline2160 and tex2html_wrap_inline2162

displaymath2123


3. Neighbor lists
-Buffer region is created between tex2html_wrap_inline2162 and tex2html_wrap_inline2178
-A neighbor list is kept for all pairs of atoms closer than tex2html_wrap_inline2178
-Neighbor list is not updated every step
-Only those pairs of atoms in neighbor list are considered during calculation of nonbonded interactions (saves time in calculating distances between pairs of atoms to test if closer than tex2html_wrap_inline2162 )
-Neighbor list is updated often enough to ensure that no atom outside the buffer region can become closer than tex2html_wrap_inline2162 before the neighbor list is updated
Example: update the neighbor list whenever any atom moves more than 1/2 the buffer width
-Buffer width and velocities determine time between updates
4. Charge groups and switching atoms
London-van der Waals: often reasonable to truncate at 8-10 tex2html_wrap_inline2186
Electrostatic: charge-charge interactions go as 1/r tex2html_wrap_inline1966 very long-ranged
-Since most molecules are composed of neutral fragments with dipoles and quadrupoles, the leading term in the electrostatic interaction between molecules or parts of molecules is dipole-dipole, which goes as tex2html_wrap_inline2192 (still long-ranged but falls off much faster)
-If cutoffs are applied on atom-by-atom basis, they could split dipoles, which would artificially introduce a large charge-charge interaction;
Avoid splitting dipoles by applying cutoffs over charge groups
-charge group: small group of atoms near each other that have a net charge of zero
-one atom in each group is designated the switching atom, and cutoffs are applied for distances between switching atoms in different groups
-size of charge group must be smaller than cutoff distance;
Example:
For a typical charge group of size 1-3 tex2html_wrap_inline2142 , cutoffs larger than 7 or 8 tex2html_wrap_inline2142 are reasonable
B. Boundaries
Problem: want to simulate bulk liquid but can only include a relatively small number of atoms tex2html_wrap_inline1966 large fraction of molecules would lie on the surface and thus experience different forces than those in bulk liquid
1. Periodic boundary conditions
-Replicate cubic box to form an infinite 3-dimensional lattice
-As a molecule moves in the original box, all its periodic images move in exactly the same way
-If a molecule leaves the central box, an image enters through the opposite face
-Only need to store coordinates of central box
-Assumes that the properties of a small, infinitely periodic system are the same as those of a macroscopic system
tex2html_wrap_inline1966 must check by performing simulation with different box lengths L
Minimum image model:
-each molecule can interact only with the molecule or molecular image closest to it
-must truncate potentials so that cutoff distances are less than L/2 (i.e. tex2html_wrap_inline2206 )
-maximum of N(N-1)/2 pairwise-additive interaction terms
Explicit image model:
-if tex2html_wrap_inline2210 must generate periodic images of molecules in central box (ghost moleculeses) to as great a distance as necessary for cutoff (i.e. include all ghost molecules that are within tex2html_wrap_inline2138 of a real particle)
-ghost molecules shadow their respective particles in central box
-the molecules that are ghost molecules can change during the simulation
2. Stochastic forces at spherical boundary
-System is represented by a sphere, and the boundary potential mimics the potential that a solvent or lattice would produce
-Outer shell ( tex2html_wrap_inline2214 ) of the sphere is designated as a buffer region, and molecules in this region obey a stochastic (Langevin) equation of motion rather than the standard Newtonian equation of motion
-To prevent evaporation either
1. include an additional shell ( tex2html_wrap_inline2216 ) of frozen molecules
2. include an effective boundary potential
C. Long-range forces
Long-range forces: potentials that fall off no faster than tex2html_wrap_inline2218
(i.e. charge-charge ( tex2html_wrap_inline2220 ); charge-dipole ( tex2html_wrap_inline2222 ); and dipole-dipole ( tex2html_wrap_inline2224 ))
Problem:
conditionally convergent sums tex2html_wrap_inline1966 different answer depending on method of grouping terms
Disadvantage of minimum image method
-Cube will be electrically neutral, but periodic structure will be imposed on the liquid (i.e. similary charged ions will occupy positions in opposite corners of the cube)
1. Ewald sums
-Technique for summing interaction between an ion and all its periodic images
-Convert one conditionally convergent sum into two absolutely convergent sums
-Surround the sphere with a conducting material in order to cancel net dipole
Physical description:
-Surround each ion with a Gaussian charge distribution of opposite sign and equal magnitude to screen the interactions so that they are now short-ranged (and the sum of interactions is absolutely convergent)
-Add a cancelling charge distribution so that the overall potential is identical to the original one
-Sum this cancelling distribution in reciprocal space so that it is absolutely convergent
Mathematical description:
-Multiply lattice sum by a convergence function tex2html_wrap_inline2228 to make the sum absolutely convergent
-Add a term equal to the product of the lattice sum and tex2html_wrap_inline2230 to preserve equality
-Fourier transform this second term so that it is also absolutely convergent
General lattice sum:

equation432

where the first summation is over all simple cubic lattice points tex2html_wrap_inline2232 and the second two sums are over the N molecules in the central box. The prime on the first sum indicates that i=j is omitted for tex2html_wrap_inline2236 .
Ewald sum:

equation444


tex2html_wrap_inline2238 is the complementary error function: tex2html_wrap_inline2240
first term: real space sum
second term: reciprocal space sum
last term: cancels the self term that is included in the reciprocal space sum
-Sums in first and second terms are absolutely convergent
- tex2html_wrap_inline2242 is chosen so that both sums converge appropriately
2. Reaction field method
-The particles within a cutoff sphere are treated directly, and the charges outside the cavity are treated as if they form a dielectric continuum that produces a reaction field within the cavity
-Assumes that the interaction from molecules beyond a cutoff distance can be handled in an average way, using macroscopic electrostatics


next up previous
Next: IV. Statistical Mechanics as Up: No Title Previous: II. Molecular Mechanical Potentials

College of Science WWW Development
Tue Nov 26 09:48:20 EST 1996