introduction to NVE, NVT, NpT ensembles

The first term with which to get acquainted is “ensemble”. In the context of statistical mechanics, an ensemble is a collection of configurations of a system. These configurations will have a set of common constraints such as temperature, volume and pressure, but will differ in the positions and velocities of the individual components (atom/molecule/fragment) across the different configurations. So, even though these configurations differ microscopically, they possess the same common macroscopic properties (based on the ensemble they belong to). This means that ensemble averages for these properties can be computed. The simplest such ensemble is the microcanonical ensemble, in which the constraints are number of particles, volume and total energy, i.e., these quantities are constrained to be the same across all the configurations belonging to the microcanonical ensemble, and hence this ensemble is also known as the “NVE” ensemble. The common nomenclature of the ensemble provides information about the constraints in that ensemble. Hence, the NVE ensemble contains all possible combinations of positions and velocities for the particles which satisfy an energy constraint, provided all these configurations have the same particle density. The total number of configurations for an ensemble is given by its “partition function”, and for the microcanonical ensemble, the microcanonical partition function (also called density of states) is represented as Ω(N,V,E).
The fundamental importance of the density of states is based on the relationship
S(N,V,E) = kB ln Ω(N,V,E) (2.2)
where S is the entropy and kB is Boltzmann’s constant. Eq. 2.2 shows how a macroscopic property like entropy is directly related to the density of states for a system at constant density. A more convenient and commonly used ensemble is the canonical (NVT) ensemble, where the density and temperature are constrained across the configuration.
Similar to the microcanonical ensemble, a thermodynamic quantity can be related to the partition function. In this case, it is the Helmholtz free energy A = −kBT ln Q(N,V,T) (2.6) A commonly used ensemble in MC is the grand canonical ensemble, in which the chemical potential, volume and temperature are constrained. This means that the number of particles is allowed to fluctuate to satisfy the chemical potential constraint. The NpT ensemble, also called the isothermal-isobaric ensemble is also a commonly used ensemble. Here, the number of particles, temperature and pressure are constrained.
The fundamental relationship for the NpT ensemble is the expression relating the Gibbs free energy to the partition function: G = −kBT ln ∆(N,P,T) (2.8) where G is the Gibbs free energy, a thermodynamic quantity of fundamental importance to the problem of phase equilibrium. In all the ensembles listed above, the partition function integrals grow rapidly with increasing complexity and system sizes. Hence, trying to solve these integrals analytically is impractical. Molecular simulations aid in sampling the configuration space for these ensembles subject to the probability distribution for a given ensemble, and ensemble averages can be obtained using expressions similar to Eq. 2.5.


Ref: http://etd.nd.edu/ETD-db/theses/available/etd-08182009-110710/unrestricted/JayaramanS082009D.pdf

Molecular Dynamic simulation of heating of gold nanorod in water (nature 2015)

Molecular Dynamic simulation of heating of gold nanorod in water

Potential Model for Gold

Molecular dynamics (MD) simulations of the gold nano rods in presence of water were
performed using embedded atom method (EAM)
[7. Daw, M. & Baskes, M. Semiempirical, Quantum Mechanical Calculation of Hydrogen Embrittlement in Metals. Phys. Rev. Lett. 50, 1285–1288 (1983).
8. Daw, M. & Baskes, M. Embedded-atom method: Derivation and application to impurities, surfaces,and other defects in metals. Phys. Rev. B 29, 6443–6453 (1984). ] .
For the EAM, the total energy (Etot) for a system of N atoms can be written as equa S14.

To model the interactions of water molecules, we have employed flexible SPC/Fw water
model to conduct simulations at various temperatures9
. The potential interaction in the flexible
SPC/Fw, as shown in equation (S15), is a sum of pair-wise interactions for bonded and nonbonded
terms.

Potential model and parameters for gold-water interaction
We used the Lorentz-Berthelot rule to estimate the interaction between water and gold
that is modeled using standard 6-12 potential.

Melting temperatures as a function of nanorod aspect ratios:
We have further calculated the bulk melting point as well as the melting of nanorods with
different aspect ratios. The melting point of bulk gold predicted by the EAM potential used in this work is 1281 K. The experimental value is 1337 K and thus the prediction of our potential model is very good.
The melting points for nanorods with different aspect ratios (AR) are summarized below. It can be seen that nanorods with higher aspect ratios have lower melting points compared to bulk. This is not surprising since the surface area (Table S3) increases with increasing AR.


Ref: http://www.nature.com/srep/2015/150130/srep08146/extref/srep08146-s1.pdf

Gold and Lennard-Jones potential

Gold
e/k=5123°K
sigma=2.637 Angström
[T. HALICIOGLU and G. M. POUND: Calculation of Potential Energy Parameters ; phys. stat. sol. (a) 30, 619 (1975) ]

---------------------------------
Lennard-Jones potential:

A = 4εσ6 and B = 4εσ12 

Using the continuous approach, where the atoms at discrete locations on the mol- ecule are averaged over a surface, the molecular interatomic energy is obtained by calculating integrals over the surfaces of each molecule, given by

 where η1 and η2 represent the mean surface density of atoms on each molecule. Further, we may define the integral In in the form of

and therefore, E = η1η2(−AI3 + BI6).
The Lennard-Jones parameters for gold nanoparticles are taken from the work
of Pu et al. [Q. Pu, Y. Leng, X. Zhao, P.T. Cummings, Molecular simulations of stretching gold nanowires in solvents. Nanotechnology 18, 424007 (2007)], who use 
ε = 0.039kcal mol−1 and σ = 2.934Å 
corresponding to A = 7.2465eV ×Å6 and B = 4,622.6273eV ×Å12. 

From the fact that gold adopts a face-centred cubic (fcc) crystal structure, the triangular arrangement is employed to determine a mean surface density of gold nanoparticles. Assuming that each atom of gold is coordinated with six other atoms, the area of a unit cell is given by A⋆ = √3d2/2Å2 where d denotes the equilibrium spacing which can be obtained as d = 21/6σ = 3.2933Å. Consequently, the mean surface density η of the gold layers is taken to be 1/A⋆ = 0.1065Å−2. 

Ref: http://link.springer.com/article/10.1007%2Fs10910-010-9796-x

example of NVT equilibration in LAMMPS

All the simulation was carried out in the NVT ensemble.
The initial equilibration was done by heating the system from 20 K to 300 K with discrete steps of 20 K/50 ps.
After the equilibration, the systems were run at a temperature of 300 K for additional 5 ns.
Dynamics of transfer, and variance in distributions  were analyzed by heating equilibrated systems from 300 K to 400 K followed by cooling from 400 K to 300 K in discrete steps of 10 K/ 200 ps.
For accurate analyses at each temperature, the last snapshot from the heating cycles were run for more than 2 ns and data collected in the equilibrated phase as seen through block averages of
total energy.
All simulations were carried out using the LAMMPS simulation package.
The equations of motion were integrated using the velocity-Verlet algorithm using a time step
of 1 fs. For NVT ensembles, a Nosé-Hoover thermostat with a damping time of 100 fs
was used to maintain the temperature of the system.

Ref: http://www.rsc.org/suppdata/nr/c3/c3nr05671f/c3nr05671f.pdf

triple point: gold (Au) in both the solid and liquid phases; binodal, spinodal.

The Liquid-vapor critical point in a pressure–temperature phase diagram 

is at the high-temperature extreme of the liquid–gas phase boundary 

(the dotted green line shows the anomalous behavior of water).

The figure shows the schematic PT diagram of a pure substance (as opposed to mixtures, which have additional state variables and richer phase diagrams, discussed below). The commonly known phases solid, liquid and vapor are separated by phase boundaries, i.e. pressure-temperature combinations where two phases can coexist.

 At the triple point, even all three phases coexist. However, the liquid-vapor boundary terminates in an endpoint at some critical temperature Tc and critical pressure pc. This is the critical point.

In water, the critical point occurs at around 647 K and 22.064 MPa.

- Gold 7250 K 510MPa.

At the critical point, only one phase exists. The heat of vaporization is zero. 

----------------------

Van der Walls-like phase-diagram of gold

Temperature 0.5ev<->0.5ev/kB=5802°K

kB=8.6173 10^-5eV/°K

DOI: 10.1016/j.apsusc.2005.03.019 ; [Vitiello,2005]

--------------------------------------------

GOLD

Electron configuration [Xe] 4f14 5d10 6s1

melting point: 1336.15°K

Heat of fusion (10^3 J/kg): 64.5

or 12.55 kJ·mol^−1

boiling point: 3243 °K

Heat of vaporization (10^3 J/kg): 1578

or 342 kJ·mol^−1

Density near r.t.    19.30 g·cm−3

when liquid, at m.p.    17.31 g·cm−3

Molar heat capacity   25.418 J·mol−1·K−1

---------------------------------------------

Problem:

1) Using the CRC thermodynamic data from the handout in class, make a plot of the Gibbs free energy (in [J/mole]) as a function of temperature for gold (Au) in both the solid and liquid phases over the temperature range of 1000-2000K. Both free energies must be on the same plot. Determine the equilibrium temperature between liquid and solid Au from your calculations.

2) Yet another approach to calculating the Gibbs free energy is using empirical polynomial fitting equations from thermodynamic databases. One such example is the CALPHAD database. As an example of the fitting equations I have listed the equations for solid Au in the FCC crystal structure below. Make a plot of the Gibbs free energy (in [J/mole]) as a function of temperature over the range of 1000-2000 [K] using both the CRC data from problem 1 and the CALPHAD data. Both calculations must be on the same plot. Comment as to the agreement or lack thereof between the two calculations.

3) Calculate the triple point (as T, P) for gold. (NOTE that all the thermodynamic information you need to do this calculation can be found from the CRC Handbook data. I'm not looking for the whole phase diagram, just the triple point. There are a variety of ways to approach this problem, but I'm looking for a pretty simple one.)

Solution

1) Using the CRC thermodynamic data from the handout in class, make a plot...

The CRC data for Au gives the following equations for the Gibbs free energy as a function of temperature:

Plugging into EXCEL and plotting between 1000-2000 [K] gives the following plot:

With this calculation, the crossing point is at about T = 1341 [K], which is the equilibrium melting temperature (a bit higher than the stated temperature of 1336.16 [K].)

NOTE: The CRC is pretty confusing by calling the function for Gibbs free energy "F_T-H₂₉₈." In fact, the enthalpy at 298 [K] is included in the equations and fitting parameters. What it means is that the thermodynamic functions are determined relative to a temperature of 298 [K].

2) Yet another approach to calculating the Gibbs free energy...

I used EXCEL again to make the calculation and the plot.

3) Calculate the triple point (as T, P) for gold.

The simple trick approach uses the fact that we discussed in class that the slope of the S-L equilibrium line on the P-T phase diagram is VERY STEEP (a nearly vertical line with pressure) so that the equilibrium melting temperature changes very little with pressure. 

This means that, within a small error, the triple point temperature is given by the equilibrium melting temperature at 1 [atm], which for Au is found in the CRC handout as T_M = 1336.16 [K]. The problem is half done!

To find the triple point pressure, we need to find the equilibrium pressure for the temperature of 1336.16 [K] for either the S-V or the L-V equilibrium. I decided to use the "shortcut" approach (see problem 2 on themidterm for F08 for more information.) The Clausius-Clapeyron equation gives:

In this equation, we can either use the enthalpy of vaporization for ∆H, and the vaporization temperature for T₀ (which would give us the L-V equilibrium), or we could use the enthalpy and temperature of sublimation (which would give us the S-V equilibrium.) The CRC has the data for the vaporization, but not the sublimation information, so...

There's our answer for the triple point of Au: T = 1336 [K], P = 2.44 x 10^-7 [atm].

Ref: http://oregonstate.edu/instruct/me581/Homework/F08/ME581Hmwk3.html

Comparison of time- and length-scales, accessible to different types of simulation techniques (quantum simulations (QM), molecular dynamics (MD), Brownian Dynamics (BD) and hydrodynamics/fluid dynamics (HD))

Computing every force on each atom is a demanding part of a MD simulation. In order to limit the number of force calculations, the smallest possible system – in terms of number of atoms – should be simulated. The minimum system size necessary to compute a given desired property is not always obvious and a common practice is to analyze how increasing system size influences obtained results. A number of algorithmic approaches exist that help minimize system size such as applying periodic boundary conditions.
In addition to how large a system must be, one must determine how long a simulation must be run to obtain reliable answers (i.e. with sufficient accuracy and precision). Related to this is that the MD simulation time step should be small compared to the period associated with the highest frequency atomic event. A typical time step for MD simulations of metals is of order 1 fs (i.e. 1 E-15 s). The very small time step is the main drawback of MD method compared to other commonly used numerical simulations. Processes with time scale of order microseconds and larger are outside the reasonable realm for even highly parallel, efficient MD codes. Conversely, MC techniques cannot probe dynamics but they can probe long time equilibrium properties.

Figure 2.1 Schematic comparison of time- and length-scales, accessible to different types of simulation techniques (quantum simulations (QM), molecular dynamics (MD), Brownian Dynamics (BD) and hydrodynamics/fluid dynamics (HD))

Accessible length- and time- scales of microscopic simulation should be considered. Figure 2.1 shows a schematic representation for different types of simulations in a length-time-diagram. It is clear that the more detailed a simulation technique addresses, the smaller is the accessibility of times and length scales. Hence, quantum simulations (QS), where fast motion of electrons are considered, are located in the lower left corner of the diagram and typical length time scales are of order of Å (Å) and ps. MD method puts fixed partial charges on interaction sites or adds a dynamics model for polarization effects to simulate electronic distribution.  Hence, the time scale of the system is determined by the time of collision between atoms and atomic vibration frequencies. Consequently, the accessible time scale is of order ns and corresponding accessible length scale becomes of order 10-10000 Å. Note Figure 2.1 was originally presented in 2002; accessible length scales for MD have increased dramatically. Brownian Dynamics (BD) can be applied to deal with particles in a solvent medium, where one is not interested in a detailed description of the solvent. In BD the effect of the solvent is hidden in average quantities. If one is not interested in a molecular picture but in macroscopic quantities, the concepts of hydrodynamics (HD) may be applied, where the system properties are hidden in effective numbers, like density and sound velocity.

introduction to Molecular dynamics, part 2: The MD Algorithm; why ∆t equal to 1 femtosecond?

2.2 The MD Algorithm

2.2.1 Verlet’s Algorithm and Velocity Verlet’s Algorithm

The governing equations of motions are all ordinary differential equations and a
standard method for the solutions is the finite difference method. Of a large number of
finite-difference approaches devised, the most commonly used one is the Runge-Kutta
(RK) method. [33,34] The beauty of the RK method is that the order of a finite
difference method can often be increased by using positions and velocities from
several points in time instead of only the current time. However, if a fourth-order RK
method is applied to MD simulation, four interatomic force evaluations per atom
would be required per step, which will make the simulation too slow. Thus, reduction
of order is desired; especially, the time expense of calculating intermolecular force
more than once per atom per step should be also avoided. However, the algorithm
must be stable.
The most widely used method of integrating the equations of motions in molecular
dynamics is Verlet‘s method or algorithm. It uses positions and velocities from
previously evaluated steps to provide a direct solution of equation (2.4). In this study,
the simulation code LAMMPS is used and its default algorithm is the velocity Verlet.
[35] Hence, the basic idea of Verlet‘s algorithm, and specifically the velocity Verlet,
will be presented here. Information on other integration schemes can be found in texts
on atomic scale simulation. [32]
This standard, or basic, Verlet algorithm is a combination of two Taylor‘s
expansions as following. The Taylor series for position from time t forward to 𝑡 + ∆𝑡

 And the series form t backward to 𝑡 − ∆𝑡 is

 Adding equations (2.6) and (2.7) with some simplifying produces

 (2.8)

This provides Verlet‘s algorithm for positions. The truncation error of the algorithm
when evolving the system by ∆t is of the order of ∆t^4. The acceleration needed in equation (2.8) is obtained from the intermolecular forces and Newton‘s second law. Velocities are not explicitly solved in equation (2.8), and they are calculated from first-order central difference.

A more commonly used, related algorithm is the velocity Verlet algorithm. The
approach is similar with the standard one but starts with position and velocity
expansions explicitly

The implement scheme of this algorithm is following:
1. Calculate velocities at mid-step,
Calculate the positions at next step 𝑡 + ∆t
2. Calculate the accelerations at 𝑡 + ∆𝑡, 𝑎(𝑡 + ∆𝑡) from the interaction potential using 𝑥(𝑡 + ∆𝑡).
3. Calculate the velocities at 𝑡 + ∆t

The advantage of this algorithm compared to the standard one is that we never need to simultaneously store the values at two different times for any of these
quantities. The velocity Verlet algorithm gives good feature of time reversibility,
simplicity and stability for even larger time steps compared to standard algorithm.



Ref: http://preserve.lehigh.edu/cgi/viewcontent.cgi?article=2257&context=etd
Molecular Dynamics Simulation of Gold Nanoparticles and Surface Stress Effect
Siming Zhang
Lehigh University
2011


---------------------------------------------------
Verlet is an efficient algorithm, and its benefits are stability, time-reversibility and
energy conservation.
The time step, ∆t, should be chosen in such manner that the simulation system
is stable and energy is conserved, that is to say that the provided ∆t should be
short enough. Usually, in classical MD simulations time steps of order 1–2 fs are
being used (especially when using constraints), but in some cases time steps can
be even 4–5 fs (when using virtual site algorithms). The length of the time
step is dominated by the fastest degrees of motion in the system, and usually theses
are bond vibrations. Most often constraints are used in simulations to reduce the
degrees of freedom, removing high-frequency vibrations by restricting the movement
of covalent bonds. This is implemented by fixing the covalent bonds to constant
angles and lengths. In addition to constraints, also restraints can be used to fix

atomic positions, bond lengths, bond angles and dihedrals.

introduction to Molecular dynamics, part 1

 Computer simulation methods have been applied to solve many-body problems in
physics, physical chemistry and biophysics; these methods are essentially experiments
in statistical mechanics. [30] The systems which statistical mechanics addresses are
composed of a large (i.e. statistical) number of particles; this branch of science allows
one to relate the collective properties of interesting atoms and molecules to the
macroscopic bulk (i.e. thermodynamic) properties of materials. In many cases, it is
impossible to study specific aspects of how atomic interactions manifest continuum
properties unless via computer simulation. Among simulation techniques, it is natural
to classify them according to the length scale at which the model is resolved. For
atomic scale models, there are two main families of simulation technique, Molecular
Dynamics (MD) method and Monte Carlo (MC) method [31, 32].
Molecular Dynamics simulation is a computational tool used to describe how
positions, velocities, and orientations of interacting particles change over time.
Particles may represent atoms, groups of atoms, or even larger entities; herein, focus
is on the use of MD as an atomic scale simulation technique. A set of models (i.e.
mathematic functions) that describe atomic scale interaction are foundations of any
MD simulation. Dynamics are governed by the system‘s Hamiltonian, and
Hamilton‘s equations of motion:

where H is the Hamiltonian, T is the kinetic energy which is a function of momenta, T=p*p/(2m), and U is the potential energy which is a function of position in the atomic
ensemble, U=U(q). The simulation model for atomic interactions is defined by an
interatomic potential energy function, or potential for brevity. [33] The potential can
be represented as 𝑈 = 𝑈(𝑟⃗ 𝑁), where 𝑟⃗ 𝑁 is the set of location vectors of an ensemble
of N atoms. For simplicity, in many simulations only pair interactions among atoms
are taken into account, so the potential energy of a system of N atoms is:

where 𝑢(𝑟𝑖𝑗) is a pair potential function having a known form, and 𝑟𝑖𝑗 is the
absolute value of distance between atom i and j. Since the interatomic forces are
derived from potential, they are all conservative, and can be related to the potential as
the gradient of U, or

In Molecular Dynamics, positions of atoms are attained by seeking the solution of
Newton‘s equation of motion:

where 𝐹𝑖
⃗⃗ is the force induced by the interactions between atom i and the other N-1
atoms, and the double dots mean second order time derivative while m is the mass of
each atom. Substitution of (2.3) into (2.4) relates the force to the potential function,
which is a known input to a given MD simulation. Integration once yields the
momenta of atoms, and the atomic position can be found by integrating an additional
time.
Repeating the integration procedure many times will reveal the trajectories of
atoms, based on which time average of some desired thermodynamic value, , can
be computed

When a state of equilibrium is reached, the average will be independent from the
initial time 𝑡0. Because there are many times of integration involved to achieve
equilibrium, it is of great importance that an appropriate integration scheme and a
time steps to the algorithm are selected. Poor integration schemes or too large a time
step will introduce systematic errors in the thermodynamic ensemble and therefore in
the computed value of .
MD simulation can be used to realize different sorts of thermodynamic ensembles,
which dictate the state thermodynamic quantities that will remain constant throughout
the simulation. Once the ensemble is established, time averages may be obtained, say,
for the thermodynamic state variables that are not held constant. For example, a
commonly employed MD ensemble is the microcanonical, [30] where number of
atoms N, volume V, and total energy E (sum of kinetic and potential energy) are held
constant (i.e. NVE ensemble); in such an ensemble, one may seek, for instance, time
averages of the system pressure P or temperature T. Other often employed ensembles
include the canonical, i.e. constant-NVT, ensemble, and the isothermal/isobaric, i.e. constant-NPT, ensemble. In order to implement a given ensemble in a MD simulation,
certain thermodynamic quantities must be held constant (either absolutely or in a time
averaged sense). Various numerical algorithms have been formulated to fix the
corresponding properties to a desired value. More theoretically detailed information
on ensembles can be found in [32] [33] and texts of statistical mechanics. The key
point is that properly implemental MD algorithms maintain thermodynamic
ensembles, form which methods of statistical mechanics reveal desired
thermodynamic quantities (i.e. material properties).
In short, the mechanism of the MD simulation process can be summarized as
following
1. Select the potential energy function that will be used to describe the
interatomic interactions
2. Set the initial position and velocity of each atom in a system and the
dimensions of the simulation space
3. Calculate all the forces on total N atoms via the governing interaction model
4. Find accelerations of each atom
5. Obtain each atom‘s velocity and position after a given timestep through
numerical integration by a chosen algorithm
6. Update all atomic positions
7. Repeat steps 3 - 6
8. Post process, or analyze, generated thermodynamic and atomic trajectory data


MD simulation was first introduced by B. Alder and T. Wainright in 1957 to study
the phase transition of hard spheres system, where particle propagation is determined
by each successive binary collision, rather than an integration of the equations of
motion. In 1964 a model system of liquid Argon was simulated by A. Rahman. In this
simulation not only the collision events between two particles but also interactions
modeled by Lennard-Jones potential were taken into account. This was the first time
that MD simulation was applied to atoms interacting through a continuous potential.
During 1970s, MD method continued to mature. Rahman and Stillinger (1971)
performed a simulation of 216 water molecules using an effective pair potential and
Ewald summation to calculate long-ranged forces triggered by electrostatic (i.e.
Coulomb) interactions. Several years later, the first protein simulation was carried out
by McCammon for bovine pancreatic tyrpsin inhibitor. [33] After that, MD
simulations began to be operated under different conditions and ensembles and
became more sophisticated. For instance, using MD to study non equilibrium
processes was first advanced in 1983 by Gillan and Dixon to calculate the
autocorrelation function of the microscopic heat current and hence the thermal
conductivity of a liquid. The Nose-Hoover thermostat algorithm, which is used to
obtain isothermal ensembles, was developed in 1984. The Car-Parrinello method to
carry out ab initio quantum mechanical MD (in which quantum mechanical effects for
the electronic degrees of freedom are taken into account) was developed in 1985. The
progress of available computer facilities is a major factor driving the development of
MD simulations. [31,32] Today MD methods are applied to a huge class of subjects,
e.g. surface properties, molecular clusters, biomolecules, and properties of liquids and
solids, and many different simulation codes have been developed.
Another well-developed simulation methodology is Monte Carlo (MC) method.
MC simulation is a computer simulation method that explores the phase space of a
system in a time-independent and stochastic manner. It uses a sequence of random
numbers to render and sample the phase space, from which diverse physical
properties of the system can be extracted. The main advantage of Monte Carlo
simulation over many analytical techniques is that it is a simple and straightforward
method that can provide an approximate solution to the problems encountered in
complex systems. Specifically, MC is able to reveal long time equilibrium properties.
The main advantage of MD method, with respect to MC method, is that MD
probes equilibrium configurations as well as dynamical properties of the system; this
makes MD well suited to compute time correlation functions and associated transport
coefficients. MD is therefore more efficient to evaluate some thermodynamic
properties like heat capacity, compressibility, and interfacial transport properties.
There are also computational advantages of MD method, particularly as compared to
ab initio methods. For example, the dynamic behavior of the atomic system is
described empirically without having to solve Schrodinger‘s equation at each time
step. Furthermore, it has been extensively discussed that many MD models are
inherently able to be executed in parallel. On the other hand, parallel implementations
of MC or ab initio based MD are highly non-trivial to realize.
Computing every force on each atom is a demanding part of a MD simulation. In
order to limit the number of force calculations, the smallest possible system – in terms
of number of atoms – should be simulated. The minimum system size necessary to
compute a given desired property is not always obvious and a common practice is to
analyze how increasing system size influences obtained results. A number of
algorithmic approaches exist that help minimize system size such as applying periodic
boundary conditions; these will be discussed in greater detail in the following section.
In addition to how large a system must be, one must determine how long a simulation
must be run to obtain reliable answers (i.e. with sufficient accuracy and precision).
Related to this is that the MD simulation time step should be small compared to the
period associated with the highest frequency atomic event. A typical time step for MD
simulations of metals is of order 1 fs (i.e. 1 E-15 s); this is roughly 1/1000 of the
period for an atomic oscillation. The very small time step is the main drawback of
MD method compared to other commonly used numerical simulations. Processes with
time scale of order microseconds and larger are outside the reasonable realm for even
highly parallel, efficient MD codes. Conversely, MC techniques cannot probe
dynamics but they can probe long time equilibrium properties.

Figure 2.1 Schematic comparison of time- and length-scales, accessible to different

types of simulation techniques (quantum simulations (QM), molecular dynamics

(MD), Brownian Dynamics (BD) and hydrodynamics/fluid dynamics (HD))

Accessible length- and time- scales of microscopic simulation should be
considered. Figure 2.1 shows a schematic representation for different types of
simulations in a length-time-diagram. It is clear that the more detailed a simulation
technique addresses, the smaller is the accessibility of times and length scales. Hence,
quantum simulations (QS), where fast motion of electrons are considered, are located
in the lower left corner of the diagram and typical length time scales are of order of Å
(Å) and ps. MD method puts fixed partial charges on interaction sites or adds a
dynamics model for polarization effects to simulate electronic distribution. [31]
Hence, the time scale of the system is determined by the time of collision between
atoms and atomic vibration frequencies. Consequently, the accessible time scale is of
order ns and corresponding accessible length scale becomes of order 10-10000 Å.
Note Figure 2.1 was originally presented in 2002; accessible length scales for MD
have increased dramatically. Brownian Dynamics (BD) can be applied to deal with
particles in a solvent medium, where one is not interested in a detailed description of
the solvent. In BD the effect of the solvent is hidden in average quantities. If one is
not interested in a molecular picture but in macroscopic quantities, the concepts of
hydrodynamics (HD) may be applied, where the system properties are hidden in
effective numbers, like density and sound velocity. These example methods, or levels
of theory, are geared towards liquids but analogous examples exist for solids.






Ref: http://preserve.lehigh.edu/cgi/viewcontent.cgi?article=2257&context=etd
Molecular Dynamics Simulation of Gold Nanoparticles and Surface Stress Effect
Siming Zhang
Lehigh University
2011

gold nanoparticles. Surface stress

Among nanoparticles, gold nanoparticles (Au NPs) have been extensively studied;
this is partly because gold is the subject of one of the most ancient themes of
investigation in science. Au NPs are among the most stable metal nanoparticles [5]
since Au is relatively chemically inert. They present fascinating aspects such as
low-symmetry structures at some geometric magic numbers [9] and size-related
electronic, magnetic [5, 7] and optical properties (due to quantum size effect).
Moreover, their applications to catalysis, owing to the high surface-to-volume ratio,
and biology are also significant. For example, conjugates of Au NPs-oligonucleotides
are of great interest owing to the potential use of the programmability of DNA
base-pairing to organize nanocrystals in space. This gives multiple ways of providing
a signature for the detection of precise DNA sequences to develop biosensors, disease
diagnosis, and gene expression, etc [5]. Besides, research is expanding on Au NP
catalytic effects associated with CO oxidation, NO reduction, and the water-gas shift
reaction, i.e. the chemical reaction in which carbon monoxide reacts with water vapor
to form carbon dioxide and hydrogen; new Au NP based catalytic systems are now
being explored. [4][5]
Applications exploiting the optical properties of Au NPs utilize functionalization
of NPs with chromophores. Chromophore means the part of a molecule responsible for its color. These have diverse applications due to a range of options for the
chromophore. In 2003, K. Thomas et al. reported that gold nanoparticles associated
with fluorophores were utilized in photocurrent generation and fluorescent display
devices [10]. Furthermore, they also showed that gold nanoparticles can bind and
release amino acids when linked with appropriate chromophores. All in all, it appears
safe to assume that Au NPs will be a key building block for nano-science and
-technology in the 21st century.
Given the many applications of Au NPs, as well as the extensive literature on
them, this thesis will focus on Au NPs as a canonical metal NP system.


The structure of gold nanoparticles has been perhaps the most investigated aspect
about them.[4] It has been found, theoretically and experimentally that gold
nanoparticles have crystallographic structures different from the bulk material.[4,11]
In general, gold nanoparticles have icosahedral or decahedral motifs depending on
their size while bulk gold has face-centered cubic (fcc) crystal structure [4].
For Au NPs larger than 2 nm, the truncated-octahedral structure (fcc) is the dominant shape [4, 12].
-----------------
Because NPs in general - and Au NPs specifically - possess significant surface
area to volume ratio, their surface thermodynamic properties can be very important in
determining a NP‘s properties. An important surface thermodynamic quantity is the
surface stress. Surface stress, f, is a thermodynamic quantity that describes the amount
of energy or reversible work per unit area required to elastically deform a solid
surface. It differs from another fundamental thermodynamic parameter, i.e. surface
free energy γ, which represents the energy needed to form a new surface by a process
like cleavage.[24] In liquid, these two values are identical, as the configuration of
fluid surface remains constant owing to bulk atoms or molecules moving exteriorly to
the surface when a fluid surface is stretched. Put differently, a liquid cannot support a
shear stress so the only way to create/eliminate new surface is by adding/removing
atoms from the surface, rather than elastically deforming existing atoms at the surface.
In contrast, when a solid surface is put in tension (within an elastic limit), the total
number of surface atoms is conserved; therefore, the number of atoms per unit area
changes and consequently, f ≠ γ.


Ref: http://preserve.lehigh.edu/cgi/viewcontent.cgi?article=2257&context=etd
Molecular Dynamics Simulation of Gold Nanoparticles and Surface Stress Effect
Siming Zhang
Lehigh University
2011

Au-Au Metallic Bonding Potential: The second-moment approximation to the tight-binding potential (TB-SMA)

The second-moment approximation to the tight-binding potential (TB-SMA) 89 is applied throughout all chapters of this dissertation to describe Au-Au interactions. Simple pairwise potentials such as the 12-6 Lennard-Jones potential fail to properly describe many of the properties (e.g., vacancy formation energies, surface structure, and relaxation properties) of transition metals. 89 Semi-empirical potentials, whose functional forms are derived from electronic structure considerations and then fit to experimental data, are better suited for simulations of transition metals. For instance, TB-SMA contains a many-body term that is modeled after the square-root dependence of the band energy on the second moment electron density of state:

where EiB is the many-body energy of atom i.
TB-SMA also contains a pairwise repulsive term given by

The total TB-SMA energy is then

Values for the parameters A, ξ , p, q, and r0 for Au are shown in Table 3.1, and are obtained from fits to the Au experimental cohesive energy, lattice parameter, and elastic constant [Cleri, F.; Rosato, V. Phys. Rev. B 1993, 48, 22–33].
An energy cutoff, rcut, of 5.8 A, is applied such that any pair of Au atoms separated by a distance greater than rcut do not interact. Differentiating equation 3.1 yields an expression for force that depends on the electron density, ρ, of atoms i and j, where

Thus, the total force acting on each atom is calculated in two stages, with each stage looping over atom i’s neighbors. This amounts to an additional computational cost compared to pairwise interaction models, where only a single loop over atom i’s neighbors is performed.

Ref: http://etd.library.vanderbilt.edu/available/etd-03212013-122215/unrestricted/French.pdf

Why does Intel CPU display 8 cores instead of true 4 cores? Hyper-threading is the answer.

Hyper-Threading Technology HTT is Intel's proprietary simultaneous multithreading (SMT) implementation used to improve parallelization of computations performed on x86 microprocessors.

For each processor core that is physically present, the operating system addresses two virtual or logical cores, and shares the workload between them when possible. The main function of hyper-threading is to increase the number of independent instructions in the pipeline; it takes advantage of superscalar architecture, in which multiple instructions operate on separate data in parallel. With HTT, one physical core appears as 2 processors to the operating system, which can use each core to schedule 2 processes at once. In addition, 2 or more processes can use the same resources: if resources for one process are not available, then another process can continue if its resources are available.
Architecturally, a processor with HTT consists of 2 logical processors per core, each of which has its own processor architectural state. Each logical processor can be individually halted, interrupted or directed to execute a specified thread, independently from the other logical processor sharing the same physical core.

Hyper-threading is a way of taking that 2nd process of task sharing away from the operating system, at least for a little bit. The reason is that the processor can task switch much faster than the operating system can tell it to. So by presenting two logical processors when there is in fact only one, the operating system has to do only half as many task switching operations, but more CPU scheduling. But the net result is supposedly a faster machine in multi-tasking operations.

The actual benefits of this vary greatly depending on the workload that you're doing. For most people, there is no harm in leaving it on or turning it off. 4 vs 8 threads is like choosing between a ferrari or a lamborghini for home users.

However, in server scenarios, it can make a large difference. For example, hypervisors can often get a large increase in speed through the use of hyperthreading, as they have very strict CPU scheduling requirements.

pair potential style and LAMMPS (pairwise interactions and many-body effects)

In a script we have:
#potential

pair_style eam
pair_coeff * * Au_u3.eam

The command "pair_style" is very important.


Set the formula(s) LAMMPS uses to compute pairwise interactions. In LAMMPS, pair potentials are defined between pairs of atoms that are within a cutoff distance and the set of active interactions typically changes over time. See the bond_style command to define potentials between pairs of bonded atoms, which typically remain in place for the duration of a simulation.

In LAMMPS, pairwise force fields encompass a variety of interactions, some of which include many-body effects, e.g. EAM, Stillinger-Weber, Tersoff, REBO potentials. They are still classified as "pairwise" potentials because the set of interacting atoms changes with time (unlike molecular bonds) and thus a neighbor list is used to find nearby interacting atoms.

In the formulas listed for each pair style, E is the energy of a pairwise interaction between two atoms separated by a distance r. The force between the atoms is the negative derivative of this expression:

-∂E/∂r is a force

If the pair_style command has a cutoff argument, it sets global cutoffs for all pairs of atom types.

There are also additional accelerated pair styles (not listed on the pair_style doc page) included in the LAMMPS distribution for faster performance on CPUs and GPUs. The list of these with links to the individual styles are given in the pair section of this page.

The list of pair_style (LAMMPS 2015 january):

• pair_style none - turn off pairwise interactions
• pair_style hybrid - multiple styles of pairwise interactions
• pair_style hybrid/overlay - multiple styles of superposed pairwise interactions
• pair_style adp - angular dependent potential (ADP) of Mishin
• pair_style airebo - AIREBO potential of Stuart
• pair_style beck - Beck potential
• pair_style body - interactions between body particles
• pair_style bop - BOP potential of Pettifor
• pair_style born - Born-Mayer-Huggins potential
• pair_style born/coul/long - Born-Mayer-Huggins with long-range Coulombics
• pair_style born/coul/msm - Born-Mayer-Huggins with long-range MSM Coulombics
• pair_style born/coul/wolf - Born-Mayer-Huggins with Coulombics via Wolf potential
• pair_style brownian - Brownian potential for Fast Lubrication Dynamics
• pair_style brownian/poly - Brownian potential for Fast Lubrication Dynamics with polydispersity
• pair_style buck - Buckingham potential
• pair_style buck/coul/cut - Buckingham with cutoff Coulomb
• pair_style buck/coul/long - Buckingham with long-range Coulombics
• pair_style buck/coul/msm - Buckingham long-range MSM Coulombics
• pair_style buck/long/coul/long - long-range Buckingham with long-range Coulombics
• pair_style colloid - integrated colloidal potential
• pair_style comb - charge-optimized many-body (COMB) potential
• pair_style comb3 - charge-optimized many-body (COMB3) potential
• pair_style coul/cut - cutoff Coulombic potential
• pair_style coul/debye - cutoff Coulombic potential with Debye screening
• pair_style coul/dsf - Coulombics via damped shifted forces
• pair_style coul/long - long-range Coulombic potential
• pair_style coul/msm - long-range MSM Coulombics
• pair_style coul/wolf - Coulombics via Wolf potential
• pair_style dpd - dissipative particle dynamics (DPD)
• pair_style dpd/tstat - DPD thermostatting
• pair_style dsmc - Direct Simulation Monte Carlo (DSMC)
• pair_style eam - embedded atom method (EAM)
• pair_style eam/alloy - alloy EAM
• pair_style eam/fs - Finnis-Sinclair EAM
• pair_style eim - embedded ion method (EIM)
• pair_style gauss - Gaussian potential
• pair_style gayberne - Gay-Berne ellipsoidal potential
• pair_style gran/hertz/history - granular potential with Hertzian interactions
• pair_style gran/hooke - granular potential with history effects
• pair_style gran/hooke/history - granular potential without history effects
• pair_style hbond/dreiding/lj - DREIDING hydrogen bonding LJ potential
• pair_style hbond/dreiding/morse - DREIDING hydrogen bonding Morse potential
• pair_style kim - interface to potentials provided by KIM project
• pair_style lcbop - long-range bond-order potential (LCBOP)
• pair_style line/lj - LJ potential between line segments
• pair_style lj/charmm/coul/charmm - CHARMM potential with cutoff Coulomb
• pair_style lj/charmm/coul/charmm/implicit - CHARMM for implicit solvent
• pair_style lj/charmm/coul/long - CHARMM with long-range Coulomb
• pair_style lj/charmm/coul/msm - CHARMM with long-range MSM Coulombics
• pair_style lj/class2 - COMPASS (class 2) force field with no Coulomb
• pair_style lj/class2/coul/cut - COMPASS with cutoff Coulomb
• pair_style lj/class2/coul/long - COMPASS with long-range Coulomb
• pair_style lj/cut - cutoff Lennard-Jones potential with no Coulomb
• pair_style lj/cut/coul/cut - LJ with cutoff Coulomb
• pair_style lj/cut/coul/debye - LJ with Debye screening added to Coulomb
• pair_style lj/cut/coul/dsf - LJ with Coulombics via damped shifted forces
• pair_style lj/cut/coul/long - LJ with long-range Coulombics
• pair_style lj/cut/coul/msm - LJ with long-range MSM Coulombics
• pair_style dipole/cut - point dipoles with cutoff
• pair_style dipole/long - point dipoles with long-range Ewald
• pair_style lj/cut/tip4p/cut - LJ with cutoff Coulomb for TIP4P water
• pair_style lj/cut/tip4p/long - LJ with long-range Coulomb for TIP4P water
• pair_style lj/expand - Lennard-Jones for variable size particles
• pair_style lj/gromacs - GROMACS-style Lennard-Jones potential
• pair_style lj/gromacs/coul/gromacs - GROMACS-style LJ and Coulombic potential
• pair_style lj/long/coul/long - long-range LJ and long-range Coulombics
• pair_style lj/long/dipole/long - long-range LJ and long-range point dipoles
• pair_style lj/long/tip4p/long - long-range LJ and long-range Coulomb for TIP4P water
• pair_style lj/smooth - smoothed Lennard-Jones potential
• pair_style lj/smooth/linear - linear smoothed Lennard-Jones potential
• pair_style lj96/cut - Lennard-Jones 9/6 potential
• pair_style lubricate - hydrodynamic lubrication forces
• pair_style lubricate/poly - hydrodynamic lubrication forces with polydispersity
• pair_style lubricateU - hydrodynamic lubrication forces for Fast Lubrication Dynamics
• pair_style lubricateU/poly - hydrodynamic lubrication forces for Fast Lubrication
                   with polydispersity
• pair_style meam - modified embedded atom method (MEAM)
• pair_style mie/cut - Mie potential
• pair_style morse - Morse potential
• pair_style nb3b/harmonic - nonbonded 3-body harmonic potential
• pair_style nm/cut - N-M potential
• pair_style nm/cut/coul/cut - N-M potential with cutoff Coulomb
• pair_style nm/cut/coul/long - N-M potential with long-range Coulombics
• pair_style peri/eps - peridynamic EPS potential
• pair_style peri/lps - peridynamic LPS potential
• pair_style peri/pmb - peridynamic PMB potential
• pair_style peri/ves - peridynamic VES potential
• pair_style reax - ReaxFF potential
• pair_style rebo - 2nd generation REBO potential of Brenner
• pair_style resquared - Everaers RE-Squared ellipsoidal potential
• pair_style snap - SNAP quantum-accurate potential • pair_style soft - Soft (cosine) potential
• pair_style sw - Stillinger-Weber 3-body potential
• pair_style table - tabulated pair potential
• pair_style tersoff - Tersoff 3-body potential
• pair_style tersoff/mod - modified Tersoff 3-body potential
• pair_style tersoff/zbl - Tersoff/ZBL 3-body potential
• pair_style tip4p/cut - Coulomb for TIP4P water w/out LJ
• pair_style tip4p/long - long-range Coulombics for TIP4P water w/out LJ
• pair_style tri/lj - LJ potential between triangles
• pair_style yukawa - Yukawa potential
• pair_style yukawa/colloid - screened Yukawa potential for finite-size particles
• pair_style zbl - Ziegler-Biersack-Littmark potential

potential files, directories and LAMMPS

In a script we have:
#potential

pair_style eam
pair_coeff * * Au_u3.eam

The command "pair_coeff" 
When a pair_coeff command using a potential file is specified, LAMMPS looks for the potential file in 2 places.

• First it looks in the location specified. E.g. if the file is specified as "Au_u3.eam", it is looked for in the current working directory. 
• If it is specified as "../potentials/niu3.eam", then it is looked for in the potentials directory, assuming it is a sister directory of the current working directory. 
• If the file is not found, it is then looked for in the directory specified by the LAMMPS_POTENTIALS environment variable. 
• Thus if this is set to the potentials directory in the LAMMPS distro, then you can use those files from anywhere on your system, without copying them into your working directory. Environment variables are set in different ways for different shells. Here are example settings for

csh, tcsh:
% setenv LAMMPS_POTENTIALS /path/to/lammps/potentials
bash:
% export LAMMPS_POTENTIALS=/path/to/lammps/potentials
Windows:
% set LAMMPS_POTENTIALS="C:\Path to LAMMPS\Potentials

-------------------------------
Styles with a cuda, gpu, intel, kk, omp, or opt suffix are functionally the same as the corresponding style without the suffix. They have been optimized to run faster, depending on your available hardware. The accelerated styles take the same arguments and should produce the same results, except for round-off and precision issues.
These accelerated styles are part of the USER-CUDA, GPU, USER-INTEL, KOKKOS, USER-OMP and OPT packages, respectively. They are only enabled if LAMMPS was built with those packages.
You can specify the accelerated styles explicitly in your input script by including their suffix, or you can use the -suffix command-line switch when you invoke LAMMPS, or you can use the suffix command in your input script.
-------------------------------
As an example, imagine a file SiC.sw has Stillinger-Weber values for Si and C. If your LAMMPS simulation has 4 atoms types and you want the 1st 3 to be Si, and the 4th to be C, you would use the following pair_coeff command:

pair_coeff * * SiC.sw Si Si Si C

The 1st 2 arguments must be * * so as to span all LAMMPS atom types. 
The first three Si arguments map LAMMPS atom types 1,2,3 to the Si element in the SW file. The final C argument maps LAMMPS atom type 4 to the C element in the SW file. If a mapping value is specified as NULL, the mapping is not performed. This can be used when a sw potential is used as part of the hybrid pair style. The NULL values are placeholders for atom types that will be used with other potentials.
Stillinger-Weber files in the potentials directory of the LAMMPS distribution have a ".sw" suffix. Lines that are not blank or comments (starting with #) define parameters for a triplet of elements. The parameters in a single entry correspond to the two-body and three-body coefficients in the formula above:
• element 1 (the center atom in a 3-body interaction)
• element 2
• element 3
• epsilon (energy units)
• sigma (distance units)
•a
• lambda
• gamma
• costheta0
•A
•B
•p
•q
• tol
The A, B, p, and q parameters are used only for two-body interactions. The lambda and costheta0 parameters are used only for three-body interactions. The epsilon, sigma and a parameters are used for both two-body and three-body interactions. gamma is used only in the three-body interactions, but is defined for pairs of atoms. The non-annotated parameters are unitless.
LAMMPS introduces an additional performance-optimization parameter tol that is used for both two-body and three-body interactions. In the Stillinger-Weber potential, the interaction energies become negligibly small at atomic separations substantially less than the theoretical cutoff distances. LAMMPS therefore defines a virtual cutoff distance based on a user defined tolerance tol. The use of the virtual cutoff distance in constructing atom neighbor lists can significantly reduce the neighbor list sizes and therefore the computational cost. LAMMPS provides a tol value for each of the three-body entries so that they can be separately controlled. If tol = 0.0, then the standard Stillinger-Weber cutoff is used.
The Stillinger-Weber potential file must contain entries for all the elements listed in the pair_coeff command. It can also contain entries for additional elements not being used in a particular simulation; LAMMPS ignores those entries.
For a single-element simulation, only a single entry is required (e.g. SiSiSi). For a two-element simulation, the file must contain 8 entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that specify SW parameters for all permutations of the two elements interacting in three-body configurations. Thus for 3 elements, 27 entries would be required, etc.
As annotated above, the first element in the entry is the center atom in a three-body interaction. Thus an entry for SiCC means a Si atom with 2 C atoms as neighbors. The parameter values used for the two-body interaction come from the entry where the 2nd and 3rd elements are the same. Thus the two-body parameters for Si interacting with C, comes from the SiCC entry. The three-body parameters can in principle be specific to the three elements of the configuration. In the literature, however, the three-body parameters are usually defined by simple formulas involving two sets of pair-wise parameters, corresponding to the ij and ik pairs, where i is the center atom. The user must ensure that the correct combining rule is used to calculate the values of the threebody parameters for alloys. Note also that the function phi3 contains two exponential screening factors with parameter values from the ij pair and ik pairs. So phi3 for a C atom bonded to a Si atom and a second C atom will depend on the three-body parameters for the CSiC entry, and also on the two-body parameters for the CCC and CSiSi entries. Since the order of the two neighbors is arbitrary, the threebody parameters for entries CSiC and CCSi should be the same. Similarly, the two-body parameters for entries SiCC and CSiSi should also be the same. The parameters used only for two-body interactions (A, B, p, and q) in entries whose 2nd and 3rd element are different (e.g. SiCSi) are not used for anything and can be set to 0.0 if desired. This is also true for the parameters in phi3 that are taken from the ij and ik pairs (sigma, a, gamma)