/*
* serdy.c
*
* Example of a SERial molecular DYnamics simulation in C
*
* Input (from standard input or filename given on command line):
* N = <number of particles>
* rho = <density of the system>
* T = <initial temperature (standard deviation of the velocities)>
* runtime = <runtime>
* dt = <time step>
* seed = <random number generator seed>
* equil = <equilibration time>
* usecells
* (when the latter flag is not present, cells will not be used
* within processes)
*
* Output (from standard output): lines containing:
* time, energy E, pot. en. U, kin. en. K, temperature T, fluctuations, walltime-per-step, walltime-overall
*
* Notes:
* - To compile: make
* - To run: ./pardy input.ini
* where input parameters are listed in the file "input.ini", and
* output is sent to stdout.
* - Appropriate values for the equilibrium time are best found by
* doing a short run and seeing when the potential energy has reach
* a stationary value.
* - All reported energies values are divided by the number of particles N.
* - Fluctuations are the root mean square of E-<E> with <E> the mean energy E.
*
* Ramses van Zon, 13 November 2008
* - In Nov 25 2009:
* + condensed code
* + OpenMP version
* - In Aug 2016:
* + Renamed to ljhpc, included mpi header (not used yet)
* - In Sep 2016:
* + Added cells to determine interacting pairs.
* + Added MPI
* + Renamed to pardy
* + Change input format and modularized
*/
#include <assert.h>
#include <string.h>
#include <stdbool.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ndmalloc.h"
#include "lcg.h"
#include "lattice.h"
#include "debug.h"
#include "system.h"
#include "atom.h"
#include "inifile.h"
#include "global.h"
#include "estimates.h"
#include "forces.h"
#include "cells.h"
#include "elapsed.h"
void sortParticlesAndComputeForces(int *N, atom_t atoms[], interaction_pairs_t* p, system_t* sys)
{
/* sort particles into cells */
cellDivide(atoms, sys);
/* find neighbors */
findPairsFromCells(p, sys);
/* compute forces */
sys->U = computeForces(*N, atoms, p, sys->L, false);
}
double periodic(double u, double L)
{
while (u < -0.5*L)
u += L;
while (u >= 0.5*L)
u -= L;
return u;
}
void exchangeParticles(atom_t atoms[], system_t* sys)
{
for (int i =0; i < sys->N; ++i) {
atoms[i].rx = periodic(atoms[i].rx, sys->L);
atoms[i].ry = periodic(atoms[i].ry, sys->L);
atoms[i].rz = periodic(atoms[i].rz, sys->L);
}
}
void initialize(atom_t atoms[], interaction_pairs_t* p, system_t* sys)
{
/* generate positions */
/* initialize lattice position generator, then pick at most Ntot (will be less if parallel) */
lattice_t lat = init_partial_lattice(sys->L, sys->Ntot,
(double[3]){-0.5*sys->L,-0.5*sys->L,-0.5*sys->L},
(double[3]){sys->L,sys->L,sys->L});
int i;
for (i = 0; i < sys->Ntot; ++i) {
long long index = make_lattice_position(&lat, &(atoms[i].rx), &(atoms[i].ry), &(atoms[i].rz));
if (index != NO_LATTICE_POSITION) {
atoms[i].index = index;
} else {
break; /* no more points */
}
}
sys->N = i;
/* generate momenta */
lcg_t rng = lcg_init_rand48(sys->seed); /* initialize random number
generator used in gaussian
instead of "srand48(seed)" */
double scale = sqrt(sys->T0);
sys->K = 0;
/* must do some skipping (4 random numbers per particle)*/
lcg_skip(&rng, 4*atoms[0].index);
for (int i = 0; i < sys->N; ++i){
atoms[i].px = scale*lcg_normal(&rng);
atoms[i].py = scale*lcg_normal(&rng);
atoms[i].pz = scale*lcg_normal(&rng);
lcg_normal(&rng); /* must have even number of normal random numbers */
sys->K += atoms[i].px*atoms[i].px + atoms[i].py*atoms[i].py + atoms[i].pz*atoms[i].pz;
if (i<(sys->N)-1) {
int index_diff = atoms[i+1].index - atoms[i].index;
if (index_diff > 1)
lcg_skip(&rng, 4*(index_diff-1));
}
}
/* compute instantaneous kinetic temperature and energy */
sys->K *= 0.5;
sortParticlesAndComputeForces(&(sys->N), atoms, p, sys);
double Utot = sys->U, Ktot = sys->K;
/* report results */
printf("# step time E U K T <[E-<E>]^2> walltime(s) cum.walltime(s)\n");
printf("%7d %7.3f %7.3f %7.3f %7.3f %7.3f ..... 0.000 0.000\n",
0, 0., (Ktot+Utot) / (sys->Ntot), Utot / (sys->Ntot), Ktot / (sys->Ntot), (2*Ktot)/(DIM*(sys->Ntot)));
}
/* Verlet integration step */
void integrateStep(system_t* sys, atom_t atoms[], interaction_pairs_t* p)
{
/* first momentum half step (forces assume to be computed in
initialization or in the last step) */
for (int i = 0; i < sys->N; ++i) {
atoms[i].px += 0.5*sys->dt*atoms[i].fx;
atoms[i].py += 0.5*sys->dt*atoms[i].fy;
atoms[i].pz += 0.5*sys->dt*atoms[i].fz;
}
/* full free motion step */
for (int i = 0; i<sys->N; ++i) {
atoms[i].rx = atoms[i].rx + sys->dt*atoms[i].px;
atoms[i].ry = atoms[i].ry + sys->dt*atoms[i].py;
atoms[i].rz = atoms[i].rz + sys->dt*atoms[i].pz;
}
/* send atoms that have digressed too far to other side (pbc)*/
exchangeParticles(atoms, sys);
/* positions were changed, so recompute the forces */
sortParticlesAndComputeForces(&(sys->N), atoms, p, sys);
/* final momentum half-step */
double K2 = 0;
for (int i = 0; i < sys->N; ++i) {
atoms[i].px+=0.5*sys->dt*atoms[i].fx;
atoms[i].py+=0.5*sys->dt*atoms[i].fy;
atoms[i].pz+=0.5*sys->dt*atoms[i].fz;
K2 += atoms[i].px*atoms[i].px + atoms[i].py*atoms[i].py + atoms[i].pz*atoms[i].pz;
}
/* finish computing K */
sys->K = 0.5*K2;
}
/* integration and measurement */
void run(system_t* sys)
{
/* setup cell stuff */
for (int d = 0; d < DIM; ++d) {
if (sys->usecells)
sys->nc[d] = (int)(sys->L/rc);
else
sys->nc[d] = 1;
sys->totnc[d] = sys->nc[d];
sys->minc[d] = 0;
sys->maxc[d] = sys->nc[d] - 1;
sys->cellsize[d] = sys->L/sys->nc[d];
}
sys->ncprod = sys->nc[0] * sys->nc[1] * sys->nc[2];
sys->start_of_cell = ndmalloc(sizeof(sys->start_of_cell[0]), 1, sys->ncprod);
sys->n_in_cell = ndmalloc(sizeof(sys->n_in_cell[0]), 1, sys->ncprod);
/* allocate buffers */
const int Nmax = estimateNmax(sys->N, sys->rho, (double[3]){sys->L, sys->L, sys->L}, sys->cellsize);
const long long Npairsmax = estimateNpairsmax(sys->N, sys->rho, (double[3]){sys->L, sys->L, sys->L}, sys->cellsize);
interaction_pairs_t p;
atom_t* atoms = ndmalloc(sizeof(atom_t), 1, Nmax);
int numSteps = (int)(sys->runtime/sys->dt+0.5);
int equilSteps = (int)(sys->equil/sys->dt+0.5);
int count; /* counts time steps */
int numPoints = 0; /* counts measurements */
double sumE = 0; /* total energy accumulated over steps */
double sumE2 = 0; /* total energy squared accumulated */
double avgE, avgE2, fluctE; /* average energy, its square, and its fluctuations */
interaction_pairs_alloc(&p, Npairsmax);
/* draw initial conditions */
initialize(atoms, &p, sys);
double walltimestart = elapsed_etime();
double walltimelast = walltimestart;
/* perform equilibration steps */
for (count = 0; count < equilSteps; ++count) {
integrateStep(sys, atoms, &p);
double t = elapsed_etime();
printf("%7d %7.3f %7.3f %7.3f %7.3f %7.3f ..... %10.3lf %10.3lf\n",
count+1, (count+1)*sys->dt, (sys->K+sys->U)/sys->Ntot, sys->U/sys->Ntot, sys->K/sys->Ntot, (2*sys->K)/(DIM*sys->Ntot), t-walltimelast, t-walltimestart );
walltimelast = t;
}
/* perform rest of time steps */
for (count = equilSteps; count<numSteps; ++count) {
integrateStep(sys, atoms, &p);
double Etot = sys->K + sys->U;
sumE += Etot;
sumE2 += Etot * Etot;
++numPoints;
avgE = sumE/numPoints;
avgE2 = sumE2/numPoints;
if (avgE2 > avgE*avgE) {
fluctE = sqrt(avgE2-avgE*avgE);
} else { /* correct for round-off */
fluctE = 0.0;
}
double t = elapsed_etime();
printf("%7d %7.3f %7.3f %7.3f %7.3f %7.3f %10.2e %10.3lf %10.3lf\n",
count+1, (count+1)*sys->dt, Etot/sys->Ntot, sys->U/sys->Ntot, sys->K/sys->Ntot, (2*sys->K)/(DIM*sys->Ntot), fluctE/sys->Ntot, t-walltimelast, t-walltimestart);
walltimelast = t;
}
interaction_pairs_free(&p);
ndfree(atoms);
ndfree(sys->n_in_cell);
ndfree(sys->start_of_cell);
}
/* program start */
int main(int argc, char* argv[])
{
system_t sys;
FILE* file = (argc>1)?fopen(argv[1],"r"):stdin;
if (file == NULL)
exit(1);
key_value_table_t* ini = calloc(1,sizeof(*ini));
kvt_read(ini, file);
sys.Ntot = kvt_lookup_long_long(ini, "N");
sys.rho = kvt_lookup_double(ini, "rho");
sys.T0 = kvt_lookup_double(ini, "T");
sys.runtime = kvt_lookup_double(ini, "runtime");
sys.dt = kvt_lookup_double(ini, "dt");
sys.seed = kvt_lookup_long(ini, "seed");
sys.equil = kvt_lookup_double(ini, "equil");
sys.usecells = kvt_lookup_entry(ini, "usecells") != NULL;
sys.L = cbrt(sys.Ntot/sys.rho);
if (argc > 1 && file != NULL)
fclose(file);
free(ini);
printf("#Ntot=%lld rho=%f L=%f T=%f runtime=%f dt=%f seed=%ld equil=%f usecells=%d\n",
sys.Ntot,sys.rho,sys.L,sys.T0,sys.runtime,sys.dt,sys.seed,sys.equil,sys.usecells);
if (sanity_check(&sys, rc, true))
run(&sys);
else
exit(1);
return 0;
}
