Discontinuous molecular dynamics for semiflexible and rigid bodies

A general framework for performing event-driven simulations of systems with semiflexible or rigid bodies interacting under impulsive forces is outlined. The method consists of specifying a means of computing the free evolution of constrained motion, evaluating the times at which interactions occur, and determining the consequences of interactions on subsequent motion. Algorithms for computing the times of interaction events and carrying out efficient event-driven simulations are discussed. The semiflexible case and the rigid case differ qualitatively in that the free motion of a rigid body can be computed analytically and need not be integrated numerically.     

Isotropic-nematic spinodals of rigid long thin rodlike colloids by event-driven Brownian dynamics simulations

The isotropic-nematic spinodals of solutions of rigid spherocylindrical colloids with various shape anisotropies L/D in a wide range from 10 to 60 are investigated by means of Brownian dynamics simulations. To make these simulations feasible, we developed a new event-driven algorithm that takes the excluded volume interactions between particles into account as instantaneous collisions, but neglects the hydrodynamic interactions. This algorithm is applied to dense systems of highly elongated rods and proves to be efficient. The calculated isotropic-nematic spinodals lie between the previously established binodals in the phase diagram and extrapolate for infinitely long rods to Onsager's [Ann. N. Y. Acad. Sci. 51, 627 (1949)] theoretical predictions. Moreover, we investigate the shear induced shifts of the spinodals, qualitatively confirming the theoretical prediction of the critical shear rate at which the two spinodals merge and the isotropic-nematic phase transition ceases to exist.     

Periodic orientational motions of rigid liquid-crystalline polymers in shear flow

The collective periodic motions of liquid-crystalline polymers in a nematic phase in shear flow have, for the first time, been simulated at the particle level by Brownian dynamics simulations. A wide range of parameter space has been scanned by varying the aspect ratio L/D between 10 and 60 at three different scaled volume fractions Lphi/D and an extensive series of shear rates. The influence of the start configuration of the box on the final motion has also been studied. Depending on these parameters, the motion of the director is either characterized as tumbling, kayaking, log-rolling, wagging, or flow-aligning. The periods of kayaking and wagging motions are given by T=4.2(Lphi/D)gamma(-1) for high aspect ratios. Our simulation results are in agreement with theoretical predictions and recent shear experiments on fd viruses in solution. These calculations of elongated rigid rods have become feasible with a newly developed event-driven Brownian dynamics algorithm.     

Brownian dynamics simulations of the self- and collective rotational diffusion coefficients of rigid long thin rods

Recently a microscopic theory for the dynamics of suspensions of long thin rigid rods was presented, confirming and expanding the well-known theory by Doi and Edwards [The Theory of Polymer Dynamics (Clarendon, Oxford, 1986)] and Kuzuu [J. Phys. Soc. Jpn. 52, 3486 (1983)]. Here this theory is put to the test by comparing it against computer simulations. A Brownian dynamics simulation program was developed to follow the dynamics of the rods, with a length over a diameter ratio of 60, on the Smoluchowski time scale. The model accounts for excluded volume interactions between rods, but neglects hydrodynamic interactions. The self-rotational diffusion coefficients D(r)(phi) of the rods were calculated by standard methods and by a new, more efficient method based on calculating average restoring torques. Collective decay of orientational order was calculated by means of equilibrium and nonequilibrium simulations. Our results show that, for the currently accessible volume fractions, the decay times in both cases are virtually identical. Moreover, the observed decay of diffusion coefficients with volume fraction is much quicker than predicted by the theory, which is attributed to an oversimplification of dynamic correlations in the theory.     

Coarse-grained rigid blob model for soft matter simulations

We have developed a coarse-grained multiscale molecular simulation method for soft matter systems that directly incorporates stereochemical information. We divide the material into disjoint groups of atoms or particles that move as separate rigid bodies; we call these groups "rigid blobs," hence the name coarse-grained rigid blob model. The method is enabled by the construction of transferable interblob potentials that approximate the net intermolecular interactions, as obtained from ab initio electronic structure calculations, other all-atom empirical potentials, experimental data, or any combination of the above. We utilize a multipolar expansion to obtain the interblob potential-energy functions. The series, which contains controllable approximations that allow us to estimate the errors, approaches the original intermolecular potential as the number of terms increases. Using a novel numerical algorithm, we can calculate the interblob potentials very efficiently in terms of a few interaction moment tensors. This reduces the labor well beyond what is required in standard molecular-dynamics calculations and allows large-scale simulations for temporal scales commensurate with characteristic times of nano- and mesoscale systems. A detailed derivation of the formulas is presented, followed by illustrative applications to several systems showing that the method can effectively capture realistic microscopic details and can easily extend to large-scale simulations.     

Molecular alignment of rigid rods in nonrigid spherical pores

We have investigated the orientation ordering of two shish-kebab chains confined by spherically harmonic potentials through Monte Carlo simulations and asymptotic analysis. The rigid rod is modeled as shish-kebab chains consisting of tangent hard spheres aligned in the same axis, and the harmonic potential is chosen to model nonrigid cavities. We first show that the interactions between a rod and the spherically harmonic potential are independent of chain orientation, indicating that the alignment of two confined rods arises from the excluded volume interactions alone. In the strong fields, the order parameter of two confined rods converges to different values, depending on the parity of chain length. From asymptotic order parameters, we find that the rods of odd-number beads rotate more freely even under the limiting strong confinement. However, the two rods of even-number beads are essentially trapped in a configuration of perpendicular alignment through intercalation of their central grooves. We attribute the dependence of the parity of chain length to the different locations of the center-of-mass in a rod for these two cases. Furthermore, we compare the shish-kebab chains with different rod models in the simulations, and utilize these models to explore the effect of the local rod smoothness on molecular alignment. Our findings suggest that increasing local rod smoothness enhances the rotational degree of freedom for confined rods, and the effect of local rod roughness emerges under strong enough applied potentials.     

Performance of a rigid rod statistical mechanical treatment to predict monolayer ordering: a study of chain interactions and comparison with molecular dynamics simulation

A statistical mechanical model that treats hydrocarbon self-assembled monolayer (SAM) chains as rigid rods is examined to interrogate the mechanisms involved in monolayer ordering. The statistical mechanical predictions are compared to fully atomistic molecular dynamics simulations of SAMs with different packing densities. The monolayer chain order is examined as a function of surface coverage, chain-surface interactions, and chain–chain interactions. Reasonable interaction potentials are deduced from ab initio electronic structure calculations of small model systems. It is found that the chain-surface interaction is the most important parameter in formation of flat-lying monolayer phases, while formation of standing phase monolayers is driven most importantly by increased density of molecules at the surface. A brief discussion of the utility and validity of the rigid rod treatment is given in light of the molecular dynamics results.     

Multiple grid methods for classical molecular dynamics

Presented in the context of classical molecular mechanics and dynamics are multilevel summation methods for the fast calculation of energies/forces for pairwise interactions, which are based on the hierarchical interpolation of interaction potentials on multiple grids. The concepts and details underlying multigrid interpolation are described. For integration of molecular dynamics the use of different time steps for different interactions allows longer time steps for many of the interactions, and this can be combined with multiple grids in space. Comparison is made to the fast multipole method, and evidence is presented suggesting that for molecular simulations multigrid methods may be superior to the fast multipole method and other tree methods.     

Discontinuous molecular dynamics for rigid bodies: applications

Event-driven molecular dynamics simulations are carried out on two rigid-body systems which differ in the symmetry of their molecular mass distributions. First, simulations of methane in which the molecules interact via discontinuous potentials are compared with simulations in which the molecules interact through standard continuous Lennard-Jones potentials. It is shown that under similar conditions of temperature and pressure, the rigid discontinuous molecular dynamics method reproduces the essential dynamical and structural features found in continuous-potential simulations at both gas and liquid densities. Moreover, the discontinuous molecular dynamics approach is demonstrated to be between 3 and 100 times more efficient than the standard molecular dynamics method depending on the specific conditions of the simulation. The rigid discontinuous molecular dynamics method is also applied to a discontinuous-potential model of a liquid composed of rigid benzene molecules, and equilibrium and dynamical properties are shown to be in qualitative agreement with more detailed continuous-potential models of benzene. The few qualitative differences in the angular dynamics of the two models are related to the relatively crude treatment of variations in the discontinuous repulsive interactions as one benzene molecule rotates by another.     

Collective motion artifacts arising in long‐duration molecular dynamics simulations

We tested a variety of molecular dynamics simulation strategies in long‐duration (up to several nanoseconds) constant‐temperature simulations of liquid water under periodic boundary conditions. Such long durations are necessary to achieve adequate conformational sampling in simulations of membrane assemblies and other large biomolecular systems. Under a variety of circumstances, serious artifacts arise in the form of spurious collective behavior that becomes obvious only after the simulation has gone at least several hundred picoseconds. The potential energy of the system drops and the system changes from a liquid to an icy or glassy state. The underlying cause is accumulated center‐of‐mass motion of the system, coupled with velocity rescaling associated with constant‐temperature control. The velocity rescaling in the constant‐temperature algorithm reduces the thermal velocity as the net center‐of‐mass velocity grows, effectively causing the kinetic energy of the system to drain from thermal motions into coordinated motions. We found that the incidence and magnitude of the underlying artifactual motion leading to the spurious transition is mediated by: choice of method for computing electrostatic interactions; choice of ensemble; size of the simulation cell; SHAKE tolerance; frequency of nonbonded pairlist updating; and closeness of coupling to the temperature bath. The appearance of the spurious transition can be avoided by periodically subtracting net center‐of‐mass motion during the dynamics, or by improving the accuracy of the simulation by means of tightening SHAKE tolerance and updating nonbonded pairlists every timestep. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 121–131, 2000     

Depletion potentials in colloidal mixtures of hard spheres and rods

The depletion potential between a hard sphere and a planar hard wall, or two hard spheres, imposed by suspended rigid spherocylindrical rods is computed by the acceptance ratio method through the application of Monte Carlo simulation. The accurate results and ideal-gas approximation results of the depletion potential are determined with the acceptance ratio method in our simulations. For comparison, the depletion potentials are also studied by using both the density functional theory and Derjaguin approximations. The density profile as a function of positions and orientations of rods, used in the density functional theory, is calculated by Monte Carlo simulation. The potential obtained by the acceptance ratio method is in good agreement with that of density functional theory under the ideal-gas approximation. The comparison between our results and those of other theories suggests that the acceptance ratio method is the only efficient method used to compute the depletion potential induced by nonspherical colloids with the volume fraction beyond the ideal-gas approximation.     

Computer simulations of molecular motion in liquid crystals by the method of Brownian dynamics

Computer simulations of the molecular motion of polymer chains in the presence of a strong nematic field were carried out by the method of Brownian dynamics. Two models were studied: the first model (linear liquid crystal) is a freely jointed chain with rigid bonds, the second model (comb-like liquid crystal) is a chain with fixed bond angles and rigid side groups. The influence of ordering on chain conformations, orientational and translational mobility and spectra of relaxation times was investigated.     

Stress tensor in model polymer systems with periodic boundaries

The calculation of the stress tensor from molecular simulations of atomistic model polymer systems employing periodic boundary conditions is discussed. Starting from the dynamical equations governing the motion of sites, correct double summation forms of the atomic and the molecular virial equations are derived, which are valid for flexible, infinitely stiff and rigid chain models even in the presence of interactions between different images of the same parent macromolecule. A new expression for the true instantaneous stress (flux of momentum through the faces of the simulation box) is derived and shown to exhibit large fluctuations when applied in molecular dynamics simulations. A new equation for the thermodynamic stress, cast exclusively in terms of intermolecular forces on interaction sites, is also derived. Application to Monte Carlo simulations shows that the molecular virial expression exhibits the smallest fluctuations among all stress expressions discussed, and thus allows computation of the thermodynamic stress with least uncertainty. A scheme is developed for the calculation of surface tension from intermolecular forces only.     

Dynamics of attractive vesicles in shear flow

A nonequilibrium molecular dynamics (NEMD) method is employed to study the dynamics of two identical vesicles with attractive interactions immersed in shear flow. The dynamics behaviors of attractive vesicles depend on the attractive interactions and the shear rates simultaneously. There are four motion types for attractive vesicles in shear flow: a coupled-tumbling (CTB) motion, a coupled-trembling (CTR) motion, a collision/rotation mixture (CRM) motion and a separated-tank-treading (STT) motion, which are determined by the competition between the shear flow and the attractive interactions. Furthermore, the dynamics behavior of an individual vesicle shows three main motion types such as tumbling, trembling and tank-treading motions, and relies mainly on the shear rates. Meanwhile, comparisons with rigid vesicles for the dynamics behaviors are made, and the collision/rotation mixture (M) motion isn’t observed for rigid vesicles.     

A polarizable ion model for the structure of molten AgI

The results are reported of the molecular dynamics simulations of the coherent static structure factor of molten AgI at 923 K using a polarizable ion model. This model is based on a rigid ion potential, to which the many body interactions due to the anions induced polarization are added. The calculated structure factor is in better agreement with recent neutron diffraction data than that obtained by using simple rigid ion pair potentials. The Voronoi-Delaunay method has been applied to study the relationship between voids in the spatial distribution of cations and the prepeak of the structure factor.     

An anchor-dependent molecular docking process for docking small flexible molecules into rigid protein receptors

A molecular docking method designated as ADDock, anchor- dependent molecular docking process for docking small flexible molecules into rigid protein receptors, is presented in this article. ADDock makes the bond connection lists for atoms based on anchors chosen for building molecular structures for docking small flexible molecules or ligands into rigid active sites of protein receptors. ADDock employs an extended version of piecewise linear potential for scoring the docked structures. Since no translational motion for small molecules is implemented during the docking process, ADDock searches the best docking result by systematically changing the anchors chosen, which are usually the single-edge connected nodes or terminal hydrogen atoms of ligands. ADDock takes intact ligand structures generated during the docking process for computing the docked scores; therefore, no energy minimization is required in the evaluation phase of docking. The docking accuracy by ADDock for 92 receptor-ligand complexes docked is 91.3%. All these complexes have been docked by other groups using other docking methods. The receptor-ligand steric interaction energies computed by ADDock for some sets of active and inactive compounds selected and docked into the same receptor active sites are apparently separated. These results show that based on the steric interaction energies computed between the docked structures and receptor active sites, ADDock is able to separate active from inactive compounds for both being docked into the same receptor.     

Sorting of brownian rods by the use of an asymmetric potential

We present here a method for sorting nanometer scale brownian rods by using a switching asymmetric periodic potential. A two stage sorting process is used to isolate particles with specific dimensions, with acceptable sorting times as well as realizable potential barrier lengths. The method was tested using computer simulations. The ability to sort the nanometer scale anisotropic particles, such as gold nanorods, portends important applications in large scale data recording, photothermal surgery, and bioimaging.     

The short-time self-diffusion coefficient of a sphere in a suspension of rigid rods

The short-time self-diffusion coefficient of a sphere in a suspension of rigid rods is calculated in first order in the rod volume fraction phi. For low rod concentrations, the correction to the Einstein diffusion constant of the sphere due to the presence of rods is a linear function of phi with the slope alpha proportional to the equilibrium averaged mobility diminution trace of the sphere interacting with a single freely translating and rotating rod. The two-body hydrodynamic interactions are calculated using the so-called bead model in which the rod of aspect ratio p is replaced by a stiff linear chain of touching spheres. The interactions between spheres are calculated using the multipole method with the accuracy controlled by a multipole truncation order and limited only by the computational power. A remarkable accuracy is obtained already for the lowest truncation order, which enables calculations for very long rods, up to p=1000. Additionally, the bead model is checked by filling the rod with smaller spheres. This procedure shows that for longer rods the basic model provides reasonable results varying less than 5% from the model with filling. An analytical expression for alpha as a function of p is derived in the limit of very long rods. The higher order corrections depending on the applied model are computed numerically. An approximate expression is provided, valid for a wide range of aspect ratios.