提问人：user366312 提问时间：9/10/2023 最后编辑：user366312 更新时间：9/13/2023 访问量：93
跨越式集成以更新位置和速度
Leapfrog integration to update position and velocity
问：
我写了这个最小的例子来研究 Leapfrog 集成算法。
但是，我不确定它是否是正确的算法，并且列表给出了正确的输出。
这是 Leapfrog 算法吗？
注意：我正在使用 Lennard-Jones 系统。
#include "stdafx.h"
#include <cmath>
#include <iostream>
class Vec3 {
public:
Vec3() {}
Vec3(double x, double y, double z) : x(x), y(y), z(z) {}
Vec3(const Vec3& other) : x(other.x), y(other.y), z(other.z) {} // Copy constructor
Vec3& operator=(const Vec3& other) {
if (this != &other) {
x = other.x;
y = other.y;
z = other.z;
}
return *this;
} // Copy assignment operator
Vec3 operator+(const Vec3& b) const {
return Vec3(x + b.x, y + b.y, z + b.z);
}
Vec3 operator*(double scalar) const {
return Vec3(scalar * x, scalar * y, scalar * z);
}
// Friend function for scalar multiplication with scalar on the left-hand side
friend Vec3 operator*(double scalar, const Vec3& vector) {
return vector * scalar;
}
double getX() const { return x; }
double getY() const { return y; }
double getZ() const { return z; }
static Vec3 Zero() {
return Vec3(0.0, 0.0, 0.0);
}
private:
double x, y, z;
};
class Atom {
public:
Atom() {}
Atom(const Vec3& position, const Vec3& velocity)
: position(position), velocity(velocity) {}
Atom(const Atom& other)
: position(other.position), velocity(other.velocity) {}
// Copy constructor
Atom& operator=(const Atom& other) {
if (this != &other) {
position = other.position;
velocity = other.velocity;
}
return *this;
}
// Copy assignment operator
Vec3 getPosition() const { return position; }
Vec3 getVelocity() const { return velocity; }
void setPosition(const Vec3& position) { this->position = position; }
void setVelocity(const Vec3& velocity) { this->velocity = velocity; }
private:
Vec3 position;
Vec3 velocity;
};
class ArgonGasSimulation {
private:
static constexpr double Sigma = 3.4;
static constexpr double Epsilon = 1.0;
static Vec3 ComputeForce(const Vec3& position) {
double distance = sqrt(pow(position.getX(), 2) + pow(position.getY(), 2) + pow(position.getZ(), 2));
double r6 = pow(distance, 6);
double r12 = pow(r6, 2);
double magnitude = 48 * Epsilon * (pow((Sigma / distance), 6) - 0.5 * pow((Sigma / distance), 4)) / distance;
return position * magnitude;
}
static void Leapfrog(Atom& atom, double timeStep) {
Vec3 force = ComputeForce(atom.getPosition()); // Compute force based on the current position
Vec3 newPosition = atom.getPosition() + atom.getVelocity() * timeStep + 0.5 * force * pow(timeStep, 2);
Vec3 newForce = ComputeForce(newPosition);
Vec3 newVelocity = atom.getVelocity() + (force + newForce) * timeStep * 0.5;
// Update state
atom.setPosition(newPosition);
atom.setVelocity(newVelocity);
}
public:
static void Main() {
// Initial state
Vec3 position(1.0, 2.0, 3.0); // initial position
Vec3 velocity = Vec3::Zero(); // initial velocity
double timeStep = 0.01;
Atom atom(position, velocity);
// Leapfrog integration
Leapfrog(atom, timeStep);
// Output new position and velocity
std::cout << "New Position: (" << atom.getPosition().getX() << ", " << atom.getPosition().getY() << ", " << atom.getPosition().getZ() << ")" << std::endl;
std::cout << "New Velocity: (" << atom.getVelocity().getX() << ", " << atom.getVelocity().getY() << ", " << atom.getVelocity().getZ() << ")" << std::endl;
}
};
int main() {
ArgonGasSimulation::Main();
return 0;
}
输出：
New Position: (1.00014244382518, 2.00028488765035, 3.00042733147553)
New Velocity: (0.028470372373851, 0.0569407447477019, 0.0854111171215529)
Press any key to continue . . .
编辑：以下内容是根据 Lutz Lehman 的回答用 C++ 重写的。
#include "stdafx.h"
#include <iostream>
#include <vector>
#include <cmath>
#include <random>
// Constants for Argon
constexpr double epsilon = 119.8; // Depth of the potential well (in K)
constexpr double sigma = 3.405; // Distance for zero potential (in Angstrom)
constexpr double mass = 39.948; // Mass of Argon (in amu)
struct Vector3D {
double x, y, z;
};
// Lennard-Jones potential function
double lj_potential(const Vector3D& r)
{
double r_mag = std::sqrt(r.x * r.x + r.y * r.y + r.z * r.z);
double s_over_r = sigma / r_mag;
double s_over_r6 = pow(s_over_r, 6);
return 4.0 * epsilon * (s_over_r6 * s_over_r6 - s_over_r6);
}
// Derivative of the Lennard-Jones potential
Vector3D lj_force(const Vector3D& r)
{
// Define a small distance for the derivative approximation
double dr = 1e-6;
Vector3D force;
std::vector<Vector3D> r_plus_dr = { r, r, r };
r_plus_dr[0].x += dr;
r_plus_dr[1].y += dr;
r_plus_dr[2].z += dr;
// The force is the negative derivative of the potential energy
force.x = -(lj_potential(r_plus_dr[0]) - lj_potential(r)) / dr;
force.y = -(lj_potential(r_plus_dr[1]) - lj_potential(r)) / dr;
force.z = -(lj_potential(r_plus_dr[2]) - lj_potential(r)) / dr;
return force;
}
// Update the 'accel' function
void accel(std::vector<Vector3D>& a, const std::vector<Vector3D>& x)
{
int n = x.size(); // number of points
for (int i = 0; i < n; i++)
{
auto force = lj_force(x[i]);
// use Lennard-Jones force law
a[i].x = -force.x / mass;
a[i].y = -force.y / mass;
a[i].z = -force.z / mass;
}
}
void leapfrog_step(std::vector<Vector3D>& x, std::vector<Vector3D>& v, double dt) {
int n = x.size(); // number of points
std::vector<Vector3D> a(n);
accel(a, x);
for (int i = 0; i < n; i++)
{
v[i].x = v[i].x + 0.5 * dt * a[i].x; // advance vel by half-step
v[i].y = v[i].y + 0.5 * dt * a[i].y; // advance vel by half-step
v[i].z = v[i].z + 0.5 * dt * a[i].z; // advance vel by half-step
x[i].x = x[i].x + dt * v[i].x; // advance pos by full-step
x[i].y = x[i].y + dt * v[i].y; // advance pos by full-step
x[i].z = x[i].z + dt * v[i].z; // advance pos by full-step
}
accel(a, x);
for (int i = 0; i < n; i++)
{
v[i].x = v[i].x + 0.5 * dt * a[i].x; // and complete vel. step
v[i].y = v[i].y + 0.5 * dt * a[i].y; // and complete vel. step
v[i].z = v[i].z + 0.5 * dt * a[i].z; // and complete vel. step
}
}
// Initialize positions and velocities of particles
void initialize(std::vector<Vector3D>& x, std::vector<Vector3D>& v, int n_particles, double box_size, double max_vel) {
// Create a random number generator
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(-0.5, 0.5);
// Resize the vectors to hold the positions and velocities of all particles
x.resize(n_particles);
v.resize(n_particles);
// Initialize positions and velocities
for (int i = 0; i < n_particles; i++)
{
// Assign random initial positions within the box
x[i].x = box_size * distribution(generator);
x[i].y = box_size * distribution(generator);
x[i].z = box_size * distribution(generator);
// Assign random initial velocities up to max_vel
v[i].x = max_vel * distribution(generator);
v[i].y = max_vel * distribution(generator);
v[i].z = max_vel * distribution(generator);
}
}
int main() {
int n_particles = 100; // number of particles
double box_size = 10.0; // size of the simulation box
double max_vel = 0.1; // maximum initial velocity
double dt = 0.01; // time step
int n_steps = 10; // number of time steps
// Positions and velocities of the particles
std::vector<Vector3D> x, v;
// Initialize the particles
initialize(x, v, n_particles, box_size, max_vel);
// Run the simulation
for (int step = 0; step < n_steps; step++) {
leapfrog_step(x, v, dt);
}
return 0;
}
编辑-2：在阅读了 Lutz Lehman 帖子下的评论后，我写了以下列表。
#include "stdafx.h"
#include <iostream>
#include <vector>
#include <cmath>
#include <random>
// Constants for Argon
constexpr double epsilon = 119.8; // Depth of the potential well (in K)
constexpr double sigma = 3.405; // Distance for zero potential (in Angstrom)
constexpr double mass = 39.948; // Mass of Argon (in amu)
struct Vector3D {
double x, y, z;
};
// Derivative of the Lennard-Jones potential
Vector3D lj_force(const Vector3D& r)
{
double r_mag = std::sqrt(r.x * r.x + r.y * r.y + r.z * r.z);
double s_over_r = sigma / r_mag;
double s_over_r6 = s_over_r * s_over_r * s_over_r * s_over_r * s_over_r * s_over_r;
double s_over_r12 = s_over_r6 * s_over_r6;
double factor = 24.0 * epsilon * (2.0 * s_over_r12 - s_over_r6) / (r_mag * r_mag * r_mag);
Vector3D force;
force.x = factor * r.x;
force.y = factor * r.y;
force.z = factor * r.z;
return force;
}
// Update the 'accel' function
void accel(std::vector<Vector3D>& a, const std::vector<Vector3D>& x)
{
int n = x.size(); // number of points
// Reset the acceleration to zero
for (int i = 0; i < n; i++)
{
a[i].x = 0.0;
a[i].y = 0.0;
a[i].z = 0.0;
}
// Compute the acceleration due to each pair of particles
for (int i = 0; i < n; i++)
{
for (int j = i + 1; j < n; j++)
{
Vector3D r;
r.x = x[j].x - x[i].x;
r.y = x[j].y - x[i].y;
r.z = x[j].z - x[i].z;
Vector3D force = lj_force(r);
// use Lennard-Jones force law
a[i].x += force.x / mass;
a[i].y += force.y / mass;
a[i].z += force.z / mass;
a[j].x -= force.x / mass;
a[j].y -= force.y / mass;
a[j].z -= force.z / mass;
}
}
}
void leapfrog_step(std::vector<Vector3D>& x, std::vector<Vector3D>& v, double dt) {
int n = x.size(); // number of points
std::vector<Vector3D> a(n);
accel(a, x);
for (int i = 0; i < n; i++)
{
v[i].x = v[i].x + 0.5 * dt * a[i].x; // advance vel by half-step
v[i].y = v[i].y + 0.5 * dt * a[i].y; // advance vel by half-step
v[i].z = v[i].z + 0.5 * dt * a[i].z; // advance vel by half-step
x[i].x = x[i].x + dt * v[i].x; // advance pos by full-step
x[i].y = x[i].y + dt * v[i].y; // advance pos by full-step
x[i].z = x[i].z + dt * v[i].z; // advance pos by full-step
}
accel(a, x);
for (int i = 0; i < n; i++)
{
v[i].x = v[i].x + 0.5 * dt * a[i].x; // and complete vel. step
v[i].y = v[i].y + 0.5 * dt * a[i].y; // and complete vel. step
v[i].z = v[i].z + 0.5 * dt * a[i].z; // and complete vel. step
}
}
// Initialize positions and velocities of particles
void initialize(std::vector<Vector3D>& x, std::vector<Vector3D>& v, int n_particles, double box_size, double max_vel) {
// Create a random number generator
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(-0.5, 0.5);
// Resize the vectors to hold the positions and velocities of all particles
x.resize(n_particles);
v.resize(n_particles);
// Initialize positions and velocities
for (int i = 0; i < n_particles; i++)
{
// Assign random initial positions within the box
x[i].x = box_size * distribution(generator);
x[i].y = box_size * distribution(generator);
x[i].z = box_size * distribution(generator);
// Assign random initial velocities up to max_vel
v[i].x = max_vel * distribution(generator);
v[i].y = max_vel * distribution(generator);
v[i].z = max_vel * distribution(generator);
}
}
int main() {
int n_particles = 100; // number of particles
double box_size = 10.0; // size of the simulation box
double max_vel = 0.1; // maximum initial velocity
double dt = 0.01; // time step
int n_steps = 10; // number of time steps
// Positions and velocities of the particles
std::vector<Vector3D> x, v;
// Initialize the particles
initialize(x, v, n_particles, box_size, max_vel);
// Run the simulation
for (int step = 0; step < n_steps; step++) {
leapfrog_step(x, v, dt);
}
return 0;
}
答：
1赞
Lutz Lehmann
9/11/2023
#1
关于为什么代码没有实现 transfrog Verlet 以及它应该是什么的问题
标准 Verlet 方法变换中心二阶差分商
(x(t+h)-2x(t)+x(t-h))/h^2 = x''(t) + O(h^2) = a(x(t)) + O(h^2)
变成时间传播器
x[i+1] = 2*x[i]-x[i-1] + h^2*a[i]
固有的双精度求和可以使用辅助变量拆分为两个链式求和
h*v[i+0.5] = x[i+1]-x[i]
这些起初只是上述序列的点序列的差商或割线斜率。它们最接近解的导数，它们是中心差商，即在中点。简化的 Verlet 递归为t[i+0.5]=t[i]+0.5*h
h*a[i] = v[i+0.5] - v[i-0.5]
或作为算法
a[i] = acc(x[i])
v[i+0.5] = v[i-0.5] + h*a[i]
x[i+1] = x[i] + h*v[i+0.5]
将这种交错网格中的值组合在一起的方法被赋予了绰号“跳跃”或“棋盘”。
为了获得原始网格点处位置和速度的完整状态，下一个最佳中心差商是t[i]
2*h*v[i] = x[i+1]-x[i-1] = h*(v[i+0.5]+v[i-0.5])
现在有几种方法可以获得这些速度的阶跃方程。它们都以
v[i+1]-v[i] = 0.5*h*(a[i+1]+a[i])
作为显式步骤中的算法，这被实现为
v[i+0.5] = v[i] + 0.5*h*a[i]
x[i+1] = x[i] + h*v[i+0.5]
a[i+1] = acc(x[i+1])
v[i+1] = v[i+0.5] + 0.5*h*a[i+1]
这就是现在完全形成的辛方法。作为一个简称，它被称为“速度”Verlet 方法。
通过应用半步的时移来获得它的变体，因此现在主状态序列是以前半步的状态。这给出了不同的点序列，但差异低于方法误差，因此没有实际意义
x[i+0.5] = x[i] + 0.5*h*v[i]
a[i+0.5] = acc(x[i+0.5])
v[i+1] = v[i] + h*a[i+0.5]
x[i+1] = x[i+0.5] + 0.5*h*v[i+1]
加速度计算和应用现在可以看作是一个单元，加速度矢量只与当前步骤相关，不需要为下一步保存。
评论
0赞
user366312
9/13/2023
请检查编辑。
0赞
Lutz Lehmann
9/13/2023
是的，这也是正确的。如果在步骤上保持加速度矢量，则第一个加速度计算将是多余的。请注意，这为您提供了在原点周围的一个固定中心场中具有 LJ 势的独立粒子。这不是分子模型的工作方式。请参阅带有图形输出的 stackoverflow.com/questions/29374247/...。n
0赞
user366312
9/13/2023
请注意，这为您提供了 n 个独立粒子，位于原点周围的一个固定中心场中，具有 LJ 势。这不是分子模型的工作方式。--- 在这方面，我怎样才能修复我的C++代码？（注意：我不精通 Python。
0赞
Lutz Lehmann
9/13/2023
您必须遍历所有货币对才能计算所有相互作用力。力总是取决于位置差异，而不是位置本身。
0赞
user366312
9/13/2023
它正确吗？lj_force(()
评论