Try a view of this blog with HTML5 AJAX CSS3 (firefox3.6+, Chrome, Safari): http://ex-ample.blogspot.com/view/flipcardhttp://ex-ample.blogspot.com/view/flipcard
The left-sidebar and the bottom of this blog show small javascripts (js) examples; bottom: interactive services.

## Wednesday, May 20, 2015

### 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

## Monday, May 18, 2015

### 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

## Wednesday, May 13, 2015

### 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 mol1 and σ = 2.934Å
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/2Å2 where d denotes the equilibrium spacing which can be obtained as d = 21/6σ = 3.2933Å. Consequently, the mean surface density η of the gold layers is taken to be 1/A= 0.1065Å2.

### 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

## Thursday, May 7, 2015

### 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 "FT-H298." 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 TM = 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 T0 (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].

## Monday, April 27, 2015

### 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 Å (Å) 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 Å. 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.

### 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.

## Wednesday, April 22, 2015

### 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

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 Eiis 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