Constrained molecular dynamics in the isothermal-isobaric ensemble and its adaptation for adiabatic free energy dynamics

The implementation of holonomic constraints within measure-preserving integrators for molecular dynamics simulations in the isothermal-isobaric ensemble is considered. We review the basic methodology of generating measure-preserving integrators for the microcanonical, canonical, and isothermal-isobaric ensembles and proceed to show how the standard SHAKE and RATTLE algorithms must be modified for the isothermal-isobaric ensemble. Comparison is made between constrained and unconstrained simulations employing multiple time scale integration techniques. Finally, we describe a temperature accelerated version of the isothermal-isobaric molecular dynamics approach, in which the cell matrix is adiabatically decoupled from the particles and maintained at a high temperature as a means of exploring polymorphism in molecular crystals. We demonstrate that constraints can be easily adapted for this new approach and, again, we compare the performace of this temperature-accelerated scheme with and without bond constraints.     

Molecular dynamics of polymers with explicit but frozen hydrogens

A new form of holonomic constraint, called a conic constraint, is introduced for the purpose of eliminating the fast vibrations of hydrogen atoms in molecular simulations of systems of aliphatic chains. It can easily be combined with bond constraints in SHAKE/RATTLE algorithms for which a unified tolerance criterion is defined. The new form of constraint allows the use of rather large time steps (in the 2–3 fs range). The procedure is illustrated for a full atomic model of polypropylene.     

Brownian dynamics with constraints

Systems possessing degrees of freedom operating on widely separated timescales, where the effects of those operating on the smaller timescales are relatively unimportant, may be modelled by the use of the Langevin equation. In order to study such systems containing complex polyatomic particles, holonomic constraints may be used. Though there is no lack of published algorithms for the numerical solution of the Langevin equation, few of them have been developed with sufficient rigour to ensure their precision, nor to demonstrate their compatibility with constraints. This study recapitulates an approach based upon Runge-Kutta equations which has the advantage of being perfectible to any desired order in the time-step, and shows how it may be combined with the SHAKE method in order to perform constrained Brownian dynamics simulations. Results are presented for some simple systems with a third order algorithm, and it is found that the correct dynamic and statistical behaviour is recovered.     

Parallelization in time of numerical simulations of fully-developed plasma turbulence using the parareal algorithm

It is shown that numerical simulations of fully-developed plasma turbulence can be successfully parallelized in time using the parareal algorithm. The result is far from trivial, and even unexpected, since the exponential divergence of Lagrangian trajectories as well as the extreme sensitivity to initial conditions characteristic of turbulence set these type of simulations apart from the much simpler systems to which the parareal algorithm has been applied to this day. It is also shown that the parallel gain obtainable with this method is very promising (close to an order of magnitude for the cases and implementations described), even when it scales with the number of processors quite differently to what is typical for spatial parallelization.     

Editorial – Topical Review

Starting from a family of equations of motion for the dynamics of extended systems whose trajectories sample constant pressure and temperature ensemble distributions (Ferrario, M., 1993, in Computer Simulation in Chemical Physics, edited by M. P. Allen and D. J. Tildesley (Dordrecht: Kluwer)), explicit time reversible integration schemes are derived through a straightforward Trotter factorization of the dynamic Liouville propagator, along the lines first described by Tuckerman, M., Martyna, G. J., and Berne, B. J., 1992, J. chem. Phys., 97, 1990. The original Andersen's constant-pressure dynamics are recovered in the limit of zero coupling with the Nose thermostat. Reversible integration schemes are derived as a generalization of the velocity Verlet algorithm, with direct handling of the velocity dependent forces in such a way that both predictions and relative iterative corrections are not required. For the sake of clarity both the equations of motion and the Trotter factorization are kept to the basic level. The proposed structure can accommodate easily, when needed, complications such as multiple timesteps and more effective thermostats (Nose-Hoover-chain). Finally, an application is made to a model molecular system subjected to holonomic constraints by means of the SHAKE algorithm. In the constant pressure case it is no longer possible to avoid using a prediction for the constraint contribution to the volume acceleration; however, recourse to a minimal iteration scheme still achieves excellent overall behaviour for the proposed integration algorithm, with no perceptible difference from the unconstrained case.     

面向结构网格自适应并行计算的矩形区域求差集快速算法

结构网格自适应程序需要使用矩形区域求差集算法计算网格层间数据依赖关系和网格层嵌套关系.原有的矩形区域求差集算法时间复杂度较高,成为该类应用大规模并行计算可扩展性能瓶颈.本文利用分而治之的方法,构造近似线性时间复杂度的矩形区域求差集快速算法,并利用区域分解实现该算法的并行计算.分别针对规则矩形区域和多层自适应网格的非规则矩形区域求差集问题,验证该算法的效率.结果表明,该算法具有近似线性计算复杂度,对于大规模计算问题,加速效果显著.     

Small-angle Coulomb collision model for particle-in-cell simulations

We construct and investigate a set of stochastic differential equations that incorporate the physics of velocity-dependent small-angle Coulomb collisions among the plasma particles in a particle-in-cell simulation. Each particle is scattered stochastically from all the other particles in a simulation cell modeled as one or more Maxwellians. Total energy and momentum are conserved by linear transformation of the velocity increments. In two test simulations the proposed “particle-moment” collision algorithm performs well with time steps as large as 10% of the relaxation time – far larger than a particle-pairing collision algorithm, in which pairs of particles are scattered from one another, requires to achieve the same accuracy.     

Calculation of thermal effects in the photodesorption of NO from NiO(100)

The influence of temperature in the photodesorption of NO from a NiO(100)-surface is studied with a two-dimensional quantum wave packet approach. The complete process including laser-induced excitation and subsequent relaxation is treated explicitly from first principles. The electronic quenching caused by the interaction of the excited adsorbate–substrate system with electron–hole-pairs in the surface is treated with the surrogate Hamiltonian approach. We have implemented a parallelization scheme of the wave packet propagation based on a one-dimensional data decomposition to perform simulations in a reasonable computing time. The results are compared with mixed quantum-classical simulations and with experimental measurements. Both desorption yield and mean velocity of the desorbing molecules were enhanced with increasing temperature. The calculated rotational temperatures are consistent with experimental results. PACS  79.20.La; 68.43.Tj     

亚乙基硫脲分子激发拉曼虚态的研究

本文从拉曼峰强入手,求得了亚乙基硫脲(ETU)分子的"时间分辨键极化率",并讨论了该分子的激发拉曼虚态性质,发现了该分子"激发虚态电子向分子外围键流动"、"电子弛豫后的键极化率分布与基态电子的密度分布相似"、"不同激发波长下的键极化率衰减时间满足不确定关系"等特点。本文还研究了该分子的表面增强拉曼光谱,指出"电荷转移机制"的极化率弛豫时间长于"电磁增强机制"等与表面增强拉曼效应相关的结论。     

Hydrogen-bond dynamics of interior and surface molecules in a water nanocluster: temperature and size effects

Molecular dynamics simulations are employed to investigate the effects of temperature and size on the hydrogen-bond dynamics of interior molecules and surface molecules in a water nanocluster. The flexible three-centred (F3C) water model is invoked in the simulations. To inspect the dynamics of the interior hydrogen bonds and the surface hydrogen bonds, a spherical water nanocluster is modelled and then divided into interior molecules and surface molecules according to the density profile of the water nanocluster. It is observed that at higher temperatures the average number of hydrogen bonds decreases and yields faster hydrogen-bond relaxation for both interior molecules and surface molecules of the water nanocluster. Furthermore, the surface molecules have a lower average number of hydrogen bonds than the interior molecules. The lifetime of the surface hydrogen bonds is slightly longer than that of the interior hydrogen bonds, whereas the hydrogen-bond structural relaxation time of the surface molecules is more obviously slower than that of the interior molecules. Regarding the size effect, a larger water nanocluster is seen to have a larger average number of hydrogen bonds and a longer hydrogen-bond structural relaxation time.     

The relaxations of temporal bond polarizabilities of methylviologen adsorbed on the Ag electrode by 514.5 nm excitation: a Raman intensity study

An algorithm to elucidate the temporal bond polarizabilities from the surface enhanced Raman (SERS) intensities was employed to the case of methylviologen (MV) adsorbed on the Ag electrode. This enables us to obtain the properties of its SERS mechanisms and the effect of its adsorption. The analysis shows that the charge transfer and electromagnetic mechanisms involving in this MV SERS system possess different relaxation times for its various temporal bond polarizabilities. The physics is that the process involved in the charge transfer mechanism will take longer time than that involved in the electromagnetic mechanism since it needs more time to redistribute the charges during relaxation. The time division between these two mechanisms is figured out to be around 3 ps for this system. Adsorption also enhances the relaxation of the temporal bond polarizabilities, in general. The adsorption effect is indicated by the temporal bond polarizabilities close to the final stage of relaxation. They are, in fact, the quantities parallel to the bond electronic densities in the molecular orbital (MO) concept. For comparison, the case of MV solid was also analyzed. Copyright © 2008 John Wiley & Sons, Ltd.     

Algorithm for molecular dynamics simulations of spin liquids

A new symplectic time-reversible algorithm for numerical integration of the equations of motion in magnetic liquids is proposed. It is tested and applied to molecular dynamics simulations of a Heisenberg spin fluid. We show that the algorithm exactly conserves spin lengths and can be used with much larger time steps than those inherent in standard predictor-corrector schemes. The results obtained for time correlation functions demonstrate the evident dynamic interplay between the liquid and magnetic subsystems.     

Advanced flicker spectroscopy of fluid membranes

The bending elasticity of a fluid membrane is characterized by its modulus and spontaneous curvature. We present a new method, advanced flicker spectroscopy of giant nonspherical vesicles, which makes it possible to simultaneously measure both parameters for the first time. Our analysis is based on the generation of a large set of reference data from Monte Carlo simulations of randomly triangulated surfaces. As an example of the potential of the procedure, we monitor thermal trajectories of vesicle shapes and discuss the elastic response of zwitterionic membranes to transmembrane pH gradients. Our technique makes it possible to easily characterize membrane curvature as a function of environmental conditions.     

Compression of carbon nanotubes filled with C60, CH4, or Ne: predictions from molecular dynamics simulations

The effect of filling nanotubes with C60, CH4, or Ne on the mechanical properties of the nanotubes is examined. The approach is classical molecular dynamics using the reactive empirical bond order (REBO) and the adaptive intermolecular REBO potentials. The simulations predict that the buckling force of filled nanotubes can be larger than that of empty nanotubes, and the magnitude of the increase depends on the density of the filling material. In addition, these simulations demonstrate that the buckling force of empty nanotubes depends on temperature. Filling the nanotube disrupts this temperature effect so that it is no longer present in some cases.     

Role of protein-water hydrogen bond dynamics in the protein dynamical transition

The role of water in protein dynamics has been investigated using molecular dynamics simulations of crystals and a dehydrated powder. On the 100 ps time scale, the anharmonic and diffusive motions involved in the protein structural relaxation are correlated with the protein-water hydrogen bond dynamics. The complete structural relaxation of the protein requires relaxation of the hydrogen bond network via solvent translational displacement. Inhibiting the solvent translational mobility, and therefore the protein-water hydrogen bond dynamics, has an effect on the protein relaxation similar to dehydration.     

An iterative method for improving the quality of reconstruction of a three-dimensional surface

A complex image with constraints imposed on the amplitude and phase image components is processed using the Gerchberg iterative algorithm for the first time. The use of the Gerchberg iterative algorithm makes it possible to improve the quality of a three-dimensional surface profile reconstructed by the previously proposed method that is based on the multiangle projection of fringes and the joint processing of the obtained images by Fourier synthesis.     

亚乙基硫脲分子键伸缩模式与全耦合模式键极化率的对比及其拉曼激发虚态的相关研究

从拉曼峰强着手,得到了键伸缩模式与全耦合模式两种不同计算条件下,亚乙基硫脲(ethylene thiourea,ETU)分子的键极化率,并比较分析了两种计算结果的异同.研究表明:在键伸缩模式的算法中,仅考虑势能分布中键对称伸缩比重相对较大的部分拉曼峰参与极化率计算,显然忽略了键伸缩与键弯曲间的相互耦合,造成了与弯曲振动耦合较强的部分键伸缩极化率值存在一定误差.因此,此方法虽能使问题简化,但却丢失了一些信息.而全耦合模式算法考虑了所有振动模式(键伸缩与键弯曲)相互影响的情况,能更全面的反映键电荷的分布情况,但却会使问题的求解过程变得复杂.同时,在拉曼激发虚态的相关研究中,两种计算方法却得到几乎一致的键极化率弛豫方式(呈单指数规律衰减)及相同的衰减特征时间.     

A classical hamiltonian method for calculating IR-Resonant photodesorption rates

Summary In this work, classical nonlinear dynamical methods are used to study photodesorption induced by IR-laser resonant excitation of an internal vibrational mode of a physisorbed molecule. Starting from the large disagreement between the experimental values of the photodesorption rates and the theoretical quantum perturbative ones for CH₃F/NaCl, the Elastic-Resonant-State Decay (ERSD) model is analysed focusing attention on the nonexponential time decay behaviour of the survival probability of the system, that shows two different time scales. This is obtained by a suitable average of classical trajectories calculated by integrating Hamilton’s equations of motion. The Poincaré sections are used to clarify the details of the microscopic dynamics of the IR-photodesorption, pointing out the dynamical role of cantori as responsible for desorption occurring on a longer time scale. A method to build up the separatrix for the ERSD system is proposed, stressing that such a separatrix acts as a bottleneck against the energy flux exchanged between the physisorption and the intramolecular bond of the adsorbate. A calculation of the photodesorption rate when a condition of fast relaxation between the two bonds is verified, is performed using the Transport Theory in Hamiltonian System. The values of the photodesorption rate obtained by this approach are in good agreement with those obtained by integrating numerically the equations of motion.     

并行MapReduce PLS算法及其在光谱分析中的应用

偏最小二乘(PLS)算法是常用的光谱建模算法,然而对于海量光谱处理情形,在单台计算机上建模及优化时间开销很大。基于MapReduce编程模式,提出了并行MapReduce PLS回归算法,包括并行数据标准化和并行主成分提取两个过程。在多台普通计算机上搭建Hadoop云计算集群平台,以近红外光谱处理为例,开展了算法验证实验。实验结果表明,基于MapReduce编程模式的并行PLS算法对海量近红外光谱数据集进行回归建模时,能有效提高建模速度,随计算机台数的增多可得到接近线性的加速比,并具有良好的扩展性。     

Reconstruction of the fine structure of an acoustic scatterer against the distorting influence of its large-scale inhomogeneities

In the ultrasonic diagnostics of small-size neoplasms of biological tissues at the earliest stage of their development, an efficient way to eliminate the distorting influence of high-contrast or large inhomogeneities of the biological medium is to apply the iterative technique. A simple approach is proposed, which makes it possible with only two iteration steps to achieve an efficient focusing of the tomograph array. At the first step, the unknown distribution of the large-scale inhomogeneities of sound velocity and absorption over the scatterer is reconstructed, where the large-scale inhomogeneities are those whose size exceeds several wavelengths. At the second step, the fine structure of the scatterer is reconstructed against the large-scale background, which can be performed with a high accuracy owing to the evaluation of the background at the first step. The possibility of simultaneous reconstruction of the large-scale and fine structures by the noniterative Grinevich-Novikov algorithm is considered as an alternative. This algorithm reconstructs in an explicit form two-dimensional refractive-absorbing acoustic scatterers of almost arbitrary shape and strength. Taking into account the effects of multiple scattering, this algorithm provides resolution of the fine structure almost as good as that achieved in reconstructing the same structure against an undistorting homogeneous background. The results of numerical simulations of both algorithms are presented.