Molecular Dynamics Simulation Code Research Articles

Dust Condensation of SiC, SiO in Asymptotic Giant Branch Stellar Winds-SiC Spectrum

We have chosen the Large Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code to calculate the coalescence of silicon carbide (SiC), silicon oxide dust (SiO) in the AGB stellar wind. LAMMPS is a classical molecular dynamics simulation code. At the same time, we consider the effect of temperature on the evolution of molecular dynamics. We also calculated the temperature change of non-spherical SiC, SiO dust coalescence. The condensation temperature range of SiC dust in the AGB stellar wind is [300–500]k and [900–1100]k for SiO. Finally, the infrared spectrum of SiC was calculated using Gaussian 16 software. The 77SiC, 70Si3C3, and 121Si3C3 models have clear characteristic peaks of infrared spectra responding at 5, 8.6, 11.3, 15, 19, and 37 μm.

Efficient end-to-end simulation of time-dependent coherent X-ray scattering experiments.

Physical optics simulations for beamlines and experiments allow users to test experiment feasibility and optimize beamline settings ahead of beam time in order to optimize valuable beam time at synchrotron light sources like NSLS-II. Further, such simulations also help to develop and test experimental data processing methods and software in advance. The Synchrotron Radiation Workshop (SRW) software package supports such complex simulations. We demonstrate how recent developments in SRW significantly improve the efficiency of physical optics simulations, such as end-to-end simulations of time-dependent X-ray photon correlation spectroscopy experiments with partially coherent undulator radiation (UR). The molecular dynamics simulation code LAMMPS was chosen to model the sample: a solution of silica nanoparticles in water at room temperature. Real-space distributions of nanoparticles produced by LAMMPS were imported into SRW and used to simulate scattering patterns of partially coherent hard X-ray UR from such a sample at the detector. The partially coherent UR illuminating the sample can be represented by a set of orthogonal coherent modes obtained by simulation of emission and propagation of this radiation through the coherent hard X-ray (CHX) scattering beamline followed by a coherent-mode decomposition. GPU acceleration is added for several key functions of SRW used in propagation from sample to detector, further improving the speed of the calculations. The accuracy of this simulation is benchmarked by comparison with experimental data.

Unexpected enhancement of sulfuric acid-driven new particle formation by alcoholic amines: The role of ion-induced nucleation

New particle formation (NPF) contributes more than half of the global aerosol. Diethanolamine (DEA) and methyldiethanolamine (MDEA) are the most common amines used to remove CO2 and H2S, which are lost to the atmosphere from CO2 chemical absorbers, livestock and consumer products and are involved in sulfuric acid (SA)-driven NPF. Ion-induced nucleation (IIN) is an important nucleation pathway for NPF. We investigated the role of IIN on DEA and MDEA enhancing SA-driven NPF using density functional method (DFT), molecular dynamics (MD) simulation and atmospheric cluster dynamics code (ACDC). The effects of SO42−, H3O+, NH4+, HSO4−, NO3−, ammonia, methylamine, dimethylamine, trimethylamine and water (W) on the nucleation of SA-DEA were further investigated. The enhancement ability of DEA is greater than that of dimethylamine (DMA) and MDEA. Participation in SA-based NPF is a removal pathway for DEA and MDEA. DEA-SA clusters are generated that not only aggregate DEA and SA molecules, but also increase further growth of atmospheric ions. The very low Gibbs formation free energy highlights the importance of ion-induced nucleation for SA-based NPF. The order of the ability of common atmospheric ions to increase the (SA)(DEA) cluster nucleation is SO42− > H3O+ > NH4+ > HSO4− > NO3−. The addition of 20 water molecules increases the (SA)(DEA)9 cluster from 1.882 nm to 2.053 nm, promoting SA-based NPF. The atmospheric ions accelerate the aggregation rate of the (SA)5(DEA)5 cluster within 15 ns?

Understanding the Effect of Trace Solvent Content on Properties of Polymer Electrolytes through Molecular Dynamics Simulations

The rapid growth of mobile, portable, wearable and flexible electronics leads to the increasing demand for energy storage devices using solid-state polymer electrolytes (PEs), which outperform liquid electrolytes in terms of safety, mechanical properties, and simplicity of device fabrication and packaging. However, processing PEs will always introduce solvent molecules that greatly affect the ionic conductivity and mechanical properties. For example, PEs prepared through solution-casting methods always have solvent residues. A trace amount of water molecules absorbed from the air is also inevitable. Recently, we demonstrated the controlled introduction of solvent molecules to PEs to balance the ionic conductivity and mechanical stiffness for structural energy storage applications. To better understand how solvent molecules behave and interact with other components in PEs, here we present the molecular dynamics simulation of a representative polymer electrolyte system with various water content. We use simulation results to determine the effect of trace water content before forming a liquid phase on ionic conductivity and mechanical properties. The insights into the molecular interactions in the PE system will help us design and optimize Pes’ composition and processing for practical applications.The simulation model of polymer electrolyte is built with polyethylene oxide (PEO) and lithium perchlorate (LiClO4) with various water contents, in which the water molecule to lithium-ion ratio ranges from 0 to 3. The electrolyte with each water content is simulated between two graphene electrodes to determine its ionic conductivity. Uniaxial deformation has been performed on the electrolyte to obtain the mechanical properties. All simulations were performed using the molecular dynamics simulation code LAMMPS with the CHARMM force field.The results show that the ionic conductivity of the polymer electrolyte system increases significantly (up to one order of magnitude) with the increase of water content (up to 3 water molecules per lithium ion), even when the added water does not form a continuous liquid phase. The change of ionic conductivity with water content is correlated to the degree of association between different types of ions or molecules in the system, as evidenced by the evaluation of the radial distribution functions. As the association between polymer molecules and lithium ions reduces with increasing water, it becomes easier for the lithium ions to diffuse and resulting in higher ionic conductivity. It is also observed that the perchlorate ions’ interactions with polymer molecules remain the same with different water contents, which shows different roles of lithium ions and perchlorate ions in ion conduction in this system. On the other hand, the modulus of elasticity of the polymer electrolyte does not change much with the increase of water, which agrees with the previous experimental work of our group. This means that the trace amount of water is strongly associated with other solid molecules or ions and is not affecting the stiffness of the system as long as no liquid phase is formed. The results will lead to novel strategies to design polymer electrolytes with both high ionic conductivity and good mechanical properties for flexible or multifunctional energy storage applications.

Electronic structure of lattice relaxed alternating twist tNG-multilayer graphene: from few layers to bulk AT-graphite

Alternating twist (AT) multilayer graphene systems are at the heart of recent research efforts on flat band superconductivity and therefore precise descriptions of their atomic and electronic structures are desirable. We present the electronic structure of AAʹAAʹ… stacked AT N-layer (tNG) graphene for 3-10, 20 layers and bulk AT graphite systems where the atomic structure is relaxed using a molecular dynamics simulation code. The low energy bands depend sensitively on the relative sliding between the layers but we show explicitly up to N = 6 that the highly symmetric AA′AA′. stacking is energetically preferred among all interlayer sliding geometries of each added layer, justifying why experimental devices consistently show results compatible with this geometry. It is found that lattice relaxation enhances electron–hole asymmetry, and leads to small reductions of the magic angle values with respect to analytical or continuum model calculations with fixed tunneling strengths that we quantify from few layers to bulk AT-graphite. The twist angle error tolerance near the magic angles obtained by maximizing the density of states of the nearly flat bands expand progressively from 0.05∘ for twisted bilayer graphene to up to 0.2∘ for AT-graphite, hence allowing a greater twist angle flexibility in multilayers. We further comment on the role of perpendicular electric and magnetic fields in modifying the electronic structure of the system and how the decoupling of tNG multilayers bands allows mapping onto those of periodic AT-graphite’s at different k z values.

Kajian Simulasi Dinamika Molekul Adsorpsi Hidrogen pada Carbone Nanotube dengan Variasi Chirality dan Temperatur Menggunakan Kode LAMMPS

Hydrogen adsorption has been simulated on carbon nanotubes for optimum hydrogen absorption. Parameters that affect the amount of hydrogen absorbed have been studied, such as the effect of chirality and temperature on hydrogen absorption in CNTs. The simulation method of hydrogen adsorption on carbone nanotubes uses molecular dynamics simulation code LAMMPS, applies Lennard-Jones interatomic potential and hydrogen atom movement using Van Der Waals force with Microcanonical Ensemble. Data analysis is the output of LAMPS in the form of data in XYZ format. The data contains information in the form of integration steps, number of atoms, temperature, pressure, potential energy, kinetic energy, volume, van der Waals energy, total simulation time and hydrogen absorption. The simulation results show that the optimum absorption occurs at run 10000 and a temperature of 100 K, for armchair chirality of 10 atoms, chirality of 12 atoms and zigzag chrality of 5 atoms. Formation of hydrogen coordinates with Avogadro software, formation of CNT coordinates with VMD software and visualization of hydrogen adsorption on CNTs using VMD software.

FE-CLIP: A tool for the calculation of the solid–liquid interfacial free energy

In manufacturing industries, predicting the work of adhesion between complex solid and liquid surfaces has become essential. FE-CLIP offers a routine for evaluating the work of adhesion between solid and liquid surfaces by calling a subroutine from a molecular dynamics simulation code. When FE-CLIP is applied to the solid–liquid interface, liquid molecules are separated from the solid surface according to its shape using a set of spherical potentials. This is efficient when applied to the solid surfaces with complex structures such as polymer-grafted surfaces. An adaptive scheme for updating the parameters contained in the potential and automatic refinement of integration points are introduced to facilitate the application of FE-CLIP to various solid–liquid interfaces. We applied FE-CLIP to the separation of water from a polymer-grafted gold surface to demonstrate that the proposed method gives reliable results by suppressing the variation of the free energy gradient, which is important for accurate numerical integration. Program summaryProgram Title: FE-CLIPProgram Files doi:http://dx.doi.org/10.17632/7gmcvftfwj.1Licensing provisions: GPLv3Programming language: FortranNature of problem: The adhesion free energy of liquid–solid interfaces is one of the important properties for their industrial application. FE-CLIP enables quantitative evaluation of the strength of adhesion between solid and liquid surfaces using a molecular dynamics simulation. Unlike other existing methods, FE-CLIP is efficiently applicable to complex solid surfaces such as polymer-grafted surfaces.Solution method: The work of adhesion between solid and liquid surfaces is calculated using thermodynamic integration. The integration parameters of thermodynamic integration are the diameter and depth of artificial potentials of the Lennard-Jones-type, whose variation contributes to separating the liquid molecules from a solid surface. By putting a set of spherical potentials, the liquid on the solid surface is separated according to the surface shape. In addition, the parameters of the potentials are automatically updated so as to suppress the variation of the free energy gradient.

Atomistic simulations of molten carbonates: Thermodynamic and transport properties of the Li2CO3-Na2CO3-K2CO3 system.

Although molten carbonates only represent, at most, a very minor phase in the Earth's mantle, they are thought to be implied in anomalous high-conductivity zones in its upper part (70-350 km). Besides, the high electrical conductivity of these molten salts is also exploitable in fuel cells. Here, we report quantitative calculations of their properties, over a large range of thermodynamic conditions and chemical compositions, which are a requisite to develop technological devices and to provide a better understanding of a number of geochemical processes. To model molten carbonates by atomistic simulations, we have developed an optimized classical force field based on experimental data of the literature and on the liquid structure issued from ab initio molecular dynamics simulations performed by ourselves. In implementing this force field into a molecular dynamics simulation code, we have evaluated the thermodynamics (equation of state and surface tension), the microscopic liquid structure and the transport properties (diffusion coefficients, electrical conductivity, and viscosity) of molten alkali carbonates (Li2CO3, Na2CO3, K2CO3, and some of their binary and ternary mixtures) from the melting point up to the thermodynamic conditions prevailing in the Earth's upper mantle (∼1100-2100 K, 0-15 GPa). Our results are in very good agreement with the data available in the literature. To our knowledge, a reliable molecular model for molten alkali carbonates covering such a large domain of thermodynamic conditions, chemical compositions, and physicochemical properties has never been published yet.

Multiresolution molecular mechanics: Implementation and efficiency

Atomistic/continuum coupling methods combine accurate atomistic methods and efficient continuum methods to simulate the behavior of highly ordered crystalline systems. Coupled methods utilize the advantages of both approaches to simulate systems at a lower computational cost, while retaining the accuracy associated with atomistic methods. Many concurrent atomistic/continuum coupling methods have been proposed in the past; however, their true computational efficiency has not been demonstrated. The present work presents an efficient implementation of a concurrent coupling method called the Multiresolution Molecular Mechanics (MMM) for serial, parallel, and adaptive analysis. First, we present the features of the software implemented along with the associated technologies. The scalability of the software implementation is demonstrated, and the competing effects of multiscale modeling and parallelization are discussed. Then, the algorithms contributing to the efficiency of the software are presented. These include algorithms for eliminating latent ghost atoms from calculations and measurement-based dynamic balancing of parallel workload. The efficiency improvements made by these algorithms are demonstrated by benchmark tests. The efficiency of the software is found to be on par with LAMMPS, a state-of-the-art Molecular Dynamics (MD) simulation code, when performing full atomistic simulations. Speed-up of the MMM method is shown to be directly proportional to the reduction of the number of the atoms visited in force computation. Finally, an adaptive MMM analysis on a nanoindentation problem, containing over a million atoms, is performed, yielding an improvement of 6.3–8.5 times in efficiency, over the full atomistic MD method. For the first time, the efficiency of a concurrent atomistic/continuum coupling method is comprehensively investigated and demonstrated.

Coherent dynamic structure factors of strongly coupled plasmas: A generalized hydrodynamic approach

A generalized hydrodynamic fluctuation model is proposed to simplify the calculation of the dynamic structure factor S(ω, k) of non-ideal plasmas using the fluctuation-dissipation theorem. In this model, the kinetic and correlation effects are both included in hydrodynamic coefficients, which are considered as functions of the coupling strength (Γ) and collision parameter (kλei), where λei is the electron-ion mean free path. A particle-particle particle-mesh molecular dynamics simulation code is also developed to simulate the dynamic structure factors, which are used to benchmark the calculation of our model. A good agreement between the two different approaches confirms the reliability of our model.

Ion dynamics in a Paul trap driven by various radio frequency waveforms

In this study we explore single ion dynamics in a Paul trap (three-dimensional radio frequency (RF) quadrupole ion trap – 3D QIT), driven by Sinusoidal, Rectangular, Sawtooth and Triangular RF waveforms. A molecular dynamic simulation code in Python has been written to explore ion dynamics in the trap. It is shown that various radio frequency waveforms produce various shift in ion oscillation in the trap. Furthermore phase space plots and the magnitude of the potential well depth for these various waveforms are given. Possible influences of these various RF waveforms have been summarized.

Revealing the complex conduction heat transfer mechanism of nanofluids

Nanofluids are two-phase mixtures consisting of small percentages of nanoparticles (sub 1–10 %vol) inside a carrier fluid. The typical size of nanoparticles is less than 100 nm. These fluids have been exhibiting experimentally a significant increase of thermal performance compared to the corresponding carrier fluids, which cannot be explained using the classical thermodynamic theory. This study deciphers the thermal heat transfer mechanism for the conductive heat transfer mode via a molecular dynamics simulation code. The current findings are the first of their kind and conflict with the proposed theories for heat transfer propagation through micron-sized slurries and pure matter. The authors provide evidence of a complex new type of heat transfer mechanism, which explains the observed abnormal heat transfer augmentation. The new mechanism appears to unite a number of popular speculations for the thermal heat transfer mechanism employed by nanofluids as predicted by the majority of the researchers of the field into a single one. The constituents of the increased diffusivity of the nanoparticle can be attributed to mismatching of the local temperature profiles between parts of the surface of the solid and the fluid resulting in increased local thermophoretic effects. These effects affect the region surrounding the solid manifesting interfacial layer phenomena (Kapitza resistance). In this region, the activity of the fluid and the interactions between the fluid and the nanoparticle are elevated. Isotropic increased nanoparticle mobility is manifested as enhanced Brownian motion and diffusion effects

Molecular Dynamic Simulations of a Simplified Nanofluid

This study describes the methodology that was developed to run a Molecular Dynamics Simulation (MDS) code to simulate the behaviour of a single nanoparticle dispersing in a fluid with a temperature gradient. A soft disk model described by the Lennard-Jones potential is used to simulate the system. The nanoparticle is assembled via the use of four subdomains of interatomic interactions and hence presents in full resolution the transfer of energy from the fluid-to-solid-to-fluid subdomains. A cluster computing system (HTCondor) was used to perform a large scale deployment of the MDS code. The obtained showcase results were successfully evaluated using three widely documented tests from the associated literature (Randomness, Radial Distribution and Velocity Autocorrelation Distribution Functions). It was discovered that the nanoparticle travels a larger distance in the fluid than the distance travelled by a fluid molecule (recovery region). The findings were confirmed by calculating the Green-Kubo self-diffusivity coefficient halfway through the simulation at which an enhancement of 156% was discovered in favour of the Nanoparticle. This might be the physical mechanism responsible for the experimentally observed thermal performance enhancement in nanofluids.

Ls1 mardyn: The Massively Parallel Molecular Dynamics Code for Large Systems.

The molecular dynamics simulation code ls1 mardyn is presented. It is a highly scalable code, optimized for massively parallel execution on supercomputing architectures and currently holds the world record for the largest molecular simulation with over four trillion particles. It enables the application of pair potentials to length and time scales that were previously out of scope for molecular dynamics simulation. With an efficient dynamic load balancing scheme, it delivers high scalability even for challenging heterogeneous configurations. Presently, multicenter rigid potential models based on Lennard-Jones sites, point charges, and higher-order polarities are supported. Due to its modular design, ls1 mardyn can be extended to new physical models, methods, and algorithms, allowing future users to tailor it to suit their respective needs. Possible applications include scenarios with complex geometries, such as fluids at interfaces, as well as nonequilibrium molecular dynamics simulation of heat and mass transfer.

Open Source Software to Visualize Complex Data on Remote CEA's Supercomputing Facilities

In order to guaranty performances of complex systems using numerical simulation, CEA is performing advanced data analysis and scientific visualization with open source software using High Performance Computing (HPC) capability. The diversity of the physics to study produces results of growing complexity in terms of large-scale, high dimensional and multivariate data. Moreover, the HPC approach introduces another layer of complexity by allowing computation amongst thousands of remote cores accessed from sites located hundreds of kilometers away from the computing facility. This paper presents how CEA deploys and contributes to open source software to enable production class visualization tools in a high performance computing context. Among several open source projects used at CEA, this presentation will focus on Visit, VTK and Paraview. In the first part we will address specific issues encountered when deploying VisIt and Paraview in a multi-site supercomputing facility for end-users. Several examples will be given on how such tools can be adapted to take advantage of a parallel setting to explore large multi-block dataset or perform remote visualization on material interface reconstructions of billions of cells. Then, the specific challenges faced to deliver Paraview’s Catalyst capabilities to end-users will be discussed. In the second part, we will describe how CEA contributes to open source visualization software and associated software development strategy by emphasizing on two recent development projects. The first is an integrated simulation workbench providing plugins for every step required to achieve numerical simulation independently on a local or a remote computer. Embedded in an Eclipse RCP environment, VTK views allow the users to perform data input using interaction or mesh preview before running the simulation code. Contributions to VTK have been made in order to smoothly integrate these technologies. The second details how recent developments at CEA have helped to visualize and to analyze results from ExaStamp, a parallel molecular dynamics simulation code dealing with molecular systems ranging from a few millions up to a billion atoms. These developments include a GPU intensive rendering method specialized for atoms and specific parallel algorithms to process molecular data sets.

Parallel implementation of three-dimensional molecular dynamic simulation for laser-cluster interaction

The objective of this article is to report the parallel implementation of the 3D molecular dynamic simulation code for laser-cluster interactions. The benchmarking of the code has been done by comparing the simulation results with some of the experiments reported in the literature. Scaling laws for the computational time is established by varying the number of processor cores and number of macroparticles used. The capabilities of the code are highlighted by implementing various diagnostic tools. To study the dynamics of the laser-cluster interactions, the executable version of the code is available from the author.

Potential for improvement in high heat flux HyperVapotron element performance using nanofluids

HyperVapotron (HV) elements have been used extensively as high heat flux beam stopping components in nuclear fusion research facilities. These water-cooled heat exchangers use a boiling heat transfer mechanism and so are inherently limited by their critical heat flux (CHF). The use of a nanofluid as the coolant, instead of water, promises to enhance the heat transfer performance of the HV and increase the CHF by a factor of 2 or 3, which would lead to a step-change improvement in the power handling capability. This paper reports on computational and experimental analyses which have indicated mechanisms for the enhanced thermal performance of nanofluids. A molecular dynamics simulation code has been developed which has identified heat transfer augmentation mechanisms that depart from classical thermodynamics associated with the presence of nanoparticles. In addition, an experiment has been conducted which uses particle image velocimetry to measure the flow field in a full-scale HV. Past studies have yielded qualitative experimental results, but the measurements reported here provide quantitative data to aid the understanding of the initial flow field inside the HV (i.e., before a heat flux is applied). Further, the experiment is conducted using both water and Al2O3–water nanofluid as the flow medium. Thus, these velocity measurements offer a first indication for potentially enhanced heat transfer in HV devices when nanofluids are used as the coolant. The improved understanding of the HV flow regime and the cooling advantage of nanofluids could assist the design of advanced high heat flux components for future fusion machines.

Numerical Studies of the Squeeze-Film Damping of MEMS Devices with Perforations in the Non-Continuum Regime

Predicting squeeze-film air damping of resonators in rare air is crucial in the design of high-Q devices for various applications. In the past, there have been two approaches to treat the squeeze-film air damping in non-continuum regime: using effective viscosity coefficient and using the molecular dynamics method. And most of the previous work focused on devices in which the rarefaction effects of air are not significant. For such cases, continuum theory is often adequate. However, we have investigated the air damping on oscillating structures in the free molecular regime in which classical continuum theory is no longer valid. Based on this premise, Hutcherson (2004 J. Micromesh. Microeng. 14 1726-1733) has developed a molecular dynamics simulation code and used in predicting quality factors of an oscillating micro-plate at low pressures. However, his work is valid only for non-perforated micro-plate. This paper, a brief description of the molecular dynamics method is presented first. Then a molecular dynamics simulation code has been developed and used in predicting quality factors of a perforated oscillating micro-plate in free molecular regime. And we have found that the molecular dynamics simulation results have shown an excellent agreement with the experimental data of Kwok et al. Finally, the limitations of the present molecular dynamics simulation code have been reported.

Micellization Studied by GPU-Accelerated Coarse-Grained Molecular Dynamics.

The computational design of advanced materials based on surfactant self-assembly without ever stepping foot in the laboratory is an important goal, but there are significant barriers to this approach, because of the limited spatial and temporal scales accessible by computer simulations. In this paper, we report our work to bridge the gap between laboratory and computational time scales by implementing the coarse-grained (CG) force field previously reported by Shinoda et al. [Shinoda, W.; DeVane, R.; Klein, M. L. Mol. Simul. 2007, 33, 27-36] into the HOOMD-Blue graphical processing unit (GPU)-accelerated molecular dynamics (MD) software package previously reported by Anderson et al. [Anderson, J. A.; Lorenz, C. D.; Travesset, A. J. Comput. Phys. 2008, 227, 5342-5359]. For a system of 25 750 particles, this implementation provides performance on a single GPU, which is superior to that of a widely used parallel MD simulation code running on an optimally sized CPU-based cluster. Using our GPU setup, we have collected 0.6 ms of MD trajectory data for aqueous solutions of 7 different nonionic polyethylene glycol (PEG) surfactants, with most of the systems studied representing ∼1 000 000 atoms. From this data, we calculated various properties as a function of the length of the hydrophobic tails and PEG head groups. Specifically, we determined critical micelle concentrations (CMCs), which are in good agreement with experimental data, and characterized the size and shape of micelles. However, even with the microsecond trajectories employed in this study, we observed that the micelles composed of relatively hydrophobic surfactants are continuing to grow at the end of our simulations. This suggests that the final micelle size distributions of these systems are strongly dependent on initial conditions and that either longer simulations or advanced sampling techniques are needed to properly sample their equilibrium distributions. Nonetheless, the combination of coarse-grained modeling and GPU acceleration marks a significant step toward the computational prediction of the thermodynamic properties of slowly evolving surfactant systems.

Phonon dispersion measured directly from molecular dynamics simulations

A method to measure the phonon dispersion of a crystal based on molecular dynamics simulation is proposed and implemented as an extension to an open source classical molecular dynamics simulation code LAMMPS. In the proposed method, the dynamical matrix is constructed by observing the displacements of atoms during molecular dynamics simulation, making use of the fluctuation–dissipation theory. The dynamical matrix can then be employed to compute the phonon spectra by evaluating its eigenvalues. It is found that the proposed method is capable of yielding the phonon dispersion accurately, while taking into account the anharmonic effect on phonons simultaneously. The implementation is done in the style of fix of LAMMPS, which is designed to run in parallel and to exploit the functions provided by LAMMPS; the measured dynamical matrices could be passed to an auxiliary postprocessing code to evaluate the phonons. Program summaryProgram title: FixPhonon, version 1.0Catalogue identifier: AEJB_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEJB_v1_0.htmlProgram obtainable from: CPC Program Library, Queenʼs University, Belfast, N. IrelandLicensing provisions: GNU General Public licenseNo. of lines in distributed program, including test data, etc.: 105 393No. of bytes in distributed program, including test data, etc.: 3 231 800Distribution format: tar.gzProgramming language: C++Computer: AllOperating system: LinuxHas the code been vectorized or parallelized?: Yes. 1 to N processors may be usedRAM: Depends on problem, ≈1 kB to several MBClassification: 7.8External routines: MPI, FFT, LAMMPS version 15, January 2010 (http://lammps.sandia.gov/)Nature of problem: Atoms in solids make ceaseless vibrations about their equilibrium positions, and a collective vibration forms a wave of allowed wavelength and amplitude. The quantum of such lattice vibration is called the phonon, and the so-called “lattice dynamics” is the field of study to find the normal modes of these vibrations. In other words, lattice dynamics examines the relationship between the frequencies of phonons and the wave vectors, i.e., the phonon dispersion. The evaluation of the phonon dispersion requires the construction of the dynamical matrix. In atomic scale modeling, the dynamical matrices are usually constructed by deriving the derivatives of the force field employed, which cannot account for the effect of temperature on phonons, with an exception of the tedious “quasi-harmonic” procedure.Solution method: We propose here a method to construct the dynamical matrix directly from molecular dynamics simulations, simply by observing the displacements of atoms in the system thus making the constructing of the dynamical matrix a straightforward task. Moreover, the anharmonic effect was taken into account in molecular dynamics simulations naturally, the resultant phonons therefore reflect the finite temperature effect simultaneously.Restrictions: A well defined lattice is necessary to employ the proposed method as well as the implemented code to evaluate the phonon dispersion. In other words, the system under study should be in solid state where atoms vibrate about their equilibrium positions. Besides, no drifting of the lattice is expected. The method is best suited for periodic systems, although non-periodic system with a supercell approach is also possible, it will however become inefficient when the unit cell contains too many atoms.Additional comments: The readers are encouraged to visit http://code.google.com/p/fix-phonon for subsequent update of the code as well as the associated postprocessing code, so as to keep up with the latest version of LAMMPS.Running time: Running time depends on the system size, the numbers of processors used, and the complexity of the force field, like a typical molecular dynamics simulation. For the third example shown in this paper, it took about 2.5 hours on an Intel Xeon X3220 architecture (2.4G, quadcore).