2.7 local / 7 search

Merge Potential and Acceleration functions. Still need to optimize more

src/MD-original.cpp

@@ -561,7 +561,7 @@ computeAccelerations () {
 
 // returns sum of dv/dt*m/A (aka Pressure) from elastic collisions with walls
 double
-VelocityVerlet (double dt, int iter, FILE *fp) {
+VelocityVerlet (double dt, int iter, double *PE, FILE *fp) {
 	int i, j, k;
 
 	double psum = 0.;

src/MD.cpp

@@ -26,7 +26,6 @@
 */
 
 #include <math.h>
-#include <omp.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -69,7 +68,7 @@ void initialize ();
 //  update positions and velocities using Velocity Verlet algorithm
 //  print particle coordinates to file for rendering via VMD or other animation
 //  software return 'instantaneous pressure'
-double VelocityVerlet (double dt, int iter, FILE *fp);
+double VelocityVerlet (double dt, int iter, double *PE, FILE *fp);
 //  Compute Force using F = -dV/dr
 //  solve F = ma for use in Velocity Verlet
 void computeAccelerations ();
@@ -335,7 +334,7 @@ main () {
 		// Also computes the Pressure as the sum of momentum changes from
 		// wall collisions / timestep which is a Kinetic Theory of gasses
 		// concept of Pressure
-		Press = VelocityVerlet (dt, i + 1, tfp);
+		Press = VelocityVerlet (dt, i + 1, &PE, tfp);
 		Press *= PressFac;
 
 		//  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -346,7 +345,6 @@ main () {
 		//  extract the gas constant
 		mvs = MeanSquaredVelocity ();
 		KE = Kinetic ();
-		PE = Potential ();
 
 		// Temperature from Kinetic Theory
 		Temp = m * mvs / (3 * kB) * TempFac;
@@ -497,40 +495,83 @@ Kinetic () { // Write Function here!
 	// printf("  Total Kinetic Energy is %f\n",N*mvs*m/2.);
 	return kin;
 }
+double
+PotentialAndAcceleration () {
+	double quot, rnorm, term1, term2, Pot;
+	int i, j, k;
+	double f, dist, tmp = 0.;
+	double distTmp = 0;
+	double posItoJ[3]; // position of i relative to j
+	double epsilon_8 = 8. * epsilon;
+
+	memset (a, 0, sizeof (a));
+
+	Pot = 0.;
+
+	for (i = 0; i < N - 1; i++) {
+		for (j = i + 1; j < N; j++) {
+			// CALCULATE POTENTIAL ENERGY
+			dist = 0.;
+			distTmp = 0;
+
+			for (int k = 0; k < 3; k++) {
+				// POT
+				distTmp = r[i][k] - r[j][k];
+				dist += distTmp * distTmp;
+				// ACCEL
+				posItoJ[k] = distTmp;
+			}
+
+			quot = sigma / sqrt (dist);
+			term2 = quot * quot;
+			term2 = term2 * term2 * term2;
+			term1 = term2 * term2;
+
+			Pot += 8 * epsilon * (term1 - term2);
+
+			//  From derivative of Lennard-Jones with sigma and epsilon
+			//  set equal to 1 in natural units!
+			double rSqd_inv7
+				= 1.0 / (dist * dist * dist * dist * dist * dist * dist);
+			double rSqd_inv4 = rSqd_inv7 * dist * dist * dist;
+			f = 24 * (2 * rSqd_inv7 - rSqd_inv4);
+
+			//  from F = ma, where m = 1 in natural units!
+			for (k = 0; k < 3; k++) {
+				tmp = posItoJ[k] * f;
+				a[i][k] += tmp;
+				a[j][k] -= tmp;
+			}
+		}
+	}
+	return Pot;
+}
 
 // Function to calculate the potential energy of the system
 double
 Potential () {
 	double quot, r2, rnorm, term1, term2, Pot;
-	int i, j, k;
+	int i, j;
 
 	Pot = 0.;
 
-#pragma omp parallel for reduction(+ : Pot) private(j, k, r2, rnorm, quot,    \
-														term1, term2)
-	for (i = 0; i < N; i++) {
-		for (j = 0; j < N; j++) {
+	for (i = 0; i < N - 1; i++) {
+		for (j = i + 1; j < N; j++) {
+			r2 = 0.;
 
-			if (j != i) {
-				r2 = 0.;
-				for (k = 0; k < 3; k++) {
-					double tmp = r[i][k] - r[j][k];
-					r2 += tmp * tmp;
-				}
-				rnorm = sqrt (r2);
-				quot = sigma / rnorm;
-				term2 = quot * quot;
-				term2 = term2 * term2 * term2;
-				term1 = term2 * term2;
-
-				double local_energy = 4 * epsilon * (term1 - term2);
-
-#pragma omp atomic
-				Pot += local_energy;
+			for (int k = 0; k < 3; k++) {
+				double tmp = r[i][k] - r[j][k];
+				r2 += tmp * tmp;
 			}
+
+			quot = sigma / sqrt (r2);
+			term2 = quot * quot;
+			term2 = term2 * term2 * term2;
+			term1 = term2 * term2;
+
+			Pot += 8 * epsilon * (term1 - term2);
 		}
 	}
-
 	return Pot;
 }
 
@@ -540,34 +581,34 @@ Potential () {
 void
 computeAccelerations () {
 	int i, j, k;
-	double f, rSqd;
+	double f, rSqd, tmp = 0.;
 	double rij[3]; // position of i relative to j
 
-	for (i = 0; i < N; i++) { // set all accelerations to zero
-		for (k = 0; k < 3; k++) {
-			a[i][k] = 0;
-		}
-	}
+	memset (a, 0, sizeof (a));
 	for (i = 0; i < N - 1; i++) { // loop over all distinct pairs i,j
 		for (j = i + 1; j < N; j++) {
 			// initialize r^2 to zero
 			rSqd = 0;
 
+			//  component-by-componenent position of i relative
+			//  to j
 			for (k = 0; k < 3; k++) {
-				//  component-by-componenent position of i relative
-				//  to j
 				rij[k] = r[i][k] - r[j][k];
-				//  sum of squares of the components
-				rSqd += rij[k] * rij[k];
+				tmp = rij[k];
+				rSqd += tmp * tmp;
 			}
-
 			//  From derivative of Lennard-Jones with sigma and epsilon
 			//  set equal to 1 in natural units!
-			f = 24 * (2 * pow (rSqd, -7) - pow (rSqd, -4));
+			double rSqd_inv7
+				= 1.0 / (rSqd * rSqd * rSqd * rSqd * rSqd * rSqd * rSqd);
+			double rSqd_inv4 = rSqd_inv7 * rSqd * rSqd * rSqd;
+			f = 24 * (2 * rSqd_inv7 - rSqd_inv4);
+
+			//  from F = ma, where m = 1 in natural units!
 			for (k = 0; k < 3; k++) {
-				//  from F = ma, where m = 1 in natural units!
-				a[i][k] += rij[k] * f;
-				a[j][k] -= rij[k] * f;
+				tmp = rij[k] * f;
+				a[i][k] += tmp;
+				a[j][k] -= tmp;
 			}
 		}
 	}
@@ -575,7 +616,7 @@ computeAccelerations () {
 
 // returns sum of dv/dt*m/A (aka Pressure) from elastic collisions with walls
 double
-VelocityVerlet (double dt, int iter, FILE *fp) {
+VelocityVerlet (double dt, int iter, double *PE, FILE *fp) {
 	int i, j, k;
 
 	double psum = 0.;
@@ -594,7 +635,8 @@ VelocityVerlet (double dt, int iter, FILE *fp) {
 		// 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 ();
+	// computeAccelerations ();
+	*PE = PotentialAndAcceleration ();
 	//  Update velocity with updated acceleration
 	for (i = 0; i < N; i++) {
 		for (j = 0; j < 3; j++) {