Sometimes I have to put text on a path

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

Molecular Dynamics Simulation of Gold Nanoparticles and Surface Stress Effect
Siming Zhang
Lehigh University

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 ≠ Ξ³.

Molecular Dynamics Simulation of Gold Nanoparticles and Surface Stress Effect
Siming Zhang
Lehigh University

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.


Tuesday, April 21, 2015

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.