Use of Cartesian Rotation Vectors in Brownian Dynamics Algorithms: Theory and Simulation Results

Summary: We have shown that the components of Cartesian rotation vectors can be used successfully as generalized coordinates describing angular orientation in Brownian dynamics simulations of non‐spherical nanoparticles. For this particular choice of generalized coordinates, we rigorously derived the conformation‐space diffusion equations from kinetic theory for both free nanoparticles and nanoparticles interconnected by springs or holonomic constraints into polymer chains. The equivalent stochastic differential equations were used as a foundation for the Brownian dynamics algorithms. These new algorithms contain singularities only for points in the conformation‐space where both the probability density and its first coordinate derivative equal zero (weak singularities). In addition, the coordinate values after a single Brownian dynamics time step are throughout the conformation‐space equal to the old coordinate values plus the respective increments. For some parts of the conformation‐space these features represent a major improvement compared to the situation when Eulerian angles describe rotational dynamics. The presented simulation results of the equilibrium probability density for free nanoparticles are in perfect agreement with the results from kinetic theory.

Simulation of p^(eq)(Φ) for free nanoparticles.     

Singularity‐Free Brownian Dynamics Analyses of Rotational Dynamics: Non‐Spherical Nanoparticles in Solution

Summary: From kinetic theory we have rigorously derived singularity‐free Brownian dynamics analyses of nanoparticle rotational dynamics. The rigid non‐spherical nanoparticles incorporate all three rotational degrees of freedom. This was achieved by using the components of Cartesian rotation vectors as the generalized coordinates describing angular orientation. The new results constitute an important advance compared to the situation when Eulerian angles specify angular orientation. Our finding eliminates one of the main longstanding obstacles to detailed studies of nanoparticle rotational dynamics in the diffusion time domain. The described formalism is applicable to a wide range of nanoparticle systems including liquid crystals, biopolymers, and colloids.

     

Internal coordinate density of state from molecular dynamics simulation

The vibrational density of states (DoS), calculated from the Fourier transform of the velocity autocorrelation function, provides profound information regarding the structure and dynamic behavior of a system. However, it is often difficult to identify the exact vibrational mode associated with a specific frequency if the DoS is determined based on velocities in Cartesian coordinates. Here, the DoS is determined based on velocities in internal coordinates, calculated from Cartesian atomic velocities using a generalized Wilson's B ‐matrix. The DoS in internal coordinates allows for the correct detection of free dihedral rotations that may be mistaken as hindered rotation in Cartesian DoS. Furthermore, the pronounced enhancement of low frequency modes in Cartesian DoS for macromolecules should be attributed to the coupling of dihedral and angle motions. The internal DoS, thus deconvolutes the internal motions and provides fruitful insights to the dynamic behaviors of a system. © 2015 Wiley Periodicals, Inc.     

Adaptive biasing force method for scalar and vector free energy calculations

In free energy calculations based on thermodynamic integration, it is necessary to compute the derivatives of the free energy as a function of one (scalar case) or several (vector case) order parameters. We derive in a compact way a general formulation for evaluating these derivatives as the average of a mean force acting on the order parameters, which involves first derivatives with respect to both Cartesian coordinates and time. This is in contrast with the previously derived formulas, which require first and second derivatives of the order parameter with respect to Cartesian coordinates. As illustrated in a concrete example, the main advantage of this new formulation is the simplicity of its use, especially for complicated order parameters. It is also straightforward to implement in a molecular dynamics code, as can be seen from the pseudocode given at the end. We further discuss how the approach based on time derivatives can be combined with the adaptive biasing force method, an enhanced sampling technique that rapidly yields uniform sampling of the order parameters, and by doing so greatly improves the efficiency of free energy calculations. Using the backbone dihedral angles Phi and Psi in N-acetylalanyl-N'-methylamide as a numerical example, we present a technique to reconstruct the free energy from its derivatives, a calculation that presents some difficulties in the vector case because of the statistical errors affecting the derivatives.     

Use of the Rotation Vector in Brownian Dynamics Simulation of Transient Electro‐Optical Properties

We have recently developed a new singularity‐free algorithm for Brownian dynamics simulation of free rotational diffusion. The algorithm is rigorously derived from kinetic theory and makes use of the Cartesian components of the rotation vector as the generalized coordinates describing angular orientation. Here, we report on the application of this new algorithm in Brownian dynamics simulations of transient electro‐optical properties. This work serves two main purposes. Firstly, it demonstrates the integrity of the new algorithm for BD‐simulations of the most common transient electro‐optic experiments. Secondly, it provides new insight into the performance of the new algorithm compared to algorithms that make use of the Euler angles. We study the transient electrically induced birefringence in dilute solutions of rigid particles with anisotropic polarization tensor in response to external electric field pulses. The use of both one single electric pulse and two electric pulses with opposite polarity are being analyzed. We document that the new singularity‐free algorithm performs flawlessly. We find that, for these types of systems, the new singularity‐free algorithm, in general, outperforms similar algorithms based on the Euler angles. In a wider perspective, the most important aspect of this work is that it serves as an important reference for future development of efficient BD‐algorithms for studies of more complex systems. These systems include polymers consisting of rigid segments with single‐segment translational–rotational coupling, segment–segment fluid‐dynamic interactions and holonomic constraints.

     

Spherical particle in Poiseuille flow between planar walls

We study a spherical mesoparticle suspended in Newtonian fluid between plane-parallel walls with incident Poiseuille flow. Using a two-dimensional Fourier transform technique we obtain a symmetric analytic expression for the Green tensor for the Stokes equations describing the creeping flow in this geometry. From the matrix elements of the Green tensor with respect to a complete vector harmonic basis, we obtain the friction matrix for the sphere. The calculation of matrix elements of the Green tensor is done in large part analytically, reducing the evaluation of these elements to a one-dimensional numerical integration. The grand resistance and mobility matrices in Cartesian form are given in terms of 13 scalar friction and mobility functions which are expressed in terms of certain matrix elements calculated in the spherical basis. Numerical calculation of these functions is shown to converge well and to agree with earlier numerical calculations based on boundary collocation. For a channel width broad with respect to the particle radius, we show that an approximation defined by a superposition of single-wall functions is reasonably accurate, but that it has large errors for a narrow channel. In the two-wall geometry the friction and mobility functions describing translation-rotation coupling change sign as a function of position between the two walls. By Stokesian dynamics calculations for a polar particle subject to a torque arising from an external field, we show that the translation-rotation coupling induces sideways migration at right angles to the direction of fluid flow.     

Using redundant internal coordinates to optimize equilibrium geometries and transition states

A redundant internal coordinate system for optimizing molecular geometries is constructed from all bonds, all valence angles between bonded atoms, and all dihedral angles between bonded atoms. Redundancies are removed by using the generalized inverse of the G matrix; constraints can be added by using an appropriate projector. For minimizations, redundant internal coordinates provide substantial improvements in optimization efficiency over Cartesian and nonredundant internal coordinates, especially for flexible and polycyclic systems. Transition structure searches are also improved when redundant coordinates are used and when the initial steps are guided by the quadratic synchronous transit approach. © 1996 by John Wiley & Sons, Inc.     

Structural Properties of Nitromethane Molecular Crystal under High Pressure： An ab initio Investigation

The pressure-dependent structural properties under hydrostatic pressure up to 120 GPa and the decomposition under uniaxial compression along the b lattice vector up to 40 GPa of nitromethane molecular crystal using ab initio method are presented. The internal molecular bond lengths and bond angles were calculated for different pressures. All bond lengths decrease as the pressures are increased under hydrostatic compression. The obvious rotation of methyl group is 33.89° under hydrostatic pressure at 120 GPa. In addition, we observe the change of C-H bonds, which have been stretched under uniaxial compression along b lattice vector in the range of 0-40 GPa of nitromethane.     

Geometry optimization of periodic systems using internal coordinates

An algorithm is proposed for the structural optimization of periodic systems in internal (chemical) coordinates. Internal coordinates may include in addition to the usual bond lengths, bond angles, out-of-plane and dihedral angles, various "lattice internal coordinates" such as cell edge lengths, cell angles, cell volume, etc. The coordinate transformations between Cartesian (or fractional) and internal coordinates are performed by a generalized Wilson B-matrix, which in contrast to the previous formulation by Kudin et al. [J. Chem. Phys. 114, 2919 (2001)] includes the explicit dependence of the lattice parameters on the positions of all unit cell atoms. The performance of the method, including constrained optimizations, is demonstrated on several examples, such as layered and microporous materials (gibbsite and chabazite) as well as the urea molecular crystal. The calculations used energies and forces from the ab initio density functional theory plane wave method in the projector-augmented wave formalism.     

Small-amplitude protein conformational dynamics: Second-order analytic relation between cartesian coordinates and dihedral angles

An analytic expression for protein atomic displacements in Cartesian coordinate space (CCS) against small changes in dihedral angles is derived. To study time-dependent dynamics of a native protein molecule in CCS from dynamics in the internal coordinate space (ICS), it is necessary to convert small changes of internal coordinate variables to Cartesian coordinate variables. When we are interested in molecular motion, six degrees of freedom for translational and rotational motion of the molecule must be eliminated in this conversion, and this conversion is achieved by requiring the Eckart condition to hold. In this article, only dihedral angles are treated as independent internal variables (i.e., bond angles and bond lengths are fixed), and Cartesian coordinates of atoms are given analytically by a second-order Taylor expansion in terms of small deviations of variable dihedral angles. Coefficients of the first-order terms are collected in the K matrix obtained previously by Noguti and Go (1983) (see ref. 2). Coefficients of the second-order terms, which are for the first time derived here, are associated with the (newly termed) L matrix. The effect of including the resulting quadratic terms is compared against the precise numerical treatment using the Eckart condition. A normal mode analysis (NMA) in the dihedral angle space (DAS) of the protein bovine pancreatic trypsin inhibitor (BPTI) has been performed to calculate shift of mean atomic positions and mean square fluctuations around the mean positions. The analysis shows that the second-order terms involving the L matrix have significant contributions to atomic fluctuations at room temperature. This indicates that NMA in CCS involves significant errors when applied for such large molecules as proteins. These errors can be avoided by carrying out NMA in DAS and by considering terms up to second order in the conversion of atomic motion from DAS to CCS. © 1995 by John Wiley & Sons, Inc.     

Stationary points on the aminomethanol potential energy surface

Summary In this work we study surface fitting equations for a rigid rotor model of aminomethanol. The energies were obtained from the GAUSSIAN88 package using 3-21G bases and fitted on a least square equation, thus generating a Fourier series expansion of the energy as a function of two dihedral angles. The dihedral angles chosen are those that represent rotation around the C-O and N-C axes in the first case, and rotation around C-O and inversion around the amino group in the second case. Results indicate that the hydroxyl hydrogen is subject to almost free rotation around the C-O axis. Further fully relaxed 6-31G* calculations were performed in order to qualify the results obtained for the rigid rotor model.     

Optical-activity anisotropy of ordered phases in cellulose acetates

The optical activity of anisotropic solutions of cellulose diacetate in nitromethane and dimethyl sulfoxide and of diand triacetate films with vitrified ordered structure has been studied. The systems under study are characterized by high specific optical rotation [α], suggesting formation of a cholesteric mesophase. The value of [α] is found to depend on the angle of rotation of the samples relative to the direction of the polarization vector of an incident light beam in the plane perpendicular to this beam (the anisotropy of optical activity). This dependence (indicatrix) shows an irregular pattern and, when plotted in Cartesian coordinates, can be described with a distorted sinusoid. The data on the resolution of indicatrices into harmonic constituents and isolation of contributions due to isotropic components and anisotropic components, each of which is determined by the structural element with the corresponding asymmetry, are analyzed.     

Direct dynamics simulations using Hessian-based predictor-corrector integration algorithms

In previous research [J. Chem. Phys. 111, 3800 (1999)] a Hessian-based integration algorithm was derived for performing direct dynamics simulations. In the work presented here, improvements to this algorithm are described. The algorithm has a predictor step based on a local second-order Taylor expansion of the potential in Cartesian coordinates, within a trust radius, and a fifth-order correction to this predicted trajectory. The current algorithm determines the predicted trajectory in Cartesian coordinates, instead of the instantaneous normal mode coordinates used previously, to ensure angular momentum conservation. For the previous algorithm the corrected step was evaluated in rotated Cartesian coordinates. Since the local potential expanded in Cartesian coordinates is not invariant to rotation, the constants of motion are not necessarily conserved during the corrector step. An approximate correction to this shortcoming was made by projecting translation and rotation out of the rotated coordinates. For the current algorithm unrotated Cartesian coordinates are used for the corrected step to assure the constants of motion are conserved. An algorithm is proposed for updating the trust radius to enhance the accuracy and efficiency of the numerical integration. This modified Hessian-based integration algorithm, with its new components, has been implemented into the VENUS/NWChem software package and compared with the velocity-Verlet algorithm for the H(2)CO-->H(2)+CO, O(3)+C(3)H(6), and F(-)+CH(3)OOH chemical reactions.     

Quasi-Hamiltonian equations of motion for internal coordinate molecular dynamics of polymers

Conventional molecular dynamics simulations of macromolecules require long computational times because the most interesting motions are very slow compared to the fast oscillations of bond lengths and bond angles that limit the integration time step. Simulation of dynamics in the space of internal coordinates, that is, with bond lengths, bond angles, and torsions as independent variables, gives a theoretical possibility of eliminating all uninteresting fast degrees of freedom from the system. This article presents a new method for internal coordinate molecular dynamics simulations of macromolecules. Equations of motion are derived that are applicable to branched chain molecules with any number of internal degrees of freedom. Equations use the canonical variables and they are much simpler than existing analogs. In the numerical tests the internal coordinate dynamics are compared with the traditional Cartesian coordinate molecular dynamics in simulations of a 56 residue globular protein. For the first time it was possible to compare the two alternative methods on identical molecular models in conventional quality tests. It is shown that the traditional and internal coordinate dynamics require the same time step size for the same accuracy and that in the standard geometry approximation of amino acids, that is, with fixed bond lengths, bond angles, and rigid aromatic groups, the characteristic step size is 4 fs, which is 2 times higher than with fixed bond lengths only. The step size can be increased up to 11 fs when rotation of hydrogen atoms is suppressed. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18 : 1354–1364, 1997     

Torque and atomic forces for Cartesian tensor atomic multipoles with an application to crystal unit cell optimization

New equations for torque and atomic force are derived for use in flexible molecule force fields with atomic multipoles. The expressions are based on Cartesian tensors with arbitrary multipole rank. The standard method for rotating Cartesian tensor multipoles and calculating torque is to first represent the tensor with n indexes and 3ⁿ redundant components. In this work, new expressions for directly rotating the unique (n + 1)(n + 2)/2 Cartesian tensor multipole components Θ_pqr are given by introducing Cartesian tensor rotation matrix elements X( R ). A polynomial expression and a recursion relation for X( R ) are derived. For comparison, the analogous rotation matrix for spherical tensor multipoles are the Wigner functions D( R ). The expressions for X( R ) are used to derive simple equations for torque and atomic force. The torque and atomic force equations are applied to the geometry optimization of small molecule crystal unit cells. In addition, a discussion of computational efficiency as a function of increasing multipole rank is given for Cartesian tensors. © 2016 Wiley Periodicals, Inc.