moved stuff into function to remove a for loop

.clang-format

@@ -4,3 +4,4 @@ TabWidth: 4
 UseTab: Always
 BreakBeforeBraces: Attach
 ColumnLimit: 100
+AlignReturnType: Always

Makefile

@@ -6,7 +6,7 @@ ts := $(shell /usr/bin/date "+%d-%m__%H_%M_%S")
 .DEFAULT_GOAL = MD
 
 MD: $(SRC)/MD.cpp
-	$(CC) $(CFLAGS) $(SRC)MD.cpp -lm -march=native -mavx -O2 -ftree-vectorize -funroll-loops -pg -g -o ./out/MD
+	$(CC) $(CFLAGS) $(SRC)MD.cpp -lm -march=native -O2 -ftree-vectorize -funroll-loops -pg -g -o ./out/MD
 
 MDorig: $(SRC)/MD-original.cpp
 	$(CC) $(CFLAGS) $(SRC)MD-original.cpp -lm -O2 -pg -o ./out/MD-original

src/MD-original.cpp

@@ -83,8 +83,7 @@ double MeanSquaredVelocity ();
 //  Compute total kinetic energy from particle mass and velocities
 double Kinetic();
 
-int
+int main() {
-main () {
 
   //  variable delcarations
   int i;
@@ -94,11 +93,9 @@ main () {
   char trash[10000], prefix[1000], tfn[1000], ofn[1000], afn[1000];
   FILE *infp, *tfp, *ofp, *afp;
 
-	printf (
+  printf("\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
   printf("                  WELCOME TO WILLY P CHEM MD!\n");
-	printf (
+  printf("  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
   printf("\n  ENTER A TITLE FOR YOUR CALCULATION!\n");
   scanf("%s", prefix);
   strcpy(tfn, prefix);
@@ -108,11 +105,9 @@ main () {
   strcpy(afn, prefix);
   strcat(afn, "_average.txt");
 
-	printf (
+  printf("\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
   printf("                  TITLE ENTERED AS '%s'\n", prefix);
-	printf (
+  printf("  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
 
   /*     Table of values for Argon relating natural units to SI units:
    *     These are derived from Lennard-Jones parameters from the article
@@ -133,17 +128,14 @@ main () {
   //  Edit these factors to be computed in terms of basic properties in
   //  natural units of the gas being simulated
 
-	printf (
+  printf("\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+  printf("  WHICH NOBLE GAS WOULD YOU LIKE TO SIMULATE? (DEFAULT IS ARGON)\n");
-	printf (
-		"  WHICH NOBLE GAS WOULD YOU LIKE TO SIMULATE? (DEFAULT IS ARGON)\n");
   printf("\n  FOR HELIUM,  TYPE 'He' THEN PRESS 'return' TO CONTINUE\n");
   printf("  FOR NEON,    TYPE 'Ne' THEN PRESS 'return' TO CONTINUE\n");
   printf("  FOR ARGON,   TYPE 'Ar' THEN PRESS 'return' TO CONTINUE\n");
   printf("  FOR KRYPTON, TYPE 'Kr' THEN PRESS 'return' TO CONTINUE\n");
   printf("  FOR XENON,   TYPE 'Xe' THEN PRESS 'return' TO CONTINUE\n");
-	printf (
+  printf("  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
   scanf("%s", atype);
 
   if (strcmp(atype, "He") == 0) {
@@ -190,17 +182,13 @@ main () {
     timefac = 2.09618e-12;
     strcpy(atype, "Ar");
   }
-	printf (
+  printf("\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
   printf("\n                     YOU ARE SIMULATING %s GAS! \n", atype);
-	printf (
+  printf("\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
 
-	printf (
+  printf("\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"\n  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
   printf("\n  YOU WILL NOW ENTER A FEW SIMULATION PARAMETERS\n");
-	printf (
+  printf("  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
-		"  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
   printf("\n\n  ENTER THE INTIAL TEMPERATURE OF YOUR GAS IN KELVIN\n");
   scanf("%lf", &Tinit);
   // Make sure temperature is a positive number!
@@ -297,8 +285,7 @@ main () {
   Tavg = 0;
 
   int tenp = floor(NumTime / 10);
-	fprintf (
+  fprintf(ofp,
-		ofp,
           "  time (s)              T(t) (K)              P(t) (Pa)           "
           "Kinetic En. (n.u.)     Potential En. (n.u.) Total En. (n.u.)\n");
   printf("  PERCENTAGE OF CALCULATION COMPLETE:\n  [");
@@ -368,11 +355,9 @@ main () {
   Tavg /= NumTime;
   Z = Pavg * (Vol * VolFac) / (N * kBSI * Tavg);
   gc = NA * Pavg * (Vol * VolFac) / (N * Tavg);
-	fprintf (afp,
+  fprintf(afp, "  Total Time (s)      T (K)               P (Pa)      PV/nT "
-			 "  Total Time (s)      T (K)               P (Pa)      PV/nT "
                "(J/(mol K))         Z           V (m^3)              N\n");
-	fprintf (
+  fprintf(afp,
-		afp,
           " --------------   -----------        ---------------   "
           "--------------   ---------------   ------------   -----------\n");
   fprintf(afp,
@@ -407,8 +392,7 @@ main () {
   return 0;
 }
 
-void
+void initialize() {
-initialize () {
   int n, p, i, j, k;
   double pos;
 
@@ -452,8 +436,7 @@ initialize () {
 }
 
 //  Function to calculate the averaged velocity squared
-double
+double MeanSquaredVelocity() {
-MeanSquaredVelocity () {
 
   double vx2 = 0;
   double vy2 = 0;
@@ -473,8 +456,7 @@ MeanSquaredVelocity () {
 }
 
 //  Function to calculate the kinetic energy of the system
-double
+double Kinetic() { // Write Function here!
-Kinetic () { // Write Function here!
 
   double v2, kin;
 
@@ -494,8 +476,7 @@ Kinetic () { // Write Function here!
 }
 
 // Function to calculate the potential energy of the system
-double
+double Potential() {
-Potential () {
   double quot, r2, rnorm, term1, term2, Pot;
   int i, j, k;
 
@@ -524,8 +505,7 @@ Potential () {
 //   Uses the derivative of the Lennard-Jones potential to calculate
 //   the forces on each atom.  Then uses a = F/m to calculate the
 //   accelleration of each atom.
-void
+void computeAccelerations() {
-computeAccelerations () {
   int i, j, k;
   double f, rSqd;
   double rij[3]; // position of i relative to j
@@ -560,8 +540,7 @@ computeAccelerations () {
 }
 
 // returns sum of dv/dt*m/A (aka Pressure) from elastic collisions with walls
-double
+double VelocityVerlet(double dt, int iter, double *PE, FILE *fp) {
-VelocityVerlet (double dt, int iter, double *PE, FILE *fp) {
   int i, j, k;
 
   double psum = 0.;
@@ -593,13 +572,13 @@ VelocityVerlet (double dt, int iter, double *PE, FILE *fp) {
     for (j = 0; j < 3; j++) {
       if (r[i][j] < 0.) {
         v[i][j] *= -1.; //- elastic walls
-				psum += 2 * m * fabs (v[i][j])
+        psum += 2 * m * fabs(v[i][j]) /
-						/ dt; // contribution to pressure from "left" walls
+                dt; // contribution to pressure from "left" walls
       }
       if (r[i][j] >= L) {
         v[i][j] *= -1.; //- elastic walls
-				psum += 2 * m * fabs (v[i][j])
+        psum += 2 * m * fabs(v[i][j]) /
-						/ dt; // contribution to pressure from "right" walls
+                dt; // contribution to pressure from "right" walls
       }
     }
   }
@@ -617,8 +596,7 @@ VelocityVerlet (double dt, int iter, double *PE, FILE *fp) {
   return psum / (6 * L * L);
 }
 
-void
+void initializeVelocities() {
-initializeVelocities () {
 
   int i, j;
 
@@ -677,8 +655,7 @@ initializeVelocities () {
 }
 
 //  Numerical recipes Gaussian distribution number generator
-double
+double gaussdist() {
-gaussdist () {
   static bool available = false;
   static double gset;
   double fac, rsq, v1, v2;

src/MD.cpp

@@ -25,6 +25,7 @@
 
 */
 
+#include <emmintrin.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -475,22 +476,26 @@ double Kinetic() { // Write Function here!
   return kin;
 }
 
-// void
+void transposeMatrix(double mat[MAXPART][3], double matT[3][MAXPART]) {
-// transposeMatrix (double mat[MAXPART][3], double matT[3][MAXPART]) {
+  for (int i = 0; i < 3; i++) {
-//	for (int i = 0; i < 3; i++) {
+    for (int j = 0; j < MAXPART; j++) {
-//		for (int j = 0; j < MAXPART; j++) {
+      matT[i][j] = mat[j][i];
-//			matT[i][j] = mat[j][i];
+    }
-//		}
+  }
-//	}
+}
-// }
+void transposeMatrix2(double matT[MAXPART][3], double mat[3][MAXPART]) {
+  for (int i = 0; i < MAXPART; i++) {
+    for (int j = 0; j < 3; j++) {
+      matT[i][j] = mat[j][i];
+    }
+  }
+}
 
-double PotentialAndAcceleration() {
+double PotentialAndAcceleration(double dt) {
   memset(a, 0, sizeof(a));
   double Pot = 0.;
-  // double rT[3][MAXPART];
-  // transposeMatrix (r, rT);
 
-  for (int i = 0; i < N - 1; i++) {
+  for (int i = 0; i < N; i++) {
     for (int j = i + 1; j < N; j++) {
       double quot, rnorm, term1, term2;
       // CALCULATE POTENTIAL ENERGY
@@ -504,23 +509,24 @@ double PotentialAndAcceleration() {
         // ACCEL
         posItoJ[k] = distTmp;
       }
-
-      quot = sigma * sigma / dist;
-      term2 = quot * quot * quot;
-      Pot += epsilon_8 * term2 * (term2 - 1);
-
       //  From derivative of Lennard-Jones with sigma and epsilon
       //  set equal to 1 in natural units!
-      double distSqd = dist * dist;
+      double distSqd = dist * dist * dist;
-      double rSqd_inv4 = 1.0 / (distSqd * distSqd);
+      double rSqd_inv7 = distSqd * distSqd * dist;
-      double rSqd_inv7 = rSqd_inv4 / (distSqd * dist);
+      double f = (48. - (24. * distSqd)) / rSqd_inv7;
-      double f = 48.0 * rSqd_inv7 - 24 * rSqd_inv4;
       //  from F = ma, where m = 1 in natural units!
       for (int k = 0; k < 3; k++) {
         double tmp = posItoJ[k] * f;
         a[i][k] += tmp;
         a[j][k] -= tmp;
       }
+
+      quot = sigma * sigma / dist;
+      term2 = quot * quot * quot;
+      Pot += epsilon_8 * term2 * (term2 - 1.);
+    }
+    for (int j = 0; j < 3; j++) {
+      v[i][j] += 0.5 * a[i][j] * dt;
     }
   }
   return Pot;
@@ -603,21 +609,22 @@ double VelocityVerlet(double dt, int iter, double *PE, FILE *fp) {
   // printf("  Updated Positions!\n");
   for (i = 0; i < N; i++) {
     for (j = 0; j < 3; j++) {
-      r[i][j] += v[i][j] * dt + 0.5 * a[i][j] * dt * dt;
+      double tmp = 0.5 * a[i][j] * dt;
+      r[i][j] += v[i][j] * dt + tmp * dt;
 
-      v[i][j] += 0.5 * a[i][j] * dt;
+      v[i][j] += tmp;
     }
     // printf("  %i  %6.4e   %6.4e %6.4e\n",i,r[i][0],r[i][1],r[i][2]);
   }
   //  Update accellerations from updated positions
   // computeAccelerations ();
-  *PE = PotentialAndAcceleration();
+  *PE = PotentialAndAcceleration(dt);
   //  Update velocity with updated acceleration
-  for (i = 0; i < N; i++) {
+  // for (i = 0; i < N; i++) {
-    for (j = 0; j < 3; j++) {
+  //  for (j = 0; j < 3; j++) {
-      v[i][j] += 0.5 * a[i][j] * dt;
+  //    v[i][j] += 0.5 * a[i][j] * dt;
-    }
+  //  }
-  }
+  //}
 
   // Elastic walls
   for (i = 0; i < N; i++) {