Nils code after moving executable code out of the .h files.

Makefile

+.POSIX:
+
+CFLAGS = -O
+
+OBJS = md.o md-extra.o rvec.o
+INCS = md.h rvec.h
+
+md: $(OBJS) $(INCS)
+	$(CC) $(CFLAGS) rvec.o md.o md-extra.o -o md -lm
+
+clean:
+	rm -f md
+	rm -f *.o
+	rm -f *~

README.md

+# MD Uebung files
+
+This began life as Nils Uebung.
+In June 2022, I had to move executable code out of `rvec.h`. There was executable code which is not allowed in a `.h` file.
+
+I also removed some unnecessary lines from the Makefile.

md-extra.c

+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "rvec.h"
+#include "md.h"
+
+/* TODO: natom_per_dim should be dependent on box dimensions,
+   i.e. have a vector of natom_per_dim */
+void setup_fcc_lattice(rvec *x, uint natom, rvec box) {
+    uint i, j;
+    uint natom_per_dim = ceil(pow(natom, 1.0/NDIM));
+    rvec r;
+    uint c[NDIM];
+    for (i = 0; i < NDIM; i++)
+        c[i] = 0;
+    rvec_mul(&r, 1.0 / natom_per_dim, box);
+    for (i = 0; i < natom; i++) {
+        for (j = 0; j < NDIM; j++)
+            x[i][j] = c[j] * r[j];
+        c[0] += 1;
+        for (j = 0; j < NDIM - 1; j++) {
+            if (c[j] == natom_per_dim) {
+                c[j] = 0;
+                c[j+1] += 1;
+            }
+        }
+    }
+}
+
+void apply_pbc(uint natom, rvec box, rvec *x) {
+    for (uint i = 0; i < natom; i++) {
+        for (uint d = 0; d < NDIM; d++) {
+            while (x[i][d] < 0)
+                x[i][d] += box[d];
+            while (x[i][d] > box[d])
+                x[i][d] -= box[d];
+        }
+    }
+}
+
+void remove_com_velocity(uint natom, rvec *v, real *m) {
+    rvec vsum;
+    real msum = 0;
+    rvec_set(&vsum, 0);
+    for (uint i = 0; i < natom; i++) {
+        msum += m[i];
+        rvec_add(&vsum, vsum, v[i]);
+    }
+    for (uint i = 0; i < natom; i++)
+        rvec_muladd(&v[i], v[i], - m[i]/msum, vsum);
+}
+
+real com_velocity(uint natom, rvec *v, real *m) {
+    rvec vsum;
+    real msum = 0;
+    rvec_set(&vsum, 0);
+    for (uint i = 0; i < natom; i++)
+        msum += m[i];
+    for (uint i = 0; i < natom; i++)
+        rvec_muladd(&vsum, vsum, m[i]/msum, v[i]);
+    return rvec_len(vsum);
+}
+
+void rescale_temperature(uint natom, rvec *v, real *m, uint ndof, real kB, real Tnew) {
+    real Ekin = calc_kinetic_energy(natom, v, m); 
+    real Told = calc_temperature(ndof, Ekin, kB);
+    real vscale;
+    if (Told == 0.0) {
+        printf("error: couldn't rescale temperature, old temperature is zero (Told = 0)\n");
+        exit(EXIT_FAILURE);
+    }
+    vscale = sqrt(Tnew / Told);
+    for (uint i = 0; i < natom; i++)
+        rvec_mul(&v[i], vscale, v[i]);
+}
+
+void init_velocities(uint natom, rvec *v, real *m, real ndof, real kB, real Tnew) {
+    for (uint i = 0; i < natom; i++)
+        rvec_rand(&v[i]);
+    remove_com_velocity(natom, v, m);
+    rescale_temperature(natom, v, m, ndof, kB, Tnew);
+}

md.c

+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "rvec.h"
+#include "md.h"
+
+/* [ASSIGNMENT] calculate energy and forces of a WCA potential (a Lennard-Jones
+   potential truncated at the minimum) */
+real calc_energy_force(uint natom, rvec box, rvec *x, rvec *f) {
+    uint i, j;
+    real Epot, r, fscale;
+    rvec rij;
+    const real eps   = 2.0, sig = 1.5;
+    const real rcut  = pow(2.0, 1.0/6.0) * sig;
+    const real sig6  = pow(sig, 6);
+    const real sig12 = pow(sig, 12);
+    const real Eshift = eps;
+
+    rvecary_set(natom, f, 0);
+    Epot = 0;
+    for (i = 0; i < natom; i++) {
+        for (j = i + 1; j < natom; j++) {
+            /* interaction of atoms i and j */
+            /* Note: use rvec_sub_pbc(&rij, box, x[i], x[j]); to
+               calculate the distance between atoms i and j */
+        }
+    }
+    return Epot;
+}
+
+/* [ASSIGNMENT] calculate instantaneous temperature */
+real calc_temperature(uint ndof, real Ekin, real kB) {
+    return 1.0;
+}
+
+/* [ASSIGNMENT] calculate kinetic energy */
+real calc_kinetic_energy(uint natom, rvec *v, real *m) {
+    real Ekin = 0;
+    return Ekin;
+}
+
+/* [ASSIGNMENT] leapfrog integration step */
+real leapfrog_step(uint natom, rvec *x, rvec *v, rvec *f, real *m, real dt, rvec box) {
+    real Epot = 0;
+    /* x += dt * v */
+    /* calc energy, forces */
+    /* v += dt * f/m */
+    return Epot;
+}
+
+static void print_header() {
+    printf("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n",
+           "step", "time", "E/Na", "Epot/Na", "Ekin/Na", "T", "Edrift/Na", "rel.Edrift", "vcom");
+}
+
+static void print_step_info(uint step, uint natom, real dt,
+                     real E, real Epot, real Ekin, real Estart, real T, real vcom) {
+    printf("%10u\t%10.5f\t%10.5f\t%10.5f\t%10.5f\t%10.5f\t%10.5f\t%10.5f\t%10.5f\n",
+           step, step * dt, E/natom, Epot/natom, Ekin/natom, T,
+           (E - Estart)/natom, (E - Estart)/Estart, vcom);
+}
+
+int main(int argc, char **argv) {
+    uint  natom;          /* number of atoms */
+    ulong nstep,          /* number of time steps to run simulation */
+          nstep_print;    /* print simulation state every nstep_print steps */
+    rvec  *x,             /* atom positions, x[i] is vector of coordinates of atom i */
+          *v,             /* velocities */
+          *f;             /* forces */
+    real  *m;             /* atom masses */
+    rvec  box;            /* simulation box side lengths */
+    real  dt;             /* integrator timestep */
+
+    uint  ndof;           /* number of degrees of freedom */
+    real  volume,         /* simulation box volume */
+          density,        /* system density = natom / volume */
+          E,              /* total energy */
+          Epot,           /* potential energy */
+          Ekin,           /* kinetic energy */      
+          Estart,         /* total energy at start of simulation */
+          T,              /* instantaneous temperature */
+          Tstart;         /* starting temperature */
+    uint  step = 0;       /* integration step */
+    real  vcom;           /* center of mass velocity */
+
+    const real kB = 1.0;  /* Boltzmann's constant */
+    const real default_mass = 3.0;
+    const real max_rel_Edrift = 1.0;
+
+    /* parse command-line args */
+    if (argc != 7) {
+        printf("usage: %s natom density Tstart dt nstep nstep_print\n", argv[0]);
+        exit(EXIT_FAILURE);
+    }
+    natom       = atoi(argv[1]);
+    density     = atof(argv[2]);
+    Tstart      = atof(argv[3]);
+    dt          = atof(argv[4]);
+    nstep       = atol(argv[5]);
+    nstep_print = atol(argv[6]);
+
+    /* allocate memory */
+    x = calloc(natom, sizeof *x);
+    v = calloc(natom, sizeof *v);
+    f = calloc(natom, sizeof *f);
+    m = calloc(natom, sizeof *m);
+    if (x == NULL || v == NULL || f == NULL || m == NULL) {
+        printf("error: couldn't allocate memory\n");
+        exit(EXIT_FAILURE);
+    }
+
+    /* setup */
+    ndof   = NDIM * natom;
+    volume = natom / density;
+    rvec_set(&box, pow(volume, 1.0/NDIM));
+    setup_fcc_lattice(x, natom, box);
+    apply_pbc(natom, box, x);
+    for (uint i = 0; i < natom; i++)
+        m[i] = default_mass;
+    init_velocities(natom, v, m, ndof, kB, Tstart);
+
+    /* print info */
+    printf("simulation box: ");
+    rvec_print(box);
+    printf("[ x ]\n");
+    rvecary_print(x, natom);
+    printf("[ v ]\n");
+    rvecary_print(v, natom);
+    printf("[ f ]\n");
+    rvecary_print(f, natom);
+
+    /* calculate starting energy etc. */
+    Epot = calc_energy_force(natom, box, x, f);
+    Ekin = calc_kinetic_energy(natom, v, m); 
+    T = calc_temperature(ndof, Ekin, kB);
+    Estart = E = Epot + Ekin;
+    vcom = com_velocity(natom, v, m);
+
+    /* run simulation */
+    printf("START\n");
+    print_header();
+    print_step_info(step, natom, dt, E, Epot, Ekin, Estart, T, vcom);
+    for (step = 1; step <= nstep; step++) {
+        /* leapfrog integration step */
+        Epot = leapfrog_step(natom, x, v, f, m, dt, box);
+
+        /* apply periodic boundary conditions, move all atoms back
+           into box */
+        apply_pbc(natom, box, x);
+
+        /* calculate properties: kinetic energy, temperature */
+        Ekin = calc_kinetic_energy(natom, v, m);
+        T    = calc_temperature(ndof, Ekin, kB);
+        E    = Epot + Ekin;
+        vcom = com_velocity(natom, v, m);
+
+        /* periodically print out info */
+        if (step % nstep_print == 0) {
+            print_step_info(step, natom, dt, E, Epot, Ekin, Estart, T, vcom);
+        }
+
+        /* check if system has exploded */
+        if (fabs((E - Estart) / Estart) > max_rel_Edrift) {
+            printf("error: relative energy drift too large\n");
+            print_step_info(step, natom, dt, E, Epot, Ekin, Estart, T, vcom);
+            exit(EXIT_FAILURE);
+        }
+    }
+
+    /* clean up */
+    free(x);
+    free(v);
+    free(f);
+    free(m);
+
+    return EXIT_SUCCESS;
+}

md.h

+#ifndef MD_H
+#define MD_H
+
+void setup_fcc_lattice(rvec *x, uint natom, rvec box);
+void apply_pbc(uint natom, rvec box, rvec *x);
+void remove_com_velocity(uint natom, rvec *v, real *m);
+real com_velocity(uint natom, rvec *v, real *m);
+void rescale_temperature(uint natom, rvec *v, real *m, uint ndof, real kB, real Tnew);
+void init_velocities(uint natom, rvec *v, real *m, real ndof, real kB, real Tnew);
+
+/* these are the functions that have to be implemented */
+real calc_energy_force(uint natom, rvec box, rvec *x, rvec *f);
+real calc_temperature(uint ndof, real Ekin, real kB);
+real calc_kinetic_energy(uint natom, rvec *v, real *m);
+real leapfrog_step(uint natom, rvec *x, rvec *v, rvec *f, real *m, real dt, rvec box);
+
+
+#endif /* MD_H */

rvec.c

+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "rvec.h"
+
+/* u[i] = val for all i */
+void rvec_set(rvec *u, real val) {
+    for (uint i = 0; i < NDIM; i++)
+        (*u)[i] = val;
+}
+
+/* set the components of u to random values uniformly distributed in
+   the range [-0.5,0.5] */
+/* NOTE: direction of vector is not uniformly distributed */
+/* TODO: rand() is a bad RNG */
+void rvec_rand(rvec *u) {
+    for (uint i = 0; i < NDIM; i++)
+        (*u)[i] = ((real) rand() / RAND_MAX) - 0.5;
+}
+
+/* dot product of u and v */
+static real rvec_dot(rvec u, rvec v) {
+    real r = 0;
+    for (uint i = 0; i < NDIM; i++)
+        r += u[i] * v[i];
+    return r;
+}
+
+/* length of vector u */
+real rvec_len(rvec u) {
+    return sqrt(rvec_dot(u, u));
+}
+
+/* u = v + w */
+void rvec_add(rvec *u, rvec v, rvec w) {
+    for (uint i = 0; i < NDIM; i++)
+        (*u)[i] = v[i] + w[i];
+}
+
+/* u = v - w */
+static void rvec_sub(rvec *u, rvec v, rvec w) {
+    for (uint i = 0; i < NDIM; i++)
+        (*u)[i] = v[i] - w[i];
+}
+
+/* u = a * v */
+void rvec_mul(rvec *u, real a, rvec v) {
+    for (uint i = 0; i < NDIM; i++)
+        (*u)[i] = a * v[i];
+}
+
+/* u = v + r * w */
+void rvec_muladd(rvec *u, rvec v, real r, rvec w) {
+    for (uint i = 0; i < NDIM; i++)
+        (*u)[i] = v[i] + r * w[i];
+}
+
+/* calculate v - w, taking periodic boundary conditions into
+   account */
+static void rvec_sub_pbc(rvec *u, rvec box, rvec v, rvec w) {
+    for (uint i = 0; i < NDIM; i++) {
+        (*u)[i] = v[i] - w[i];
+        while ((*u)[i] < -box[i]/2)
+            (*u)[i] += box[i];
+        while ((*u)[i] > box[i]/2)
+            (*u)[i] -= box[i];
+    }
+}
+
+/* print vector _without_ trailing newline */
+static void rvec_puts(rvec u) {
+    for (uint i = 0; i < NDIM; i++)
+        printf("%10.5f ", u[i]);
+}
+
+/* print vector _with_ trailing newline */
+void rvec_print(rvec u) {
+    rvec_puts(u);
+    printf("\n");
+}
+
+
+/* functions for rvecary (an array of rvec) */
+
+void rvecary_set(uint n, rvec *a, real val) {
+    for (uint i = 0; i < n; i++)
+        rvec_set(&a[i], val);
+}
+
+/* a[i] = b[i] + r * c[i] */
+void rvecary_muladd(uint n, rvec *a, rvec *b, real r, rvec *c) {
+    for (uint i = 0; i < n; i++)
+        rvec_muladd(&a[i], b[i], r, c[i]);
+}
+
+/* print an array of rvec */
+void rvecary_print(rvec *a, uint n) {
+    for (uint i = 0; i < n; i++)
+        rvec_print(a[i]);
+}

rvec.h

+#ifndef RVEC_H
+#define RVEC_H
+
+#define NDIM 2
+typedef double real;
+typedef real rvec[NDIM];
+typedef unsigned int uint;
+typedef unsigned long ulong;
+
+void rvec_set(rvec *u, real val);
+void rvecary_set(uint n, rvec *a, real val);
+void rvec_print(rvec u);
+void rvecary_print(rvec *a, uint n);
+void rvec_mul(rvec *u, real a, rvec v);
+void rvec_add(rvec *u, rvec v, rvec w);
+void rvecary_muladd(uint n, rvec *a, rvec *b, real r, rvec *c);
+void rvec_muladd(rvec *u, rvec v, real r, rvec w);
+real rvec_len(rvec u);
+void rvec_rand(rvec *u);
+#endif /* RVEC_H */
+