Comparative structural studies of psychrophilic and mesophilic protein homologues by molecular dynamics simulation

Comparative molecular dynamics simulations of psychrophilic type III antifreeze protein from the North-Atlantic ocean-pout Macrozoarces americanus and its corresponding mesophilic counterpart, the antifreeze-like domain of human sialic acid synthase, have been performed for 10 ns each at five different temperatures. Analyses of trajectories in terms of secondary structure content, solvent accessibility, intramolecular hydrogen bonds and protein–solvent interactions indicate distinct differences in these two proteins. The two proteins also follow dissimilar unfolding pathways. The overall flexibility calculated by the trace of the diagonalized covariance matrix displays similar flexibility of both the proteins near their growth temperatures. However at higher temperatures psychrophilic protein shows increased overall flexibility than its mesophilic counterpart. Principal component analysis also indicates that the essential subspaces explored by the simulations of two proteins at different temperatures are non-overlapping and they show significantly different directions of motion. However, there are significant overlaps within the trajectories and similar directions of motion of each protein especially at 298 K, 310 K and 373 K. Overall, the psychrophilic protein leads to increased conformational sampling of the phase space than its mesophilic counterpart.Our study may help in elucidating the molecular basis of thermostability of homologous proteins from two organisms living at different temperature conditions. Such an understanding is required for designing efficient proteins with characteristics for a particular application at desired working temperatures.     

Temperature-induced unfolding pathway of a type III antifreeze protein: insight from molecular dynamics simulation

Molecular dynamics simulations of the temperature-induced unfolding reaction of a cold-adapted type III antifreeze protein (AFPIII) from the Antarctic eelpout Lycodichthys dearborni have been carried out for 10 ns each at five different temperatures. While the overall character and order of events in the unfolding process are well conserved across temperatures, there are substantial differences in the timescales over which these events take place. Plots of backbone root mean square deviation (RMSD) against radius of gyration (Rg) serve as phase space trajectories. These plots also indicate that the protein unfolds without many detectable intermediates suggestive of two-state unfolding kinetics. The transition state structures are identified from essential dynamics, which utilizes a principal component analysis (PCA) on the atomic fluctuations throughout the simulation. Overall, the transition state resembles an expanded native state with the loss of the three 3₁₀ helices and disrupted C-terminal region.Our study provides insight into the structure–stability relationship of AFPIII, which may help to engineer AFPs with increased thermal stability that is more desirable than natural AFPs for some industrial and biomedical purposes.     

分子动力学模拟的柔性对接

分子自动对接技术在过去二十年里取得很大发展和成功,但是仍然面对如何处理分子柔性这样一个难题。这篇综述概要介绍分子柔性对接技术的进展并重点介绍分子动力学模拟技术。     

A parallel algorithm for molecular dynamics simulation of branched molecules

To get an insight into the effects of molecular architecture in the behaviour of thin lubricant films we have devised an algorithm for simulation of branched molecules. We have used this algorithm successfully to simulate branched isomers of C₃₀. However the algorithm is flexible enough to be used for the simulation of more complex branched molecules. The resulting algorithm can be used in molecular dynamics simulation of branched molecules and could be helpful in designing new materials at the molecular level.     

Molecular dynamics simulation of ribosome jam

We propose a coarse-grained molecular dynamics model of ribosome molecules to study the dependence of translation process on environmental parameters. We found the model exhibits traffic jam property, which is consistent with an ASEP model. We estimated the influence of the temperature and concentration of molecules on the hopping probability used in the ASEP model. Our model can also treat environmental effects on the translation process that cannot be explained by such cellular automaton models.     

分子动力学模拟HEWL晶体在不同环境中的动力学行为

简要阐述分子动力学模拟的原理及步骤,介绍研究溶菌酶的一般方法和优缺点。在Ubuntu操作系统环境下,利用Gromacs软件和其自带的Gromos96力场,通过分子动力学模拟(MD)鸡蛋清溶菌酶晶体(chicken egg-whitelysozyme,HEWL)溶液,考察真空、水溶液和加入NaCl 3种不同环境条件对溶菌酶晶体构象动力学行为的影响,发现无论从均方根位移(rmsd)、回旋半径、还是从B因子值的轨迹图分析,HEWL在水溶液特别是加入抗衡离子(Na~+,Cl~-)的水溶液的环境下的结构更稳定、合理,与(protein data bank)数据库的真实情况相符。原因是Cl~-与溶菌酶晶体在界面处发生了吸附现象,局部形成溶菌酶-Cl~-复合物,抑制了蛋白-水合物中水分子在相邻水合位置间的跳跃,从而使单晶体在离子液态中更加稳定。模拟结果表明,在pH值6.5,等电位点13.1,总电荷7.999 6的体系下,影响HEWL的吸附位点为123号残基(色氨酸),对从分子水平上解释HEWL晶体的动力学吸附行为具有重要指导意义。     

Huge-scale molecular dynamics simulation of multibubble nuclei

We have developed molecular dynamics codes for a short-range interaction potential that adopt both the flat-MPI and MPI/OpenMP hybrid parallelizations on the basis of a full domain decomposition strategy. Benchmark simulations involving up to 38.4 billion Lennard-Jones particles were performed on Fujitsu PRIMEHPC FX10, consisting of 4800 SPARC64 IXfx 1.848 GHz processors, at the Information Technology Center of the University of Tokyo, and a performance of 193 teraflops was achieved, which corresponds to a 17.0% execution efficiency. Cavitation processes were also simulated on PRIMEHPC FX10 and SGI Altix ICE 8400EX at the Institute of Solid State Physics of the University of Tokyo, which involved 1.45 billion and 22.9 million particles, respectively. Ostwald-like ripening was observed after the multibubble nuclei. Our results demonstrate that direct simulations of multiscale phenomena involving phase transitions from the atomic scale are possible and that the molecular dynamics method is a promising method that can be applied to petascale computers.     

Scalability study of molecular dynamics simulation on Godson-T many-core architecture

Molecular dynamics (MD) simulation has broad applications, and an increasing amount of computing power is needed to satisfy the large scale of the real world simulation. The advent of the many-core paradigm brings unprecedented computing power, but it remains a great challenge to harvest the computing power due to MD’s irregular memory-access pattern. To address this challenge, this paper presents a joint application/architecture study to enhance the scalability of MD on Godson-T-like many-core architecture. First, a preprocessing approach leveraging an adaptive divide-and-conquer framework is designed to exploit locality through memory hierarchy with software controlled memory. Then three incremental optimization strategies–a novel data-layout to improve data locality, an on-chip locality-aware parallel algorithm to enhance data reuse, and a pipelining algorithm to hide latency to shared memory–are proposed to enhance on-chip parallelism for Godson-T many-core processor. Experiments on Godson-T simulator exhibit strong-scaling parallel efficiency of 0.99 on 64 cores, which is confirmed by a field-programmable gate array emulator. Also the performance per watt of MD on Godson-T is much higher than MD on a 16-cores Intel core i7 symmetric multiprocessor (SMP) and 26 times higher than MD on an 8-core 64-thread Sun T2 processor. Detailed analysis shows that optimizations utilizing architectural features to maximize data locality and to enhance data reuse benefit scalability most. Furthermore, a hierarchical parallelization scheme is designed to map the MD algorithm to Godson-T many-core cluster and a simple performance model is derived, which suggests that the optimization scheme is likely to scale well toward exascale. Certain architectural features are found essential for these optimizations, which could guide future hardware developments.     

Contributions of 3'-overhang to the dissociation of small interfering RNAs from the PAZ domain: molecular dynamics simulation study

RNA interference (RNAi) is a ‘knock-down’ reaction to reduce expression of a specific gene through highly regulated, enzyme-mediated processes. Small interfering RNAs (siRNAs) are RNA molecules that play an effector role in RNAi and can bind the PAZ domains present in Dicer and RISC. We investigated the interaction between the PAZ domain and the siRNA-like duplexes through dissociation molecular dynamics (DMD) simulations. Specifically, we focused on the response of the PAZ domain to various 3′-overhang structures of the siRNA-like duplexes. We found that the siRNA-like duplex with the 3′ UU-overhang made relatively more stable complex with the PAZ domain compared to those with 3′ CC-, AA-, and GG-overhangs. The siRNA-like duplex with UU-overhang was easily dissociated from the PAZ domain once the structural stability of the complex is impaired. Interestingly, the 3′ UU-overhang spent the least time at the periphery region of the binding pocket during the dissociation process, which can be mainly attributable to UU-overhang's smallest number of hydrogen bonds.     

Molecular docking/dynamics studies of Aurora A kinase inhibitors

The binding modes of a known 1,4,5,6-tetrahydropyrrolo[3,4-c]pyrazole, quinazoline, pyrimidine and indolinone series of Aurora A kinase inhibitors have been studied using molecular docking and molecular dynamics (MD) simulations. Crystallographic bound compound 8 was precisely predicted by our docking procedure as evident from 0.43 Å root mean square (rms) deviations. In addition compound 25 (AZ_68) has been successfully cross-docked within the Aurora A kinase active site, which was pre-organized for inhibitor 8. We found four key sites (A: solvent-exposed front pocket, B: hinge region, C: selectivity pocket and D: solvent-exposed phosphate binding region) of the Aurora A kinase contributing towards the binding of these compounds. We suggest that the small hydrophobic substituents at C-6 position of pyrrolopyrazole nucleus (in compounds 1–8); C-6 and C-7 positions of the quinazoline moiety (in compounds 9–23); C-2 position of the quinazoline and C-4 position of the pyrimidine (in compound 25) could be more effective and selective through increased hydrophobic contacts and selectivity pocket interactions with these modifications of Aurora A kinase inhibitors. Five representative complexes were subjected to 1000 ps of MD simulation to determine the stability of the predicted binding conformations. The low value of the root mean square deviations (ranging from 0.725 to 1.820 Å) between the starting complex structure and the energy minimized final average complex structure suggests that the Glide Extra Precision (XP) derived docked complexes are in a state of near equilibrium. The structure-based drug design strategy described in this study will be highly useful for the development of new inhibitors with high potency and selectivity.     

Effects of different solvents on the conformations of apoptotic cytochrome c: Structural insights from molecular dynamics simulation

Cytochrome c (cyt-c) upon binding with cardiolipin acquires peroxidase activity and is strictly connected to the pathogenesis of many human diseases including neurodegenerative and cardiovascular diseases. Interaction of cyt-c with cardiolipin mimics partial unfolding/conformational changes of cyt-c in different solvent environments. Dynamic pictures of these conformational changes of cyt-c are crucial in understanding their physiological roles in mitochondrial functions. Therefore, atomistic molecular dynamics (MD) simulations have been carried out to investigate the effect of different solvents (water, urea/water, MeOH and DMSO) on the structure and conformations of apoptotic cyt-c (Fe³⁺). Our study demonstrates that the structural changes in the protein are solvent dependent. The structural differences are observed majorly on the β-sheets and α-helical conformations and the degree of their perturbation are specific to the solvent. Although a complete loss of β-sheets (0%) is observed in MeOH and DMSO, by contrast, well preserved β-sheets (3.84%) are observed in water and urea/water. A significant decrease in the α-helical contents is observed in MeOH (41.34%) and water (42.46%), a negligible alteration in DMSO (44.25%) and well preserved α-helical (45.19%) contents in urea/water. The distances between the residues critical for electron transfer are decreased considerably for DMSO. Further, the reduction in residue flexibility and the conformational space indicate that the collective motions of cyt-c are reduced when compared to other cosolvents. Essential dynamics analysis implies that the overall motions of cyt-c in water, MeOH and urea/water are involved in three to four eigenvectors and in first eigenvector in DMSO. Overall, we believe that MD simulations of cyt-c in different solvents can provide a detailed microscopic understanding of the physiological roles, electron transport and peroxidase function in the early events of apoptosis which are hard to probe experiments.     

Coarse-grained molecular dynamics simulation of water diffusion in the presence of carbon nanotubes

Computational modeling of the translational diffusion of water molecules in anisotropic environments entails vital relevance to understand correctly the information contained in the magnetic resonance images weighted in diffusion (DWI) and of the diffusion tensor images (DTI). In the present work we investigated the validity, strengths and weaknesses of a coarse-grained (CG) model based on the MARTINI force field to simulate water diffusion in a medium containing carbon nanotubes (CNTs) as models of anisotropic water diffusion behavior. We show that water diffusion outside the nanotubes follows Ficḱs law, while water diffusion inside the nanotubes is not described by a Ficḱs behavior. We report on the influence on water diffusion of various parameters such as length and concentration of CNTs, comparing the CG results with those obtained from the more accurate classic force field calculation, like the all-atom approach. Calculated water diffusion coefficients decreased in the presence of nanotubes in a concentration dependent manner. We also observed smaller water diffusion coefficients for longer CNTs. Using the CG methodology we were able to demonstrate anisotropic diffusion of water inside the nanotube scaffold, but we could not prove anisotropy in the surrounding medium, suggesting that grouping several water molecules in a single diffusing unit may affect the diffusional anisotropy calculated. The methodologies investigated in this work represent a first step towards the study of more complex models, including anisotropic cohorts of CNTs or even neuronal axons, with reasonable savings in computation time.     

H-POSS表面润湿性的分子动力学模拟

采用分子动力学模拟方法研究了水滴在多面体低聚倍半硅氧烷(H-POSS)固体表面的润湿性能,H-POSS分子和水分子分别采用COMPASS力场和SPC力场模型。模拟得到H-POSS基体密度为1.84g/cm~3,且X射线衍射模拟发现基体具有明显衍射峰,表现出晶体特性,说明COMPASS力场适用于H-POSS基体的构建与研究。H-POSS表面水接触角的模拟值为104.9°,具有疏水性能。通过直接水解法由三氯硅烷(HSiCl_3)实验合成出H-POSS样品,傅立叶红外表征(FT-IR)发现,在2260、1142和871cm~(-1)波数位置出现吸收峰,证实了所合成的样品为H-POSS。其表面水接触角的实验值为109.3°,与模拟值的相对误差仅为4%,说明分子动力学方法可应用于计算H-POSS材料表面润湿性。模拟结果还表明体系温度影响H-POSS材料的表面润湿性,增大体系温度,表面疏水性能降低。     

Feature-preserving adaptive mesh generation for molecular shape modeling and simulation

We describe a chain of algorithms for molecular surface and volumetric mesh generation. We take as inputs the centers and radii of all atoms of a molecule and the toolchain outputs both triangular and tetrahedral meshes that can be used for molecular shape modeling and simulation. Experiments on a number of molecules are demonstrated, showing that our methods possess several desirable properties: feature-preservation, local adaptivity, high quality, and smoothness (for surface meshes). We also demonstrate an example of molecular simulation using the finite element method and the meshes generated by our method. The approaches presented and their implementations are also applicable to other types of inputs such as 3D scalar volumes and triangular surface meshes with low quality, and hence can be used for generation/improvement of meshes in a broad range of applications.     

The molecular dynamics of assembly of the ubiquitous aortic medial amyloidal medin fragment

In recent years there is an increased understanding of the molecular conformation of amyloid fibrils. However, much less is known about the early events that lead to the formation of these medically important assemblies. The clarification of these very important mechanistic details on the process may indicate directions towards the inhibition of the early stages of the assembly, where harmful species are most likely to form. Here, we study the dynamics of assembly of short amyloidogenic peptide fragments from the medin polypeptide. This polypeptide is of unique interest since amyloid deposits composed of medin are found almost in all the population above the age of 50. Twelve independent 50 ns long molecular dynamics simulations in explicit water have been run on peptide NH₂–NFGSVQFV–COOH, the minimal recognition hexapeptide element, NH₂–NFGSVQ–COOH, and several single-point mutants. In all cases a three-stranded polymeric β-sheet was used as the basic unit from which fibrils can be formed. Our results clearly indicate the need of well-defined sequence and stereochemical constraints to allow the formation of stable well-ordered aggregates. One of the key findings is the need for the presence of a phenylalanine residue, but not other hydrophobic amino acids, in specific positions within the peptide. Taken together, the results are consistent with recent high-resolution structures of amyloid assemblies and provide unique insights into the dynamics of these structures.     

非平衡分子动力学模拟在炭膜气体分离中的研究进展

叙述了近年来非平衡分子动力学在炭膜气体分离中的研究进展,包括分子模拟方法的选取及非平衡分子动力学模拟的几个主要步骤;综述了气体混合物通过炭膜时传递和分离性能的影响因素(如孔结构及孔径大小、操作压力、温度、组分组成等);并对炭膜研究现状作了简要的概述,提出目前分子模拟在炭膜研究中所存在的问题,一方面是模型建立问题．另一方面是研究膜传递和分离机理性能所考察的影响因素还不够完善,加强对膜传递分离机理的研究,对膜的设计、制备及应用都有重要的指导意义。     

Effect of allosteric molecules on structure and drug affinity of HIV-1 protease by molecular dynamics simulations

Recent experiments show that small molecules can bind onto the allosteric sites of HIV-1 protease (PR), which provides a starting point for developing allosteric inhibitors. However, the knowledge of the effect of such binding on the structural dynamics and binding free energy of the active site inhibitor and PR is still lacking. Here, we report 200 ns long molecular dynamics simulation results to gain insight into the influences of two allosteric molecules (1H-indole-6-carboxylic acid, 1F1 and 2-methylcyclohexano, 4D9). The simulations demonstrate that both allosteric molecules change the PR conformation and stabilize the structures of PR and the inhibitor; the residues of the flaps are sensitive to the allosteric molecules and the flexibility of the residues is pronouncedly suppressed; the additions of the small molecules to the allosteric sites strengthen the binding affinities of 3TL-PR by about 12–15 kal/mol in the binding free energy, which mainly arises from electrostatic term. Interestingly, it is found that the action mechanisms of 1F1 and 4D9 are different, the former behaviors like a doorman that keeps the inhibitor from escape and makes the flaps (door) partially open; the latter is like a wedge that expands the allosteric space and meanwhile closes the flaps. Our data provide a theoretical support for designing the allosteric inhibitor.     

An improved version of the Green's function molecular dynamics method

This work presents an improved version of the Green's function molecular dynamics method (Kong et al., 2009; Campañá and Müser, 2004 [1,2]), which enables one to study the elastic response of a three-dimensional solid to an external stress field by taking into consideration only atoms near the surface. In the previous implementation, the effective elastic coefficients measured at the Γ-point were altered to reduce finite size effects: their eigenvalues corresponding to the acoustic modes were set to zero. This scheme was found to work well for simple Bravais lattices as long as only atoms within the last layer were treated as Green's function atoms. However, it failed to function as expected in all other cases. It turns out that a violation of the acoustic sum rule for the effective elastic coefficients at Γ (Kong, 2010 [3]) was responsible for this behavior. In the new version, the acoustic sum rule is enforced by adopting an iterative procedure, which is found to be physically more meaningful than the previous one. In addition, the new algorithm allows one to treat lattices with bases and the Green's function slab is no longer confined to one layer.

New version program summary

Program title: FixGFC/FixGFMD v1.12Catalogue identifier: AECW_v1_1Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AECW_v1_1.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 206 436No. of bytes in distributed program, including test data, etc.: 4 314 850Distribution format: tar.gzProgramming language: C++Computer: AllOperating system: LinuxHas the code been vectorized or parallelized?: Yes. Code has been parallelized using MPI directives.RAM: Depends on the problemClassification: 7.7External routines: LAMMPS (http://lammps.sandia.gov/), MPI (http://www.mcs.anl.gov/research/projects/mpi/), FFT (http://www.fftw.org/)Catalogue identifier of previous version: AECW_v1_0Journal reference of previous version: Comput. Phys. Comm. 180 (2009) 1004Does the new version supersede the previous version?: YesNature of problem: Green's function molecular dynamics (GFMD) is a coarse-graining method that enables one to investigate the full elastic response of an interface between a semi-infinite solid and a contact while taking only the surface atoms in the solid into consideration. The effect of long-range elastic deformations on the surface atoms from the semi-infinite solid is replaced by effective elastic interactions, thus reducing the problem from three dimensions to two dimensions without compromising the physical essence of the problem.Solution method: See “Nature of problem”.Reasons for new version: The basic theory underlying the new version is essentially the same as the previous one, while the special treatment to reduce the finite size effect on effective elastic coefficients at the Γ-point is now realized in a physically meaningful manner. Finite size effects are an important issue in molecular dynamics simulations, particularly for GFMD, they result in a violation of the acoustic sum rule (ASR) for the effective elastic coefficients measured at the Γ-point (ΦΓ). In the previous implementation, the effective elastic coefficients measured at the Γ-point were altered by setting their eigenvalues corresponding to the acoustic modes to zero. This scheme was found to work well for simple Bravais lattices as long as only atoms within the last layer were treated as Green's function atoms. However, it failed to function as expected in all other cases. We therefore adopt a new algorithm to enforce the ASR for ΦΓ (Kong, 2010 [3]) which is implemented in this revision.Summary of revisions: Assuming the lattice under study consists of surface unit cells with n basis atoms labeled by k=1,2,…,n. When all atoms in the lattice are moved by the same amount, i.e., the crystal is rigidly translated, the force on any atom must be zero. This is known as the translational invariance, leading to the so-called acoustic sum rule: k^′_∑Φkα,k^′β(Γ)=0 where Φkα,k^′β(Γ) is the kα,k^′β component of the effective elastic coefficients at the Γ-point; we will denote it as ΦΓ hereafter. α and β enumerate the Cartesian directions. In addition, ΦΓ should be Hermitian (or symmetric, since at the Γ-point, the imaginary part of ΦΓ is zero.) because of the commutative nature of the force constants: . These two properties are expected for ΦΓ, yet the ASR is not satisfied during the measurement (done by FixGFC) because of the finite size effect. A scheme is therefore needed to enforce ASR on ΦΓ afterwards, while the symmetric nature of ΦΓ should be retained.We list below the detailed scheme adopted to enforce ASR implemented in the improved version of GFMD together with some other revisions to the code after the previous release.

• 1. 

  In FixGFMD, the previously employed method to rescale the effective elastic coefficients at Γ is obsoleted. Instead, an iterative procedure is adopted to enforce the acoustic sum rule on ΦΓ (Kong, 2010 [3]). (i)  k′∑Φkα,k′β=0 is enforced by subtracting each element involved by a constant term; this procedure removes the violation of the acoustic sum rule, while in turn, usually destroys the symmetry of the force constant matrix. (ii)  Symmetry is restored by replacing Φkα,k′β and Φk′β,kα with their average value; this will ensure the symmetry of the matrix, however, it will break the acoustic sum rule slightly. (iii)  The above steps are repeated for several iterations, followed by a “symmetric ASR”: similar to step (i), k′∑Φkα,k′β=0 is enforced but only elements with k′?k are subtracted by a constant value, while setting Φk′β,kα=Φkα,k′β.

• 2. 

  In FixGFC, the surface lattice vectors and the relative positions of each atom in the surface unit cell are also computed and written to the binary file, which can be used in FixGFMD to set the equilibrium positions in the Green's function slab based on their lattice indices.

• 3. 

  In FixGFMD, it is now possible to output the total forces applied on atoms in the Green's function slab before applying the elastic forces as a thermal quantity for LAMMPS (http://lammps.sandia.gov [6]). It is also possible to reset these forces to zero before applying the elastic forces.

• 4. 

  In both FixGFC and FixGFMD, the dependence on MPI-enabled FFTW 2.1.5 was lifted. The Fourier transformations are now accomplished by calling the FFT3d wrapper from standard package “kspace” of LAMMPS (Plimpton, 1995; Plimpton et al., 1997; http://lammps.sandia.gov [4,5,6]).

Restrictions: By adopting the new method to enforce the acoustic sum rule, the restriction that atoms in the Green's function slab must be in the same layer is lifted, while it is still necessary to ensure that the mean equilibrium positions of atoms in the Green's function slab satisfy the Born–von Karman boundary condition. In addition, only deformations within the harmonic regime are expected in the slab during Green's function molecular dynamics simulations.Additional comments: The new version is not compatible with the previous one: the contents in the binary file are different and therefore the effective elastic coefficients measured by the previous version of FixGFC cannot be used by the current version of FixGFMD.Running time: FixGFC varies from minutes to days, like a typical molecular dynamics simulation, depending on the system size, the number of processors used, and the complexity of the force field. FixGFMD varies from seconds to hours, depending on the system size and the number of processors used.References:

• [1] 

  L.T. Kong, G. Bartels, C. Campañá, C. Denniston, M.H. Müser, Implementation of Green's function molecular dynamics: An extension to LAMMPS, Computer Physics Communications 180 (6) (2009) 1004–1010.

• [2] 

  C. Campañá, M.H. Müser, Practical Green's function approach to the simulation of elastic semi-infinite solids, Physical Review B (Condensed Matter and Materials Physics) 74 (7) (2006) 075420.

• [3] 

  L.T. Kong, Phonon dispersion measured directly from molecular dynamics simulations, Computer Physics Communications (2010), submitted for publication.

• [4] 

  S.J. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comp. Phys. 117 (1995) 1–19.

• [5] 

  S.J. Plimpton, R. Pollock, M. Stevens, Particle-mesh Ewald and RRESPA for parallel molecular dynamics simulation, in: Proc. of the Eighth SIAM Conference on Parallel Processing for Scientific Computing, Minneapolis, MN, 1997.

• [6] 

  Large-scale Atomic/Molecular Massively Parallel Simulator, LAMMPS, available at: http://lammps.sandia.gov.