Last week, I posted a simple Lennard–Jones simulation, written in C++, that models the behavior of liquid Ar in only 1561 characters. By popular request, I'm now posting a breakdown/explanation of this program, both to explain the underlying algorithm and illustrate how it can be compressed.
Here's the most current version of the program, now condensed even further to 1295 characters:
#include <iostream>
#include <random>
#define H return
typedef double d;typedef int i;using namespace std;i N=860;d B=8.314e-7;d S=1544
.8;d q(d x){if(x>17)x-=34;if(x<-17)x+=34;H x;}d rd(){H 34*drand48()-17;}struct v
{d x,y,z;v(d a=0,d b=0,d c=0):x(a),y(b),z(c){}d n(){H x*x+y*y+z*z;}v p(){H v(q(x
),q(y),q(z));}};v operator+(v a,v b){H v(a.x+b.x,a.y+b.y,a.z+b.z);}v operator*(d
a,v b){H v(a*b.x,a*b.y,a*b.z);}struct p{v r,q,a,b;vector<i> W;p(v l):r(l){}void
F(v f){b=0.025*f+b;}};vector<p> P;void w(){cout<<N<<"\n\n";for(p l:P)cout<<"Ar
"<<l.r.x<<" "<<l.r.y<<" "<<l.r.z<<"\n";cout<<"\n";}void E(){for(i j=N;j--;){for(
i k:P[j].W){v R=(P[j].r+-1*P[k].r).p();d r=R.n();if(r<70){d O=r*r*r;R=2880*B*(2*
S*S/O/O-S/O)/r*R;P[j].F(R);P[k].F(-1*R);}}}}void e(){for(i j=0;j<N;j++){P[j].W.c
lear();for(i k=j+1;k<N;k++)if((P[j].r+-1*P[k].r).p().n()<90)P[j].W.push_back(k);
}}i main(){i A=1e3;i Y=5e3;for(i j=N;j--;){for(i a=999;a--;){v r=v(rd(),rd(),rd(
));i c=0;for(p X:P){d D=(r+-1*X.r).p().n();if(D<6.8)c=1;}if(!c){P.push_back(p(r)
);break;}}}for(i t=0;t<=3e4;t++){for(p& I:P)I.r=(I.r+I.q+0.5*I.a).p();if(t%20==0
)e();E();d K=0;for(p& I:P){I.q=I.q+0.5*(I.b+I.a);I.a=I.b;I.b=v();K+=20*I.q.n();}
d T=2*K/(3*B*N);if(t<2*Y){d C=75;if(t<Y)C=75*t/Y+(A-75)*(Y-t)/Y;for(p& I:P)I.q=s
qrt(C/T)*I.q;}if(t%100==0&&t>2*Y)w();}}
Let's add some whitespace and break this down, line by line:
#include <iostream> #include <random> // some abbreviations #define H return typedef double d; typedef int i; using namespace std; // constants i N=860; // number of particles d B=8.314e-7; // boltzmann's constant d S=1544.8; // the minimum of the Lennard-Jones potential, 3.4 Å, to the 6th power
The program begins by importing the necessary packages, defining abbreviations, and declaring some constants. Redefining double and int saves a ton of space, as we'll see. (This is standard practice in the code abbreviation world.)
// given a 1d coordinate, scale it to within [-17,17].
// (this only works for numbers within [-51,51] but that's fine for this application)
d q(d x){
if(x>17)
x-=34;
if(x<-17)
x+=34;
H x;
}
// returns a uniform random number in [-17,17]
d rd(){
H 34*drand48()-17;
}
We now define a helper function to keep coordinates within the boundaries of our cubical box, and create another function to "randomly" initialize particles' positions. (drand48 is not a particularly good random-number generator, but it has a short name and works well enough.)
// vector class
struct v{
d x,y,z;
v(d a=0,d b=0,d c=0):x(a),y(b),z(c){}
// return the squared length of the vector
d n(){
H x*x+y*y+z*z;
}
// return a vector with periodic boundary conditions applied
v p(){
H v(q(x),q(y),q(z));
}
};
// vector addition
v operator+(v a,v b){
H v(a.x+b.x,a.y+b.y,a.z+b.z);
}
// multiplication by a scalar
v operator*(d a,v b){
H v(a*b.x,a*b.y,a*b.z);
}
Here, we define a vector class to store positions, velocities, and accelerations, and define addition and multiplication by a scalar. (It would be better to pass const references to the operators, but it takes too many characters.)
// particle class
struct p{
// r is position, q is velocity, a is acceleration, b is acceleration in the next timestep
v r,q,a,b;
// neighbor list, the list of particles close enough to consider computing interactions with
vector<i> W;
p(v l):r(l){}
// apply a force to this particle.
// 0.025 = 1/40 = 1/(Ar mass)
void F(v f){
b=0.025*f+b;
}
};
// global vector of all particles
vector<p> P;
// write current coordinates to stdout
void w(){
cout<<N<<"\n\n";
for(p l:P)
cout<<"Ar "<<l.r.x<<" "<<l.r.y<<" "<<l.r.z<<"\n";
cout<<"\n";
}
Now, we define a class that represents a single argon atom. Each atom has an associated position, velocity, and acceleration, as well as b, which accumulates acceleration for the next timestep. Atoms also have a "neighbor list", or a list of all the particles close enough to be considered in force calculations. (To prevent double counting, each neighbor list only contains the particles with index larger than the current particle's index.)
We create a global variable to store all of the particles, and create a function to report the current state of this variable.
// compute forces between all particles
void E(){
for(i j=N;j--;){
for(i k:P[j].W){
// compute distance between particles
v R=(P[j].r+-1*P[k].r).p();
d r=R.n();
// if squared distance less than 70 (approximately 6 * 3.4Å**2), the interaction will be non-negligible
if(r<70){
d O=r*r*r;
// this is the expression for the lennard–jones force.
// the second lennard–jones parameter, the depth of the potential well (120 kB), is factored in here.
R=2880*B*(2*S*S/O/O-S/O)/r*R;
// apply force to each particle
P[j].F(R);
P[k].F(-1*R);
}
}
}
}
Now, we create a function to calculate the forces between all pairs of particles. For each particle, we loop over the neighbor list and see if the distance is within six minima of the adjacent particle, using squared distance to avoid the expensive square-root calculation. If so, we calculate the force and apply it to each particle.
// build neighbor lists
void e(){
for(i j=0;j<N;j++){
// clear the old lists
P[j].W.clear();
for(i k=j+1;k<N;k++)
// if squared distance between particles less than 90 (e.g. close to above cutoff), add to neighbor list
if((P[j].r+-1*P[k].r).p().n()<90)
P[j].W.push_back(k);
}
}
Finally, we create a function to build the neighbor lists, which is a straightforward double loop. We add every particle which might conceivably be close enough to factor into the force calculations within 10–20 frames.
i main(){
i A=1e3; // initial temperature (1000 K)
i Y=5e3; // time to reach final temperature (5000 fs)
// initialize the system. each particle will be randomly placed until it isn't too close to other particles.
for(i j=N;j--;){
for(i a=999;a--;){
// generate random position
v r=v(rd(),rd(),rd());
i c=0;
// check for clashes with each extant particle
for(p X:P){
d D=(r+-1*X.r).p().n();
if(D<6.8)
c=1;
}
// if no clashes, add particle to list
if(!c){
P.push_back(p(r));
break;
}
}
}
To begin the program, we randomly initialize each particle and test if we're too close to any other particles. If not, we save the position and move on to the next particle. This is crude but works well enough for such a simple system.
// run MD! this is basically just the velocity verlet algorithm.
for(i t=0;t<=3e4;t++){
// update position
for(p& I:P)
I.r=(I.r+I.q+0.5*I.a).p();
// every 20 timesteps (20 fs), update neighbor lists
if(t%20==0)
e();
// compute forces
E();
// finish velocity verlet, and sum up the kinetic energy.
d K=0; // kinetic energy
for(p& I:P){
I.q=I.q+0.5*(I.b+I.a);
I.a=I.b;
I.b=v();
K+=20*I.q.n();
}
d T=2*K/(3*B*N); // temperature
// in the first 10 ps, apply berendsen thermostat to control temperature
if(t<2*Y){
d C=75; // target temperature
if(t<Y)
C=75*t/Y+(A-75)*(Y-t)/Y;
for(p& I:P)
I.q=sqrt(C/T)*I.q;
}
// every 100 fs after the first 10 ps, write the configuration to stdout
if(t%100==0&&t>2*Y)
w();
}
}
Finally, we simply have to run the simulation. We use the velocity Verlet algorithm with a 1 fs timestep, updating neighbor lists and writing to stdout periodically. The temperature is gradually lowered from 1000 K to 75 K over 5 ps, and temperature is controlled for the first 10 ps.
Hopefully this helps to shed some light on how simple a molecular dynamics simulation can be, and highlights the wonders of obfuscated C++!