#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// input parameters
int np = 500; // number of Ar particles
double density = 0.75; // density of system
double temp = 1.0; // desired temperature
double dt = 0.001; // time step
double numstep = 1000; // number of steps for MD
// constants and variables
double r[500][3]; // positions
double v[500][3]; // velocities
double f[500][3]; // forces
double L, t, temperature; // box length
double ke, pe, total; // kinetic energy, potential energy, total energy
int nc; // number of cells in each direction (cubic box)
// functions prototype
void initialize();
void calcforce();
void update();
int main ()
{
int p, i;
FILE *fh;
fh = fopen("sim-box.xyz", "w");
initialize();
// writing initial configuration to file
fprintf(fh, "%d\n\n", np);
for (i = 0; i < np; i++) {
fprintf(fh, "%s\t%f\t%f\t%f\n","Ar", r[i][0], r[i][1], r[i][2]);
}
calcforce();
t = 0;
for (t =0; t < numstep; t++) {
update();
// writing configuration at each step to file
for (i = 0; i < np; i++) {
fprintf(fh, "%s\t%f\t%f\t%f\n","Ar", r[i][0], r[i][1], r[i][2]);
}
ke = 0;
for (p = 0; p < np; p++) {
for (i = 0; i < 3; i++) {
ke += 0.5 * v[p][i] * v[p][i];
}
}
total = ke + pe;
printf("pe = %f ke = %f Total Energy = %f\n",pe, ke, total);
}
fclose(fh);
}
void initialize () {
int i, j, k, m, n, p, nc, IX, IY, IZ, IREF, M;
double L, CELL, CELL2, sumv2, sumv3, df, fs;
double Sum[3] = {0.0, 0.0, 0.0};
// compute length of box
L = powf( np / density, 1.0 / 3.0 );
nc = powf( np / 4, 1.0 / 3.0 );
/* CALCULATE THE SIDE OF THE UNIT CELL */
CELL = L / nc;
CELL2 = 0.5 * CELL;
/* BUILD THE UNIT CELL */
/* SUBLATTICE A */
r[0][0] = 0.0;
r[0][1] = 0.0;
r[0][2] = 0.0;
/* SUBLATTICE B */
r[1][0] = 0.0;
r[1][1] = CELL2;
r[1][2] = CELL2;
/* SUBLATTICE C */
r[2][0] = CELL2;
r[2][1] = 0.0;
r[2][2] = CELL2;
/* SUBLATTICE D */
r[3][0] = CELL2;
r[3][1] = CELL2;
r[3][2] = 0.0;
/* CONSTRUCT THE LATTICE FROM THE UNIT CELL */
M = 0;
for( IZ=0; IZ<5; IZ++ )
for( IY=0; IY<5; IY++ )
for( IX=0; IX<5; IX++ ){
for( IREF=0; IREF<4; IREF++ ){
r[IREF+M][0] = r[IREF][0] + CELL * IX ;
r[IREF+M][1] = r[IREF][1] + CELL * IY ;
r[IREF+M][2] = r[IREF][2] + CELL * IZ ;
}
M = M + 4;
}
// random distributed initial velocities
for (p = 0; p < np; p++)
for (i = 0; i < 3; i++) {
v[p][i] = (double)rand()/RAND_MAX - 0.5;
Sum[i] += v[p][i];
}
// with zero total momentum
for (p = 0; p < np; p++)
for (i = 0; i < 3; i++)
v[p][i] -= Sum[i] / np;
// rescaled to desired temperature
for (p = 0; p < np; p++)
for (i = 0; i < 3; i++)
sumv2 += v[p][i] * v[p][i];
df = 3 * np - 3;
fs = sqrt(temp * df / sumv2);
printf("temperature before scaling= %f\n", sumv2/df);
sumv3 = 0.0;
for (p = 0; p < np; p++)
for (i = 0; i < 3; i++) {
v[p][i] *= fs;
sumv3 += v[p][i] * v[p][i];
}
printf("temperature after scaling= %f\n", sumv3/df);
}
void calcforce () {
int i, j, k, sign;
double r2, rij[3], mag, L;
double r2i, r6i, r12i, ff, fij[3];
for (i = 0; i < np; i++)
for (k = 0; k < 3; k++)
f[i][k] = 0.0;
pe = 0.0;
L = powf(np / density, 1.0 / 3.0 );
// loop over all pairs of particles
for (i = 0; i < np - 1; i++) {
for (j = i + 1; j < np; j++) {
r2 = 0.0;
for (k = 0; k < 3; k++) {
rij[k] = r[i][k] - r[j][k];
sign = rij[k] > 0.0 ? 1 : -1;
mag = rij[k] * sign;
if (mag > L / 2) // use nearest image
rij[k] -= L * sign;
r2 += rij[k] * rij[k];
}
r2i = 1 / r2;
r6i = r2i * r2i * r2i;
r12i = r6i * r6i;
ff = r2i * (48 * r12i - 24 * r6i);
for (k = 0; k < 3; k++) {
fij[k] = rij[k] * ff;
f[i][k] += fij[k];
f[j][k] -= fij[k];
}
pe += 4 * (r12i - r6i);
}
}
}
void update () {
int p, k;
// compute length of box
L = powf( np / density, 1.0 / 3.0 );
for (p = 0; p < np; p++)
for (k = 0; k < 3; k++) {
r[p][k] += v[p][k] * dt + 0.5 * f[p][k] * dt * dt;
// impose periodic boundary conditions
if (r[p][k] < 0.0)
r[p][k] += L;
if (r[p][k] >= L)
r[p][k] -= L;
v[p][k] += 0.5 * f[p][k] * dt;
}
calcforce();
for (p = 0; p < np; p++)
for (k = 0; k < 3; k++)
v[p][k] += 0.5 * f[p][k] * dt;
}
