Paradigms of Science

• Chinese Version
  

What is AI

Four Dimension of AI

• Systems that think like humans.
• Systems that think rationally.
• Systems that act like humans.
• Systems that act rationally.

How AI works?

• 7 Steps: Collect Data -> Prepare data -> Select model -> Train -> Evaluate -> Adjust Parameters -> Estimate

Strategic Goals and Objectives of MGI

Exam Suggestions: For answers to the problems on the exam sheet, just pick one below.

V1

• An Ecosystem of Collaborators as Opposed to Individual PI's
• Algorithms with SW Discipline
• Big data
• Advanced materials are the product
• An enhanced culture of materials design for all stakeholders - including the materials science and engineering community

Goal:To decrease the time-to-market by over 50%

V2

The success of the Materials Genome Initiative (MGI) will be achieved by meeting the following four goals:

1. Enable a Paradigm Shift in Culture
2. Integrate Experiments, Computation, and Theory
3. Facilitate Access to Materials Data
4. Equip the Next-Generation Materials Workforce

V3

Objectives:

1. Build foundational conditions and infrastructure for material innovation;
2. Enhance the enabling effect of data across the entire material value chain;
3. Comprehensively cultivate and develop a skilled workforce in the materials sector.

Computational Materials Science & The mission of Computational Materials Science

• Computational Materials Science, through atomistic simulation based on atomic-scale modeling and high-performance computing, serves as a bridge connecting theory and experiment, as well as the microscopic and macroscopic worlds.
  
  It has always been a dream of materials researchers to design a new material completely on paper, optimizing the composition and processing steps in order to achieve the properties required for a given application.

Typical Computational Simulation Methods and their underlying theory

• Quantum Mechanics -> First Principle
• Newton's 2nd Law -> Molecular Dynamics
• Statistic Mechanics -> Monte Carlo Simulation
• Meso-Landau Equation -> Phase Field
• Continuum Mechanics -> Finite Element Method

What can MD simulation be used

• Analyze Mechanical Behavior
• Calculate Dislocation Density
• Shock Compression
• Phase Transformation
• Irradiation Damage
• Predict/Understand exotic properties

Definition of Each Ensembles, Phase Space

Ensembles 系综

Definition: An ensemble is a fundamental concept in statistical mechanics used to describe the collection of all possible microscopic states of a macroscopic system under given macroscopic conditions. By studying ensembles, one can derive the macroscopic thermodynamic properties of a system (e.g., temperature, pressure, entropy).

1. Micro-canonical ensemble (NVE) 微正则系综

   • fixed number of atoms, N, a fixed volume, V, and a fixed energy, E.
   • This corresponds to an isolated system.

2. Canonical Ensemble (NVT) 正则系综

   • fixed number of atoms, N, a fixed volume, V, and a fixed temperature, T.

3. Isobaric-Isothermal Ensemble (NPT) 等压-等温系综

   • fixed number of atoms, N, a fixed pressure, P, and a fixed temperature, T.

4. Grand canonical Ensemble (mVT) 巨正则系综

   • fixed chemical potential, m, a fixed volume, V, and a fixed temperature, T.

相空间

Definition: Phase space is a space that represents all possible states of a system.

Ergodicity 各向遍历性

Definition: Ergodicity refers to the trajectory of a dynamical system that can pass through all possible states in phase space after a long enough time.

What is molecular dynamics simulations？What is its advantage and disadvantage？How to solve it?

• __Molecular dynamics (MD)__: Model the motion of some group of particles (e.g., atoms) by solving the classical equations of motion.

• Cons:

  • Validity -> Find more accurate potential
  • Limited Time Scale -> Find Suitable Cases; AMD(Accelerated Molecular Dynamics)
  • Limited Spatial Scale -> Multi-scale Calculation
• Pros have been mentioned above

Procedure of MD Simulation

Periodic Boundary Conditions

• The main reason is to remove the effects of thesurface.
• Whenever an atom leaves the simulation cell, it is replaced by another with exactly the same velocity, entering from the opposite cell face. So the number of atoms in the cell is conserved.
• Furthermore, no atom feels any surface forces, as these are now completely removed.

• $r_{\text{cut}}$ is the cutoff radius when calculating the force between two atoms.

  Coding

subroutine Refold_Positions
!
! Particles that left the box are refolded back into the box by
! periodic boundary conditions
!
use Particles
implicit none
where ( pos >  0.5d0 ) pos = pos - 1.d0
where ( pos < -0.5d0 ) pos = pos + 1.d0
end subroutine Refold_Positions

Crystal Building Through Coding

Place Ti atoms into B2 structure

/* Place Ti atoms into B2 structure */
n=0;
for(i=0;i<nsize;i++)
{
    for(j=0;j<nsize;j++)
    {
        for(k=0;k<nsize;k++)
        {
            RX[n]=(i+0.25)/nsize;
            RY[n]=(j+0.25)/nsize;
            RZ[n]=(k+0.25)/nsize;
            VX[n]=0.0; VY[n]=0.0; VZ[n]=0.0;
            AX[n]=0.0; AY[n]=0.0; AZ[n]=0.0;
            BX[n]=0.0; BY[n]=0.0; BZ[n]=0.0;
            CX[n]=0.0; CY[n]=0.0; CZ[n]=0.0;
            n=n+1;
        }
    }
}
NTi=n;
printf("Number of Ti atoms=%d \n",NTi);

Here is the combined and translated version of the two answers in English:

Explanation of the Code

This code is designed to place titanium (Ti) atoms into a B2 crystal structure and initialize their positions, velocities, accelerations, and other properties. Below is a detailed explanation of the code:

1. Variable Definitions

• n: Index used to track the current Ti atom being processed.
• i, j, k: Loop variables for traversing the three dimensions of the lattice.
• nsize: Size of the lattice in each dimension.
• RX, RY, RZ: Arrays storing the x, y, and z coordinates of the Ti atoms.
• VX, VY, VZ: Arrays storing the x, y, and z velocities of the Ti atoms.
• AX, AY, AZ: Arrays storing the x, y, and z accelerations of the Ti atoms.
• BX, BY, BZ: Arrays storing some property of the Ti atoms (specific meaning unclear).
• CX, CY, CZ: Arrays storing some property of the Ti atoms (specific meaning unclear).
• NTi: Variable storing the total number of Ti atoms.

2. Loop Structure

The code uses three nested loops to traverse each position in the lattice:

for(i=0;i<nsize;i++)
{
    for(j=0;j<nsize;j++)
    {
        for(k=0;k<nsize;k++)
        {
            // Initialize properties for each Ti atom
        }
    }
}

• i, j, k represent the indices of the lattice in the x, y, and z directions, respectively.
• Each loop runs from 0 to nsize-1, ensuring the entire lattice is traversed.

3. Initialization of Ti Atom Properties

At each lattice position, the code initializes the properties of a Ti atom:

RX[n]=(i+0.25)/nsize;
RY[n]=(j+0.25)/nsize;
RZ[n]=(k+0.25)/nsize;
VX[n]=0.0; VY[n]=0.0; VZ[n]=0.0;
AX[n]=0.0; AY[n]=0.0; AZ[n]=0.0;
BX[n]=0.0; BY[n]=0.0; BZ[n]=0.0;
CX[n]=0.0; CY[n]=0.0; CZ[n]=0.0;
n=n+1;

• RX[n], RY[n], RZ[n]: Set the position of the Ti atom to $(i+0.25)/nsize$, $(j+0.25)/nsize$, $(k+0.25)/nsize$. The $0.25$ offset places the Ti atom at a specific position in the B2 structure.
• VX[n], VY[n], VZ[n]: Initialize the velocity of the Ti atom to 0.
• AX[n], AY[n], AZ[n]: Initialize the acceleration of the Ti atom to 0.
• BX[n], BY[n], BZ[n]: Initialize some property of the Ti atom to 0.
• CX[n], CY[n], CZ[n]: Initialize some property of the Ti atom to 0.
• n=n+1: Increment the index of the Ti atom.

4. Output the Total Number of Ti Atoms

NTi=n;
printf("Number of Ti atoms=%d \n",NTi);

• NTi: Stores the total number of Ti atoms, which is the value of n.
• printf: Outputs the total number of Ti atoms.

5. Summary

This code uses three nested loops to traverse each position in the lattice and initializes the position, velocity, acceleration, and other properties of a Ti atom at each position. Finally, it outputs the total number of Ti atoms.

Coordinates Traversed When nsize = 2

When nsize = 2, the program traverses all positions in the lattice and initializes the coordinates of a Ti atom at each position. Since nsize = 2, the values of i, j, and k range from 0 to 1, resulting in a total of $2 \times 2 \times 2 = 8$ positions.

Coordinates Traversed

At each position, the coordinates of the Ti atom are calculated using the following formulas:

$$ RX[n] = \frac{i + 0.25}{nsize}, \quad RY[n] = \frac{j + 0.25}{nsize}, \quad RZ[n] = \frac{k + 0.25}{nsize} $$

Substituting nsize = 2, the formulas become:

$$ RX[n] = \frac{i + 0.25}{2}, \quad RY[n] = \frac{j + 0.25}{2}, \quad RZ[n] = \frac{k + 0.25}{2} $$

Specific Coordinates

The following table lists all the coordinates traversed by the program:

Summary

When nsize = 2, the program traverses 8 positions, and the coordinates of the Ti atoms at these positions are:

$$ (0.125, 0.125, 0.125), \\ (0.125, 0.125, 0.625), \\ (0.125, 0.625, 0.125), \\ (0.125, 0.625, 0.625), \\ (0.625, 0.125, 0.125), \\ (0.625, 0.125, 0.625), \\ (0.625, 0.625, 0.125), \\ (0.625, 0.625, 0.625) $$

These coordinates correspond to the positions of Ti atoms in the B2 structure.

This combined explanation provides a clear understanding of the code and the specific coordinates it generates when nsize = 2.

Place Ni atoms into B2 structure

/* Place Ni atoms into B2 structure */
for(i=0;i<nsize;i++)
{
    for(j=0;j<nsize;j++)
    {
        for(k=0;k<nsize;k++)
        {
            RX[n]=(i+0.75)/nsize;
            RY[n]=(j+0.75)/nsize;
            RZ[n]=(k+0.75)/nsize;
            VX[n]=0.0; VY[n]=0.0; VZ[n]=0.0;

            AX[n]=0.0; AY[n]=0.0; AZ[n]=0.0; BX[n]=0.0; BY[n]=0.0; BZ[n]=0.0;
            CX[n]=0.0; CY[n]=0.0; CZ[n]=0.0;
            n=n+1;
        }
    }
}
NNi=n-NTi;
NT=n;
NSITE=NT;
printf("Number of Ni atoms=%d \n",NNi);
printf("Total Number of Sites in B2=%d \n",NSITE);

Explanation of the Code

This code is designed to place nickel (Ni) atoms into a B2 crystal structure and initialize their positions, velocities, accelerations, and other properties. Below is a detailed explanation of the code:

1. Variable Definitions

• n: Index used to track the total number of atoms (including both Ti and Ni) being processed.
• i, j, k: Loop variables for traversing the three dimensions of the lattice.
• nsize: Size of the lattice in each dimension.
• RX, RY, RZ: Arrays storing the x, y, and z coordinates of the atoms.
• VX, VY, VZ: Arrays storing the x, y, and z velocities of the atoms.
• AX, AY, AZ: Arrays storing the x, y, and z accelerations of the atoms.
• BX, BY, BZ: Arrays storing some property of the atoms (specific meaning unclear).
• CX, CY, CZ: Arrays storing some property of the atoms (specific meaning unclear).
• NNi: Variable storing the total number of Ni atoms.
• NT: Variable storing the total number of atoms (including both Ti and Ni).
• NSITE: Variable storing the total number of sites in the lattice.

2. Loop Structure

The code uses three nested loops to traverse each position in the lattice:

for(i=0;i<nsize;i++)
{
    for(j=0;j<nsize;j++)
    {
        for(k=0;k<nsize;k++)
        {
            // Initialize properties for each Ni atom
        }
    }
}

• i, j, k represent the indices of the lattice in the x, y, and z directions, respectively.
• Each loop runs from 0 to nsize-1, ensuring the entire lattice is traversed.

3. Initialization of Ni Atom Properties

At each lattice position, the code initializes the properties of a Ni atom:

RX[n]=(i+0.75)/nsize;
RY[n]=(j+0.75)/nsize;
RZ[n]=(k+0.75)/nsize;
VX[n]=0.0; VY[n]=0.0; VZ[n]=0.0;
AX[n]=0.0; AY[n]=0.0; AZ[n]=0.0;
BX[n]=0.0; BY[n]=0.0; BZ[n]=0.0;
CX[n]=0.0; CY[n]=0.0; CZ[n]=0.0;
n=n+1;

• RX[n], RY[n], RZ[n]: Set the position of the Ni atom to $(i+0.75)/nsize$, $(j+0.75)/nsize$, $(k+0.75)/nsize$. The $0.75$ offset places the Ni atom at a specific position in the B2 structure.
• VX[n], VY[n], VZ[n]: Initialize the velocity of the Ni atom to 0.
• AX[n], AY[n], AZ[n]: Initialize the acceleration of the Ni atom to 0.
• BX[n], BY[n], BZ[n]: Initialize some property of the Ni atom to 0.
• CX[n], CY[n], CZ[n]: Initialize some property of the Ni atom to 0.
• n=n+1: Increment the total index of atoms.

4. Calculate the Total Number of Ni Atoms and Sites

NNi=n-NTi;
NT=n;
NSITE=NT;
printf("Number of Ni atoms=%d \n",NNi);
printf("Total Number of Sites in B2=%d \n",NSITE);

• NNi: Total number of Ni atoms, calculated by subtracting the total number of Ti atoms (NTi) from the current total index n.
• NT: Total number of atoms (including both Ti and Ni), which is the value of n.
• NSITE: Total number of sites in the lattice, equal to NT.
• printf: Outputs the total number of Ni atoms and the total number of sites in the lattice.

5. Coordinates Traversed When nsize = 2

When nsize = 2, the program traverses all positions in the lattice and initializes the coordinates of a Ni atom at each position. Since nsize = 2, the values of i, j, and k range from 0 to 1, resulting in a total of $2 \times 2 \times 2 = 8$ positions.

Traversed Coordinates

At each position, the coordinates of the Ni atom are calculated using the following formulas:

$$ RX[n] = \frac{i + 0.75}{nsize}, \quad RY[n] = \frac{j + 0.75}{nsize}, \quad RZ[n] = \frac{k + 0.75}{nsize} $$

Substituting nsize = 2, the formulas become:

$$ RX[n] = \frac{i + 0.75}{2}, \quad RY[n] = \frac{j + 0.75}{2}, \quad RZ[n] = \frac{k + 0.75}{2} $$

Specific Coordinates

The following table lists all the coordinates traversed by the program:

Summary

When nsize = 2, the program traverses 8 positions, and the coordinates of the Ni atoms at these positions are:

$$ (0.375, 0.375, 0.375), \\ (0.375, 0.375, 0.875), \\ (0.375, 0.875, 0.375), \\ (0.375, 0.875, 0.875), \\ (0.875, 0.375, 0.375), \\ (0.875, 0.375, 0.875), \\ (0.875, 0.875, 0.375), \\ (0.875, 0.875, 0.875) $$

These coordinates correspond to the positions of Ni atoms in the B2 structure.

6. Summary

This code uses three nested loops to traverse each position in the lattice and initializes the position, velocity, acceleration, and other properties of a Ni atom at each position. Finally, it outputs the total number of Ni atoms and the total number of sites in the lattice. When nsize = 2, the program traverses 8 positions and generates the corresponding coordinates for Ni atoms.

Single Crystal

Explanation of the Code for Creating a Single Crystal Structure

This code is designed to create a single crystal structure with one type of atom (e.g., a simple cubic structure) and initialize their positions, velocities, accelerations, and other properties. Below is a detailed explanation of the code:

1. Variable Definitions

• n: Index used to track the current atom being processed.
• i, j, k: Loop variables for traversing the three dimensions of the lattice.
• nsize: Size of the lattice in each dimension.
• RX, RY, RZ: Arrays storing the x, y, and z coordinates of the atoms.
• VX, VY, VZ: Arrays storing the x, y, and z velocities of the atoms.
• AX, AY, AZ: Arrays storing the x, y, and z accelerations of the atoms.
• BX, BY, BZ: Arrays storing some property of the atoms (specific meaning unclear).
• CX, CY, CZ: Arrays storing some property of the atoms (specific meaning unclear).
• NT: Variable storing the total number of atoms.
• NSITE: Variable storing the total number of sites in the lattice.

2. Loop Structure

The code uses three nested loops to traverse each position in the lattice:

for(i=0;i<nsize;i++)
{
    for(j=0;j<nsize;j++)
    {
        for(k=0;k<nsize;k++)
        {
            // Initialize properties for each atom
        }
    }
}

• i, j, k represent the indices of the lattice in the x, y, and z directions, respectively.
• Each loop runs from 0 to nsize-1, ensuring the entire lattice is traversed.

3. Initialization of Atom Properties

At each lattice position, the code initializes the properties of an atom:

RX[n]=i/nsize;
RY[n]=j/nsize;
RZ[n]=k/nsize;
VX[n]=0.0; VY[n]=0.0; VZ[n]=0.0;
AX[n]=0.0; AY[n]=0.0; AZ[n]=0.0;
BX[n]=0.0; BY[n]=0.0; BZ[n]=0.0;
CX[n]=0.0; CY[n]=0.0; CZ[n]=0.0;
n=n+1;

• RX[n], RY[n], RZ[n]: Set the position of the atom to $i/nsize$, $j/nsize$, $k/nsize$. This places the atom at the lattice points of a simple cubic structure.
• VX[n], VY[n], VZ[n]: Initialize the velocity of the atom to 0.
• AX[n], AY[n], AZ[n]: Initialize the acceleration of the atom to 0.
• BX[n], BY[n], BZ[n]: Initialize some property of the atom to 0.
• CX[n], CY[n], CZ[n]: Initialize some property of the atom to 0.
• n=n+1: Increment the index of the atom.

4. Calculate the Total Number of Atoms and Sites

NT=n;
NSITE=NT;
printf("Total Number of Atoms=%d \n",NT);
printf("Total Number of Sites in Single Crystal=%d \n",NSITE);

• NT: Total number of atoms, which is the value of n.
• NSITE: Total number of sites in the lattice, equal to NT.
• printf: Outputs the total number of atoms and the total number of sites in the lattice.

5. Coordinates Traversed When nsize = 2

When nsize = 2, the program traverses all positions in the lattice and initializes the coordinates of an atom at each position. Since nsize = 2, the values of i, j, and k range from 0 to 1, resulting in a total of $2 \times 2 \times 2 = 8$ positions.

Traversed Coordinates

At each position, the coordinates of the atom are calculated using the following formulas:

$$ RX[n] = \frac{i}{nsize}, \quad RY[n] = \frac{j}{nsize}, \quad RZ[n] = \frac{k}{nsize} $$

Substituting nsize = 2, the formulas become:

$$ RX[n] = \frac{i}{2}, \quad RY[n] = \frac{j}{2}, \quad RZ[n] = \frac{k}{2} $$

Specific Coordinates

The following table lists all the coordinates traversed by the program:

Summary

When nsize = 2, the program traverses 8 positions, and the coordinates of the atoms at these positions are:

$$ (0.0, 0.0, 0.0), \\ (0.0, 0.0, 0.5), \\ (0.0, 0.5, 0.0), \\ (0.0, 0.5, 0.5), \\ (0.5, 0.0, 0.0), \\ (0.5, 0.0, 0.5), \\ (0.5, 0.5, 0.0), \\ (0.5, 0.5, 0.5) $$

These coordinates correspond to the positions of atoms in a simple cubic structure.

6. Summary

This code uses three nested loops to traverse each position in the lattice and initializes the position, velocity, acceleration, and other properties of an atom at each position. Finally, it outputs the total number of atoms and the total number of sites in the lattice. When nsize = 2, the program traverses 8 positions and generates the corresponding coordinates for atoms in a simple cubic structure.

BCC Unit Cell

Here is the pseudo-code to generate atomic coordinates for a BCC (Body-Centered Cubic) unit cell with a single type of atom:

#include <stdio.h>

int main() {
    int nsize = 2; // Grid size
    int n = 0; // Atom counter
    double RX[1000], RY[1000], RZ[1000]; // Atom coordinate arrays
    int i, j, k;

    // Loop through the grid
    for (i = 0; i < nsize; i++) {
        for (j = 0; j < nsize; j++) {
            for (k = 0; k < nsize; k++) {
                // Corner atom
                RX[n] = (double)i / nsize;
                RY[n] = (double)j / nsize;
                RZ[n] = (double)k / nsize;
                n++;

                // Body-centered atom
                RX[n] = (i + 0.5) / nsize;
                RY[n] = (j + 0.5) / nsize;
                RZ[n] = (k + 0.5) / nsize;
                n++;
            }
        }
    }

    // Output atomic coordinates
    printf("Total number of atoms: %d\n", n);
    for (i = 0; i < n; i++) {
        printf("Atom %d: (%.3f, %.3f, %.3f)\n", i + 1, RX[i], RY[i], RZ[i]);
    }

    return 0;
}

Explanation of the Code

1. Grid Traversal:

   • Three nested loops iterate over an nsize x nsize x nsize grid.
   • For each grid point, a corner atom and a body-centered atom are generated.

2. Corner Atom:

   • The coordinates of the corner atom are:
     $$ RX = \frac{i}{nsize}, \quad RY = \frac{j}{nsize}, \quad RZ = \frac{k}{nsize} $$

3. Body-Centered Atom:

   • The coordinates of the body-centered atom are:
     $$ RX = \frac{i + 0.5}{nsize}, \quad RY = \frac{j + 0.5}{nsize}, \quad RZ = \frac{k + 0.5}{nsize} $$

4. Atom Counting:

   • The counter n increments by 1 for each atom generated.

5. Output:

   • The total number of atoms and their coordinates are printed.

Atomic Coordinates for nsize = 2

For nsize = 2, the grid size is 2 x 2 x 2, and the program generates the following atomic coordinates:

Summary

• The code generates a BCC unit cell with nsize = 2, resulting in 16 atoms.
• Atomic coordinates are calculated for both corner and body-centered positions.
• The output lists all atomic coordinates.

FCC

Here is the loop for placing atoms into an FCC structure, following the same format as your provided B2 structure code:

/* Place atoms into FCC structure */
for(i=0;i<nsize;i++)
{
    for(j=0;j<nsize;j++)
    {
        for(k=0;k<nsize;k++)
        {
            // Corner atom
            RX[n]=(double)i/nsize;
            RY[n]=(double)j/nsize;
            RZ[n]=(double)k/nsize;
            VX[n]=0.0; VY[n]=0.0; VZ[n]=0.0;
            AX[n]=0.0; AY[n]=0.0; AZ[n]=0.0; BX[n]=0.0; BY[n]=0.0; BZ[n]=0.0;
            CX[n]=0.0; CY[n]=0.0; CZ[n]=0.0;
            n=n+1;

            // Face-centered atoms
            RX[n]=(i+0.5)/nsize;
            RY[n]=(double)j/nsize;
            RZ[n]=(k+0.5)/nsize;
            VX[n]=0.0; VY[n]=0.0; VZ[n]=0.0;
            AX[n]=0.0; AY[n]=0.0; AZ[n]=0.0; BX[n]=0.0; BY[n]=0.0; BZ[n]=0.0;
            CX[n]=0.0; CY[n]=0.0; CZ[n]=0.0;
            n=n+1;

            RX[n]=(double)i/nsize;
            RY[n]=(j+0.5)/nsize;
            RZ[n]=(k+0.5)/nsize;
            VX[n]=0.0; VY[n]=0.0; VZ[n]=0.0;
            AX[n]=0.0; AY[n]=0.0; AZ[n]=0.0; BX[n]=0.0; BY[n]=0.0; BZ[n]=0.0;
            CX[n]=0.0; CY[n]=0.0; CZ[n]=0.0;
            n=n+1;

            RX[n]=(i+0.5)/nsize;
            RY[n]=(j+0.5)/nsize;
            RZ[n]=(double)k/nsize;
            VX[n]=0.0; VY[n]=0.0; VZ[n]=0.0;
            AX[n]=0.0; AY[n]=0.0; AZ[n]=0.0; BX[n]=0.0; BY[n]=0.0; BZ[n]=0.0;
            CX[n]=0.0; CY[n]=0.0; CZ[n]=0.0;
            n=n+1;
        }
    }
}
NNi=n-NTi;
NT=n;
NSITE=NT;
printf("Number of atoms=%d \n",NNi);
printf("Total Number of Sites in FCC=%d \n",NSITE);

Atomic Coordinates for nsize = 2

For nsize = 2, the loop generates the following atomic coordinates:

Summary

• The loop generates an FCC structure with nsize = 2, resulting in 32 atoms.
• Atomic coordinates are calculated for both corner and face-centered positions.
• The output lists all atomic coordinates.

Euler Algorithm

The Euler algorithm (or Euler method) is one of the simplest numerical integration methods used to solve ordinary differential equations (ODEs). It is often used in molecular dynamics simulations to update positions and velocities based on forces acting on particles. However, it is known for its low accuracy and lack of time reversibility.

Algorithm Steps

1. Initial Conditions: Given initial positions $r(0)$ and initial velocities $v(0)$.

2. Position Update: The position at the next time step $r(t + \Delta t)$ is calculated using the current position $r(t)$ and velocity $v(t)$:

   $$ r(t + \Delta t) = r(t) + v(t) \Delta t $$

3. Velocity Update: The velocity at the next time step $v(t + \Delta t)$ is calculated using the current velocity $v(t)$ and force $F(t)$:

   $$ v(t + \Delta t) = v(t) + \frac{F(t)}{m} \Delta t $$

Time Irreversibility of the Euler Algorithm

The Euler algorithm is time-irreversible, meaning that if we reverse the time step ($\Delta t \rightarrow -\Delta t$), the algorithm cannot accurately trace back the particle's trajectory. This is due to the first-order nature of the method, which introduces significant numerical errors over time.

Proof of Time Irreversibility

1. Forward Time Step:

   • Position update:

     $$ r(t + \Delta t) = r(t) + v(t) \Delta t $$

   • Velocity update:

     $$ v(t + \Delta t) = v(t) + \frac{F(t)}{m} \Delta t $$

2. Backward Time Step:

   • Replace $\Delta t$ with $-\Delta t$ in the position update:

     $$ r(t - \Delta t) = r(t) - v(t) \Delta t $$

   • Replace $\Delta t$ with $-\Delta t$ in the velocity update:

     $$ v(t - \Delta t) = v(t) - \frac{F(t)}{m} \Delta t $$

3. Verification of Irreversibility:

   • Substitute the forward position update into the backward position update:

     $$ r(t - \Delta t) = r(t) - v(t) \Delta t $$

     This is consistent with the backward position update.

   • However, the velocity update introduces a discrepancy:

     • From the forward step, we have:

       $$ v(t + \Delta t) = v(t) + \frac{F(t)}{m} \Delta t $$

     • If we reverse the time step, we expect:

       $$ v(t) = v(t + \Delta t) - \frac{F(t + \Delta t)}{m} \Delta t $$

     • Substituting the forward velocity update into this equation:

       $$ v(t) = \left(v(t) + \frac{F(t)}{m} \Delta t\right) - \frac{F(t + \Delta t)}{m} \Delta t $$

       Simplifying:

       $$ v(t) = v(t) + \frac{F(t) - F(t + \Delta t)}{m} \Delta t $$

     • For this to hold, we must have:

       $$ F(t) = F(t + \Delta t) $$

       This is generally not true unless the force is constant, which is not the case in most physical systems.

   Thus, the Euler algorithm fails to satisfy time reversibility unless the force is constant.

Conclusion

The Euler algorithm is a simple and easy-to-implement method for numerical integration, but it suffers from low accuracy and time irreversibility. Its first-order nature introduces significant errors over time, making it unsuitable for long-term simulations in molecular dynamics or other systems where energy conservation and time reversibility are critical.

Coding

# Euler算法
for i in range(1, len(t)):
    F = -k * r[i-1]  # 简谐振子的力
    v[i] = v[i-1] + (F / m) * dt  # 更新速度
    r[i] = r[i-1] + v[i-1] * dt    # 更新位置

Verlet Algorithm

The Verlet algorithm is a widely used numerical integration method in molecular dynamics simulations for solving Newton's equations of motion. It updates the positions of particles using their current and previous positions, offering high accuracy and stability.

Algorithm Steps

1. Initial Conditions: Given initial positions $r(0)$ and initial velocities $v(0)$.

2. Position Update: The position at the next time step $r(t + \Delta t)$ is calculated using the current position $r(t)$ and the position at the previous time step $r(t - \Delta t)$:

   $$ r(t + \Delta t) = 2r(t) - r(t - \Delta t) + \frac{F(t)}{m} \Delta t^2 $$

   Here, $F(t)$ is the force acting on the particle at time $t$, $m$ is the particle's mass, and $\Delta t$ is the time step.

3. Velocity Calculation: The velocity can be approximated using the position difference:

   $$ v(t) = \frac{r(t + \Delta t) - r(t - \Delta t)}{2 \Delta t} $$

Proof of Time Reversibility

The time reversibility of the Verlet algorithm means that if we replace the time step $\Delta t$ with $-\Delta t$, the algorithm can correctly trace back the particle's trajectory.

Proof

1. Forward Time Step: According to the Verlet algorithm, the position update for the forward time step is:

   $$ r(t + \Delta t) = 2r(t) - r(t - \Delta t) + \frac{F(t)}{m} \Delta t^2 $$

2. Backward Time Step: Replacing $\Delta t$ with $-\Delta t$, the position update for the backward time step becomes:

   $$ r(t - \Delta t) = 2r(t) - r(t + \Delta t) + \frac{F(t)}{m} \Delta t^2 $$

3. Verification of Reversibility: Substituting the forward time step expression into the backward time step expression:

   $$ r(t - \Delta t) = 2r(t) - \left(2r(t) - r(t - \Delta t) + \frac{F(t)}{m} \Delta t^2\right) + \frac{F(t)}{m} \Delta t^2 $$

   Simplifying this yields:

   $$ r(t - \Delta t) = r(t - \Delta t) $$

   This confirms that the Verlet algorithm remains valid when time is reversed, demonstrating its time reversibility.

Coding

以下是使用Verlet算法更新位置和速度的核心代码片段：

# Verlet算法
for i in range(1, len(t)):
    F = -k * r[i-1]  # 简谐振子的力
    r[i] = 2 * r[i-1] - r[i-2] + (F / m) * dt**2  # 更新位置
    v[i-1] = (r[i] - r[i-2]) / (2 * dt)           # 更新速度

Velocity-Verlet Algorithm

The Velocity-Verlet algorithm is an extension of the Verlet algorithm that explicitly incorporates velocity updates. It is widely used in molecular dynamics simulations due to its simplicity, accuracy, and time-reversibility. The algorithm updates both positions and velocities in a single time step.

Algorithm Steps

1. Initial Conditions: Given initial positions $r(0)$ and initial velocities $v(0)$.

2. Position Update: The position at the next time step $r(t + \Delta t)$ is calculated using the current position $r(t)$, velocity $v(t)$, and force $F(t)$:

   $$ r(t + \Delta t) = r(t) + v(t) \Delta t + \frac{F(t)}{2m} \Delta t^2 $$

3. Velocity Update (Half Step): The velocity is updated using the current force $F(t)$ and the force at the next time step $F(t + \Delta t)$:

   $$ v(t + \Delta t) = v(t) + \frac{F(t) + F(t + \Delta t)}{2m} \Delta t $$

Time Reversibility of the Velocity-Verlet Algorithm

The Velocity-Verlet algorithm is time-reversible, meaning that if we reverse the time step ($\Delta t \rightarrow -\Delta t$), the algorithm can accurately trace back the particle's trajectory.

Proof of Time Reversibility

1. Forward Time Step:

   • Position update:

     $$ r(t + \Delta t) = r(t) + v(t) \Delta t + \frac{F(t)}{2m} \Delta t^2 $$

   • Velocity update:

     $$ v(t + \Delta t) = v(t) + \frac{F(t) + F(t + \Delta t)}{2m} \Delta t $$

2. Backward Time Step:

   • Replace $\Delta t$ with $-\Delta t$ in the position update:

     $$ r(t - \Delta t) = r(t) - v(t) \Delta t + \frac{F(t)}{2m} \Delta t^2 $$

   • Replace $\Delta t$ with $-\Delta t$ in the velocity update:

     $$ v(t - \Delta t) = v(t) - \frac{F(t) + F(t - \Delta t)}{2m} \Delta t $$

3. Verification of Reversibility:

   • Substitute the forward position update into the backward position update:

     $$ r(t - \Delta t) = r(t) - v(t) \Delta t + \frac{F(t)}{2m} \Delta t^2 $$

     This is consistent with the backward position update.

   • Substitute the forward velocity update into the backward velocity update:

     $$ v(t - \Delta t) = v(t) - \frac{F(t) + F(t - \Delta t)}{2m} \Delta t $$

     This is consistent with the backward velocity update.

   Thus, the Velocity-Verlet algorithm satisfies time reversibility.

Conclusion

The Velocity-Verlet algorithm is a powerful and efficient method for molecular dynamics simulations. It explicitly updates both positions and velocities, ensuring high accuracy and stability. Its time-reversibility makes it particularly suitable for simulating physical systems where conservation of energy and momentum is critical.

Coding

# 更新函数
def update(frame):
global positions, velocities

# 计算所有粒子的力
forces = np.zeros_like(positions)
for i in range(num_particles):
for j in range(i + 1, num_particles):
force = lennard_jones_force(positions[i], positions[j])
forces[i] += force
forces[j] -= force

# 更新位置和速度（Velocity-Verlet）
positions += velocities * dt + 0.5 * forces / mass * dt**2
new_forces = np.zeros_like(forces)
for i in range(num_particles):
for j in range(i + 1, num_particles):
force = lennard_jones_force(positions[i], positions[j])
new_forces[i] += force
new_forces[j] -= force
velocities += 0.5 * (forces + new_forces) / mass * dt

# 周期性边界条件
positions = positions % box_size

# 更新散点图
scat._offsets3d = (positions[:, 0], positions[:, 1], positions[:, 2])

The rule to choose the time step

• 根据模拟的目的，时间步长的选择可能会有所不同。例如，对于需要精确轨迹的模拟，可能需要更小的时间步长；而对于只需要热力学性质采样的模拟，可以采用较大的时间步长。The choice of time step can vary depending on the purpose of the simulation. For instance, simulations that require precise trajectories may necessitate a smaller time step, while simulations that only need sampling of thermodynamic properties can employ a larger time step.
• 在实际模拟之前，可以通过缩短模拟时间并检查“守恒量”的数值漂移来验证所选时间步长的合理性。如果漂移过大，则可能需要减小时间步长。Before the actual simulation, the rationality of the selected time step can be verified by shortening the simulation time and checking the numerical drift of the "conserved quantities." If the drift is too large, it may be necessary to reduce the time step.
• For reversible integrators we should not have any drift in energy; for non-reversible ones we want to minimize the drift
• A fluctuation of ~1% in the kinetic energy (in NVE systems) is still acceptable.
• In standard MD with any real material, at ordinary temperatures (> 77 K) it is not possible to have time steps larger than ~ 10 fs.

Cut-off Radius 结算半径

Def
An important factor for the performance of an MD simulation code is its scalability with the number of simulated particles _N_.
If we calculated all the atomic interactions, we would have to carry out ½ N (_N_ – 1) force calculations. (Not good!)
Fortunately, some interatomic forces decrease strongly with distance (van der Waals, covalent interactions etc.) We can thus limit the interactions to be considered within a certain distance, the cut-off radius _r__cut.

Shifting Method

$$ u_{s f}(r)=\left\{\begin{array}{cc} u(r) - u\left(r_c\right) - \frac{d u}{d r}\left(r-r_c\right) & r \leq r_c \\ 0 & r>r_c \end{array}\right. $$

Note that if we do the cut-off for an already parametrized potential we have to ensure that the potential still has the desired properties.

Swift-Off Method

$$ F_s(r)=\left\{\begin{array}{cll} F_\alpha(r), & \text { if } & r<r_1 \\ F_\alpha(r)+S(r), & \text { if } & r_1<r<r_{\mathrm{cut}} \\ 0, & \text { if } & r_{\mathrm{cut}} \leq r \end{array}\right. $$

with the requirement of $S(r)$:

$$ \begin{aligned} S\left(r_1\right) & =0 \\ S^{\prime}\left(r_1\right) & =0 \\ S\left(r_{\mathrm{cut}}\right) & =-F_\alpha\left(r_{\mathrm{cut}}\right) \\ S^{\prime}\left(r_{\mathrm{cut}}\right) & =F_\alpha^{\prime}\left(r_{\mathrm{cut}}\right) \end{aligned} $$

Tersoff-type Potentials

$$ f_C(r)=\left\{\begin{array}{cll} 1, & \text { if } & r \leq R-D \\ \frac{1}{2}-\frac{1}{2} \sin \left(\frac{\pi(r-R)}{2}\right), & \text { if } & R-D<r<R+D \\ 0, & \text { if } & R+D \leq r \end{array}\right. $$

Force Function in MD

Lenard-Jones Potential

$$ U_{\mathrm{LJ}}(r)=4 \varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6\right] $$

12-6 L-J is used in:

• Properties well known (phase behavior etc.)
• soft-matter
• testing of statistical physics
• extremely intensive simulations of FCC metals

Also: 8-4 L-J

Gay-Berne Potential

A variant of L-J Potential used for elongated molecules.

$$ U_{\mathrm{GB}}\left(\vec{u}_i, \vec{u}_j, \vec{s}\right)=4 \varepsilon\left(\vec{u}_i, \vec{u}_j, \vec{s}\right)\left[\left(\frac{\sigma_0}{r-\sigma\left(\vec{u}_i, \vec{u}_j, \vec{s}\right)+\sigma_0}\right)^{12}-\left(\frac{\sigma_0}{r-\sigma\left(\vec{u}_i, \vec{u}_j, \vec{s}\right)+\sigma_0}\right)^6\right] $$

FENE Potential

FENE Potential (Finite Extensible Nonlinear Elastic Potential) is a potential energy function used to describe the interaction between adjacent monomers in polymer chains. It is commonly used to simulate the elastic behavior of polymer chains, especially in molecular dynamics simulations. The characteristic of the FENE Potential is that it is effective within a finite distance and can prevent the monomers in the polymer chain from extending infinitely.

$$ U_{\mathrm{FENE}}=-\frac{1}{2} \kappa R^2 \ln \left[1-\left(\frac{r}{R}\right)^2\right] $$

! Loop over all pairs of particles
do i = 1, N-1
   do j = i+1, N
      Sij = pos(:,i) - pos(:,j)         ! distance vector between i j
      where ( abs(Sij) > 0.5d0 )        ! (in box scaled units)
         Sij = Sij - sign(1.d0, Sij)    ! periodic boundary conditions
      end where                         ! applied where needed.
      Rij = BoxSize*Sij                 ! go to real space units
      Rsqij = dot_product(Rij, Rij)     ! compute square distance
      if ( Rsqij < Rcutoff**2 ) then    ! particles are interacting?
         ! compute Lennard-Jones potntl
         rm2 = 1.d0/Rsqij               ! 1/r^2
         rm6 = rm2**3                   ! 1/r^6
         rm12 = rm6**2                  ! 1/r^12
         phi = 4.d0 * ( rm12 - rm6 ) - phicutoff  ! 4[1/r^12 - 1/r^6] - phi(Rc)
         dphi = 24.d0*rm2*( 2.d0*rm12 - rm6 )    ! The following is dphi = -(1/r)(dV/dr)
         ene_pot(i) = ene_pot(i) + 0.5d0*phi    ! accumulate energy
         ene_pot(j) = ene_pot(j) + 0.5d0*phi    ! (i and j share it)
         virial = virial - dphi*Rsqij           ! accumul. virial=sum r(dV/dr)
         acc(:,i) = acc(:,i) + dphi*Sij         ! accumulate forces
         acc(:,j) = acc(:,j) - dphi*Sij         ! (Fji = -Fij)
      endif
   enddo
enddo

这段Fortran代码实现了一个分子动力学模拟中的粒子间相互作用计算，主要涉及Lennard-Jones势的计算和周期性边界条件的处理。以下是代码的详细解释：

1. 循环遍历所有粒子对

do i = 1, N-1
   do j = i+1, N

• 外层循环遍历所有粒子，内层循环遍历当前粒子之后的所有粒子。这样可以避免重复计算粒子对 $(i, j)$ 和 $(j, i)$。

2. 计算粒子间的距离向量

Sij = pos(:,i) - pos(:,j)

• Sij 是粒子 $i$ 和粒子 $j$ 之间的距离向量。pos(:,i) 表示粒子 $i$ 的位置向量。

3. 周期性边界条件处理

where ( abs(Sij) > 0.5d0 )
   Sij = Sij - sign(1.d0, Sij)
end where

• 如果距离向量的某个分量绝对值大于 $0.5$（在盒子缩放单位中），则应用周期性边界条件，将距离向量调整到最近的镜像位置。

4. 转换到实际空间单位

Rij = BoxSize*Sij

• 将缩放单位下的距离向量 Sij 转换为实际空间单位 Rij，BoxSize 是盒子的尺寸。

5. 计算距离的平方

Rsqij = dot_product(Rij, Rij)

• 计算粒子间距离的平方 Rsqij，即 $R_{ij}^2$。

6. 判断粒子是否在相互作用范围内

if ( Rsqij < Rcutoff**2 ) then

• 如果粒子间的距离平方小于截断距离 $R_{\text{cutoff}}^2$，则认为粒子间存在相互作用。

7. 计算Lennard-Jones势

rm2 = 1.d0/Rsqij
rm6 = rm2**3
rm12 = rm6**2
phi = 4.d0 * ( rm12 - rm6 ) - phicutoff

• rm2 是 $1/r^2$，rm6 是 $1/r^6$，rm12 是 $1/r^{12}$。

• phi 是Lennard-Jones势能，计算公式为：

  $$ \phi = 4 \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right) - \phi_{\text{cutoff}} $$

  其中 $\phi_{\text{cutoff}}$ 是截断距离处的势能。

8. 计算势能的导数

dphi = 24.d0*rm2*( 2.d0*rm12 - rm6 )

• dphi 是势能的导数，计算公式为：

  $$ \frac{d\phi}{dr} = 24 \left( \frac{2}{r^{12}} - \frac{1}{r^6} \right) $$

9. 累加势能

ene_pot(i) = ene_pot(i) + 0.5d0*phi
ene_pot(j) = ene_pot(j) + 0.5d0*phi

• 将势能 phi 的一半分别累加到粒子 $i$ 和粒子 $j$ 的势能中，因为势能是共享的。

10. 累加维里

virial = virial - dphi*Rsqij

• 累加维里项，计算公式为：

  $$ \text{virial} = \text{virial} - \frac{d\phi}{dr} \cdot r^2 $$

11. 累加力

acc(:,i) = acc(:,i) + dphi*Sij
acc(:,j) = acc(:,j) - dphi*Sij

• 将力 dphi*Sij 累加到粒子 $i$ 的加速度中，并从粒子 $j$ 的加速度中减去，因为 $F_{ji} = -F_{ij}$。

总结

这段代码实现了分子动力学模拟中粒子间Lennard-Jones势的计算，包括周期性边界条件的处理、势能和力的累加，以及维里项的计算。

Expression of L-J, F-S and Bound-order Potential

Here are the mathematical expressions of the four potential energy functions:

1. Lennard-Jones Potential (LJ Potential)

The Lennard-Jones potential is used to describe the attractive and repulsive interactions between molecules, commonly applied in simulations of non-polar molecular systems.

$$ V(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right] $$

• $r$: Distance between particles
• $\epsilon$: Depth of the potential well (energy parameter)
• $\sigma$: Distance at which the potential is zero (length parameter)

2. Embedded Atom Method Potential (EAM Potential)

The EAM potential is used to describe interatomic interactions in metallic systems, taking into account the embedding effect of electron clouds.

$$ E_i = F_\alpha \left( \sum_{j \neq i} \rho_\beta (r_{ij}) \right) + \frac{1}{2} \sum_{j \neq i} \phi_{\alpha \beta} (r_{ij}) $$

• $E_i$: Total energy of the $i$-th atom
• $F_\alpha$: Embedding function, describing the contribution of the electron cloud to the atomic energy
• $\rho_\beta$: Electron density function of the $j$-th atom
• $\phi_{\alpha \beta}$: Pair potential function, describing the interaction between atoms $i$ and $j$
• $r_{ij}$: Distance between atoms $i$ and $j$

3. Finnis-Sinclair Potential (F-S Potential)

The Finnis-Sinclair potential is a simplified form of the EAM potential, often used to describe interactions in transition metals.

$$ E_i = \frac{1}{2} \sum_{j \neq i} V(r_{ij}) - A \sqrt{\sum_{j \neq i} \rho(r_{ij})} $$

• $E_i$: Total energy of the $i$-th atom
• $V(r_{ij})$: Pair potential function
• $\rho(r_{ij})$: Electron density function
• $A$: Constant parameter
• $r_{ij}$: Distance between atoms $i$ and $j$

4. Bond-Order Potential (BOP Potential)

The Bond-Order Potential is used to describe interatomic interactions in covalent systems, taking into account the influence of bond order.

$$ E_i = \frac{1}{2} \sum_{j \neq i} \left[ V_R(r_{ij}) - b_{ij} V_A(r_{ij}) \right] $$

• $E_i$: Total energy of the $i$-th atom
• $V_R(r_{ij})$: Repulsive potential function
• $V_A(r_{ij})$: Attractive potential function
• $b_{ij}$: Bond-order parameter, describing the strength of the bond between atoms $i$ and $j$
• $r_{ij}$: Distance between atoms $i$ and $j$

These potential energy functions are widely used in molecular dynamics simulations and materials science, with the specific choice depending on the physical characteristics of the system.

Summary Comparison:

When selecting a potential function, it is necessary to weigh the specific properties of the system (e.g., metallic, covalent, or van der Waals interactions) and computational requirements (e.g., accuracy and efficiency).

Temperature & Pressure

• Pressure. Generalized equipartition theorem (Hamiltonian formalism)

$$\langle q_k \frac{\partial H}{\partial q_k} \rangle = k_B T$$

In cartesian coordinates:

$$-\frac{1}{3} \langle \sum_i \vec{r}_i \cdot (\nabla_{r_i} U) \rangle = \frac{1}{3} \langle \sum_i \vec{r}_i \cdot \vec{f}_{tot} \rangle = -N k_B T$$

Then, let us divide the total force as

$$\vec{f}_{tot} = \vec{f}_{ext} + \vec{f}_i$$

Whence we get

$$\frac{1}{3} \langle \sum_i \vec{r}_i \cdot \vec{f}_{ext} \rangle = -PV \quad \text{and} \quad W = -\frac{1}{3} \sum_i \vec{r}_i \cdot \vec{f}_i$$

Giving us the relation to pressure finally as

$$PV = N k_B T + \langle W \rangle$$

Temp

Temperature Characterization: Calculated by distribution, not averages.

Anderson

Anderson Thermostat 是一种用于分子动力学模拟的恒温方法，由 H.C. Andersen 在 1980 年提出。它的主要目的是通过随机碰撞来维持系统的温度恒定，从而模拟恒温条件下的粒子运动。

基本原理

Anderson Thermostat 的核心思想是通过随机选择粒子并重新赋予它们新的速度，使其符合给定的温度分布（通常是 Maxwell-Boltzmann 分布）。具体步骤如下：

1. 随机选择粒子：在模拟过程中，每隔一定的时间步长，随机选择一部分粒子。
2. 重新赋予速度：为选中的粒子重新赋予速度，这些速度从 Maxwell-Boltzmann 分布中随机采样，该分布对应于给定的温度 $T$。
3. 维持温度：通过这种方式，系统的温度被维持在设定值附近。

数学描述

假设系统的温度为 $T$，粒子的质量为 $m$，则 Maxwell-Boltzmann 分布下的速度分量 $v$ 的概率密度函数为：

$$ P(v) = \sqrt{\frac{m}{2\pi k_B T}} \exp\left(-\frac{m v^2}{2 k_B T}\right) $$

其中，$k_B$ 是玻尔兹曼常数。

优点

• 简单易实现：Anderson Thermostat 的实现相对简单，计算开销较小。
• 快速达到平衡：通过随机碰撞，系统可以快速达到热平衡。

缺点

• 物理不真实性：由于速度是随机重新赋予的，这种方法可能会破坏系统的动力学行为，尤其是在研究时间相关性质时。
• 参数依赖性：需要选择合适的碰撞频率，这可能会影响模拟结果。

Berenden Thermostat

The weak-coupling thermostat, commonly known as the Berendsen thermostat or temperature control, is a common way for controlling the temperature in MD simulations. (J. Chem. Phys. 81 (1984) 3684)

The basic idea is that we add a frictional term to the equations of motion,

$$ m_i \frac{dv_i}{dt} = f_i + m_i \gamma \left( \frac{T_0}{T} - 1 \right) v_i $$

which then drives the system (exponentially) toward the desired temperature. In practice, the velocities of the particles are scaled with a factor

$$ \lambda = \sqrt{1 + \frac{\delta t}{\tau_T} \left( \frac{T_0}{T} - 1 \right)} $$

$$ \frac{dT(t)}{dt} = \frac{1}{\tau} (T_0 - T(t)) $$

$$ \Delta T = \frac{\delta t}{\tau} (T_0 - T(t)) $$

The Berendsen algorithm is simple to implement and it is very efficient for reaching the desired temperature from far-from-equilibrium configurations.

... The downside is that it has been argued that this method does not produce the correct statistical ensemble.

Nosé-Hoover thermostat

Another popular thermostat is the Nosé-Hoover thermostat, which retains a Boltzmann equilibrium distribution. (Hoover, Phys. Rev. A 31 (1985) 1695)

Again, we introduce a frictional term to the equations of motion

$$ \frac{dp}{dt} = F(q, t) - \zeta(t)p(t) $$

with the dynamics of the friction coefficient by

$$ Q \frac{d\zeta}{dt} = \sum_i \left( \frac{1}{2} m_i v_i^2 - \frac{3}{2} N k_B T \right) $$

(Feedback makes kinetic energy = $3/2 k_B T$)

Q = fictional "heat bath mass". Large Q means weak coupling.

Nosé suggested $Q \sim g k_B T$, where $g$ is the number of degrees of freedom in the system ($\sim 6N$).

Properties

1. Averages
   物性参数可以根据原子的坐标和速度通过统计处理得出，在统计物理中可以利用系综微观量的统计平均值来计算物性参数值，即

   $$ \langle A \rangle_{\text{ens}} = \frac{1}{N!} \int \int dp dr A(p, r) p(p, r) $$

   式中 $ens$ 表示系综 (Ensemble)，
   $p(p, r)$ 表示分布函数 (几率密度)。分布函数因所使用的系综不同而不同。
   在图中，公式和文本如下：

   $$ P_{\text{NTV}}(p, r) = \exp \left[ -\frac{H(p, r)}{k_B T} \right] $$

   $$ P_{\text{NPT}}(p, r) = \exp \left[ -\frac{H(p, r) + PV}{k_B T} \right] $$

   $$ P_{\text{NTV}}(p, r) = \exp \left[ -\frac{H(p, r) - \mu N}{k_B T} \right] $$

   关于分子动力学中的时间平均和系综平均：时间平均等于系综平均的各态历经假设，即：

   $$ \langle A \rangle_{\text{ens}} = \langle A \rangle_{\text{time}} $$

   时间平均的计算公式为：

   $$ \langle A \rangle_{\text{time}} = \lim_{t_{\text{obs}} \to \infty} \frac{1}{t_{\text{obs}}} \int A(p, r) dt = \frac{1}{t_{\text{obs}}} \sum A(p, r) $$

2. Static Properties
   (1) 温度 ( T )

• 动能公式：

$$E_K = \sum_{i=1}^{N} \frac{|p_i|^2}{2m_i} = \frac{k_B T}{2} (3N - N_C)$$

• 温度公式：

$$T = \frac{1}{(3N - N_C) K_B} \sum_{i=1}^{N} \frac{|p_i|^2}{m_i}$$

(2) 能量

• 总能量公式：

$$E = \frac{1}{2} \sum_i m_i V_i^2 + \sum_i U_i$$

• 静态性能分析：

$$U = \sum_i U_i$$

• 动能公式：

$$E_K = \frac{1}{2} \sum_i m_i V_i^2$$

(3) 压力
压力通常通过虚功原理模拟得到。虚功定义为所有粒子坐标和作用在粒子上的力的乘积的和，通常写为：

• 虚功公式：

$$W = -3PV + \sum_{i=1}^{N} \sum_{j=i+1}^{N} r_{ij} \frac{dU(r_{ij})}{dr_{ij}} = -3Nk_B T$$

实际体系的虚功为理想气体的虚功与粒子之间相互作用部分的虚功的和。

• 压力公式：

$$P = \frac{1}{V} \left[ N k_B T - \frac{1}{3} \sum_{i=1}^{N} \sum_{j=i+1}^{N} r_{ij} f_{ij} \right]$$

(4) 径向分布函数 PCF
(5) 静态结构因子 Static SF
静态结构因子 (Static Structure Factor) 也是判断结构无序程度的物理量。它的表示式为：

$$ S(K) = \frac{1}{N} \left| \sum_{j=1}^{N} \exp(i \mathbf{K} \cdot \mathbf{r}_j) \right| $$

其中，$N$ 代表总原子数，$\mathbf{K}$ 为倒格矢，$\mathbf{r}_j$ 为原子 $j$ 的位置矢量。
对理想晶体：其静态结构因子为 1。
而对理想流体：为 0。静态结构因子在研究晶体的熔化与相变的研究中很有用。
(6) 热力学性质 Thermodynamics Properties
在相变时，比热会呈现与温度相关的特征。对一级相变点，比热呈现无限大；对二级相变点，比热呈现不连续变化。监控比热随温度的变化可以帮助我们探测到相变的发生。

$$ C_V = \left( \frac{\partial U}{\partial T} \right)_V $$

$$ C_V = \frac{\langle (E - \langle E \rangle)^2 \rangle}{k_B T^2} $$

对等温等压系综 (NPT)

$$ C_P = \frac{\langle \delta (E_K + E_p + PV)^2 \rangle}{k_B T^2} $$

$$ \beta_T = \frac{\langle \delta V^2 \rangle}{k_B T V} $$

$$ \alpha_P = \frac{\langle \delta V \delta (E_K + E_p + PV) \rangle}{k_B T^2 V} $$

常压比热、等温压缩率和热膨胀系数。
(7) 配位数 Coordination Number

5. Dynamic Perf

   1. 关联函数
   2. 运输性质
   3. 扩散系数 Diffusion Coefficient
   4. 热传导系数 Thermal conductivity coefficient
   5. 热膨胀系数 Thermal Expansion Coefficient