more testing

C/molecular_dynamics/include/dbg.h

@@ -15,9 +15,9 @@
 
 #define CLEAN_ERRNO() (errno == 0 ? "None" : strerror(errno))
 
-#define LOG_ERROR(M, ...) fprintf(stderr, "[ERROR] (%s:%d: errno: %s) " M "\n", __FILE__, __LINE__, clean_errno(), ##__VA_ARGS__)
+#define LOG_ERROR(M, ...) fprintf(stderr, "[ERROR] (%s:%d: errno: %s) " M "\n", __FILE__, __LINE__, CLEAN_ERRNO(), ##__VA_ARGS__)
 
-#define LOG_WARN(M, ...) fprintf(stderr, "[WARN] (%s:%d: errno: %s) " M "\n", __FILE__, __LINE__, clean_errno(), ##__VA_ARGS__)
+#define LOG_WARN(M, ...) fprintf(stderr, "[WARN] (%s:%d: errno: %s) " M "\n", __FILE__, __LINE__, CLEAN_ERRNO(), ##__VA_ARGS__)
 
 #define LOG_INFO(M, ...) fprintf(stderr, "[INFO] (%s:%d) " M "\n", __FILE__, __LINE__, ##__VA_ARGS__)
 

C/molecular_dynamics/runner.c

@@ -8,7 +8,7 @@
 
 int main(void) {
   LennardJonesPotential *ljp = new_LennardJonesPotential(3.40, 0.997);
-  MD_BoxParameters *mdb = new_MD_BoxParameters(20, 20, 20, 2.0 * ljp->sigma, 2.0);
+  MD_BoxParameters *mdb = new_MD_BoxParameters(5, 5, 5, 2.0 * ljp->sigma, 2.0);
   MD_RunParameters *mdp = new_MD_RunParameters(0.0, 0.01, 1200, 200, 500);
   MD_Report *report = new_Report();
   double box_length;
@@ -29,12 +29,10 @@ int main(void) {
     verlet_IteratePosition(particle, collection_size, mdp->dt, box_length, apply_Periodic);
     verlet_IterateVelocity(particle, collection_size, mdp->dt);
 
-    calculate_Forces(ljp, particle, collection_size, box_length);
+    MD_SystemEnergy *energy = calculate_Forces(ljp, particle, collection_size, box_length);
 
     verlet_IterateVelocity(particle, collection_size, mdp->dt);
 
-    MD_SystemEnergy *energy = calculate_SystemEnergy(ljp, particle, collection_size, box_length);
-
     mdp->t += mdp->dt;
 
     if (i % report_interval == 0 && i > mdp->equilibrium) {
@@ -48,6 +46,7 @@ int main(void) {
     destroy_SystemEnergy(energy);
   }
 
+  // Final positions
   for (int i = 0; i < collection_size; i++) {
     DEBUG_PRINT("%9.6f %9.6f %9.6f ", particle[i]->x, particle[i]->y, particle[i]->z);
   }

C/molecular_dynamics/src/md.c

@@ -5,6 +5,7 @@
 #include "md_math.h"
 #include "md_separation.h"
 #include "md_accel.h"
+#include "md_systemenergy.h"
 #include "lennard_jones.h"
 #include "dbg.h"
 
@@ -62,17 +63,24 @@ void apply_Periodic(Particle *particle, double length) {
   }
 }
 
-void calculate_Forces(LennardJonesPotential *ljp, ParticleCollection *particle, int collection_size, double length) {
+MD_SystemEnergy * calculate_Forces(LennardJonesPotential *ljp, ParticleCollection *particle, int collection_size, double length) {
   MD_Separation *sep;
   MD_Accel *accel;
+  MD_SystemEnergy *system = new_SystemEnergy();
+  double r, ke, pe;
+  int m;
 
   for (int i = 0; i < collection_size; i++) {
     particle[i]->ax = 0.0;
     particle[i]->ay = 0.0;
     particle[i]->az = 0.0;
   }
 
-  for (int m = 0; m < collection_size - 1; m++) {
+  ke = pe = 0.0;
+
+  for (m = 0; m < collection_size - 1; m++) {
+    ke += calculate_KineticEnergy(particle[m]);
+
     for (int n = m + 1; n < collection_size; n++) {
       sep = new_MD_Separation(particle[m], particle[n], length);
       accel = new_MD_Accel(ljp, sep); // calculates (dx/r)*force
@@ -85,16 +93,27 @@ void calculate_Forces(LennardJonesPotential *ljp, ParticleCollection *particle,
       particle[n]->ay -= accel->ay;
       particle[n]->az -= accel->az;
 
+      r = calculate_SeparationMagnitude(sep);
+      pe += LennardJones_PotentialEnergy(ljp, r);
+
       destroy_MD_Separation(sep);
       destroy_MD_Accel(accel);
     }
   }
+
+  // Account for the energy of the last particle
+  ke += calculate_KineticEnergy(particle[collection_size]);
+
+  system->pe = pe;
+  system->ke = ke;
+  system->te = pe + ke;
+
+  return system;
 }
 
 double initialize_Collection(ParticleCollection *p, MD_BoxParameters *mdb) {
-  double x, y, z;
-  double length;
   int collection_size = mdb->Nx * mdb->Ny * mdb->Nz;
+  double x, y, z;
 
   x = y = z = mdb->spacing;
 
@@ -152,9 +171,7 @@ double initialize_Collection(ParticleCollection *p, MD_BoxParameters *mdb) {
     DEBUG_PRINT("%12.6f %12.6f %12.6f", p[i]->vx, p[i]->vy, p[i]->vz);
   }
 
-  length = x;
-
-  DEBUG_PRINT("box length: %12.6f", length);
+  DEBUG_PRINT("box length: %12.6f", x);
 
-  return length;
+  return x;
 }

C/molecular_dynamics/src/md.h

@@ -3,6 +3,7 @@
 
 #include "particle.h"
 #include "lennard_jones.h"
+#include "md_systemenergy.h"
 
 typedef struct MD_BoxParameters {
   int Nx;
@@ -28,6 +29,6 @@ void destroy_MD_BoxParameters(MD_BoxParameters *);
 
 void apply_Periodic(Particle *, double);
 double initialize_Collection(ParticleCollection *, MD_BoxParameters *);
-void calculate_Forces(LennardJonesPotential *, ParticleCollection *, int, double);
+MD_SystemEnergy * calculate_Forces(LennardJonesPotential *, ParticleCollection *, int, double);
 
 #endif

C/molecular_dynamics/src/md_systemenergy.c

@@ -20,35 +20,15 @@ void destroy_SystemEnergy(MD_SystemEnergy *energy) {
   free(energy);
 }
 
-MD_SystemEnergy * calculate_SystemEnergy(LennardJonesPotential *ljp, ParticleCollection *particle, int collection_size, double length) {
-  double pe, ke, r;
-  MD_Separation *sep;
-  MD_SystemEnergy *system = new_SystemEnergy();
-
-  pe = ke = 0.0;
-
-  for (int i = 0; i < collection_size; i++) {
-    ke += 
-        particle[i]->vx * particle[i]->vx + 
-        particle[i]->vy * particle[i]->vy +
-        particle[i]->vz * particle[i]->vz;
-    if (i < collection_size - 1) {
-      for (int j = i + 1; j < collection_size; j++) {
-        sep = new_MD_Separation(particle[i], particle[j], length);
-
-        r = calculate_SeparationMagnitude(sep);
-        pe += LennardJones_PotentialEnergy(ljp, r);
-
-        destroy_MD_Separation(sep);
-      }
-    }
-  }
+double calculate_KineticEnergy(Particle *p) {
+  double ke = 0.0;
 
-  ke *= 0.5;
+  ke += 
+    p->vx * p->vx + 
+    p->vy * p->vy +
+    p->vz * p->vz;
 
-  system->pe = pe;
-  system->ke = ke;
-  system->te = pe + ke;
+  ke *= 0.5;
 
-  return system;
+  return ke;
 }

C/molecular_dynamics/src/md_systemenergy.h

@@ -12,6 +12,6 @@ typedef struct MD_SystemEnergy {
 
 MD_SystemEnergy * new_SystemEnergy();
 void destroy_SystemEnergy(MD_SystemEnergy *);
-MD_SystemEnergy * calculate_SystemEnergy(LennardJonesPotential *, ParticleCollection *, int, double);
+double calculate_KineticEnergy(Particle *p);
 
 #endif

C/molecular_dynamics/test/minunit.h

@@ -8,13 +8,13 @@
 
 #define mu_suite_start() char *message = NULL
 
-#define mu_assert(test, message) if (!(test)) { log_err(message); return message; }
-#define mu_run_test(test) debug("\n-----%s", " " #test); \
+#define mu_assert(test, message) if (!(test)) { LOG_ERROR(message); return message; }
+#define mu_run_test(test) DEBUG_PRINT("\n-----%s", " " #test); \
     message = test(); tests_run++; if (message) return message;
 
 #define RUN_TESTS(name) int main(int argc, char *argv[]) {\
     argc = 1; \
-    debug("----- RUNNING: %s", argv[0]);\
+    DEBUG_PRINT("----- RUNNING: %s", argv[0]);\
         printf("----\nRUNNING: %s\n", argv[0]);\
         char *result = name();\
         if (result != 0) {\
@@ -27,7 +27,6 @@
         exit(result != 0);\
 }
 
-
 int tests_run;
 
 #endif