diff --git a/Makefile b/Makefile
index 2fc2a85..d202571 100644
--- a/Makefile
+++ b/Makefile
@@ -1,11 +1,11 @@
CC = gcc
SRC = src/
-CFLAGS = # none
+CFLAGS = -march=native -mtune=native -mavx -O2 -ftree-vectorize -fopenmp
.DEFAULT_GOAL = MD.exe
-MD.exe: $(SRC)/MD.cpp
- $(CC) $(CFLAGS) $(SRC)MD.cpp -lm -march=native -mtune=native -mavx -O2 -ftree-vectorize -funroll-loops -o MD.exe
+MD.exe: $(SRC)/MD2.cpp
+ $(CC) $(CFLAGS) $(SRC)MD2.cpp -lm -o MD.exe
clean:
rm ./MD.exe
diff --git a/relatorio.pdf b/relatorio.pdf
new file mode 100644
index 0000000..50e2519
Binary files /dev/null and b/relatorio.pdf differ
diff --git a/src/MD2.cpp b/src/MD2.cpp
new file mode 100644
index 0000000..ffb6c61
--- /dev/null
+++ b/src/MD2.cpp
@@ -0,0 +1,837 @@
+/*
+ MD.c - a simple molecular dynamics program for simulating real gas properties
+ of Lennard-Jones particles.
+
+ Copyright (C) 2016 Jonathan J. Foley IV, Chelsea Sweet, Oyewumi Akinfenwa
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see .
+
+ Electronic Contact: foleyj10@wpunj.edu
+ Mail Contact: Prof. Jonathan Foley
+ Department of Chemistry, William Paterson University
+ 300 Pompton Road
+ Wayne NJ 07470
+
+*/
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+// Readability
+
+typedef __m256d fourD;
+typedef __m128d twoD;
+
+#define load_4d(ptr) _mm256_load_pd(ptr)
+#define store_4d(ptr, val) _mm256_store_pd(ptr, val)
+#define store_lower_1d(ptr, val) _mm_store_pd(ptr, val)
+#define add_4d(a, b) _mm256_add_pd(a, b)
+#define add_2d(a, b) _mm_add_pd(a, b)
+#define hadd_4d(a, b) _mm256_hadd_pd(a, b)
+#define sub_4d(a, b) _mm256_sub_pd(a, b)
+#define mul_4d(a, b) _mm256_mul_pd(a, b)
+#define div_4d(a, b) _mm256_div_pd(a, b)
+#define create_4d(a) _mm256_set1_pd(a)
+#define get_lower_2d(a) _mm256_castpd256_pd128(a)
+#define get_upper_2d(a) _mm256_extractf128_pd(a, 1)
+
+#define ALIGNMENT 32
+// Number of particles
+int N;
+
+// Lennard-Jones parameters in natural units!
+double sigma = 1.;
+double epsilon = 1.;
+double m = 1.;
+double kB = 1.;
+double epsilon_8 = epsilon * 8.;
+
+double NA = 6.022140857e23;
+double kBSI = 1.38064852e-23; // m^2*kg/(s^2*K)
+
+// Size of box, which will be specified in natural units
+double L;
+
+// Initial Temperature in Natural Units
+double Tinit; // 2;
+// Vectors!
+//
+const int MAXPART = 5001;
+// Position
+double r[MAXPART][3];
+// Velocity
+double v[MAXPART][3];
+// Acceleration
+double a[MAXPART][3];
+// Force
+double F[MAXPART][3];
+
+// atom type
+char atype[10];
+// Function prototypes
+// initialize positions on simple cubic lattice, also calls function to
+// initialize velocities
+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, double *PE, FILE *fp);
+// Compute Force using F = -dV/dr
+// solve F = ma for use in Velocity Verlet
+void computeAccelerations();
+// Numerical Recipes function for generation gaussian distribution
+double gaussdist();
+// Initialize velocities according to user-supplied initial Temperature
+// (Tinit)
+void initializeVelocities();
+// Compute total potential energy from particle coordinates
+double Potential();
+// Compute mean squared velocity from particle velocities
+double MeanSquaredVelocity();
+// Compute total kinetic energy from particle mass and velocities
+double Kinetic();
+
+int main() {
+
+ // variable delcarations
+ int i;
+ double dt, Vol, Temp, Press, Pavg, Tavg, rho;
+ double VolFac, TempFac, PressFac, timefac;
+ double KE, PE, mvs, gc, Z;
+ char trash[10000], prefix[1000], tfn[1000], ofn[1000], afn[1000];
+ FILE *infp, *tfp, *ofp, *afp;
+
+ printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+ printf(" WELCOME TO WILLY P CHEM MD!\n");
+ printf(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+ printf("\n ENTER A TITLE FOR YOUR CALCULATION!\n");
+ scanf("%s", prefix);
+ strcpy(tfn, prefix);
+ strcat(tfn, "_traj.xyz");
+ strcpy(ofn, prefix);
+ strcat(ofn, "_output.txt");
+ strcpy(afn, prefix);
+ strcat(afn, "_average.txt");
+
+ printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+ printf(" TITLE ENTERED AS '%s'\n", prefix);
+ printf(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+
+ /* Table of values for Argon relating natural units to SI units:
+ * These are derived from Lennard-Jones parameters from the article
+ * "Liquid argon: Monte carlo and molecular dynamics calculations"
+ * J.A. Barker , R.A. Fisher & R.O. Watts
+ * Mol. Phys., Vol. 21, 657-673 (1971)
+ *
+ * mass: 6.633e-26 kg = one natural unit of mass for
+ *argon, by definition energy: 1.96183e-21 J = one natural unit of
+ *energy for argon, directly from L-J parameters length: 3.3605e-10 m =
+ *one natural unit of length for argon, directly from L-J parameters
+ * volume: 3.79499-29 m^3 = one natural unit of volume for
+ *argon, by length^3 time: 1.951e-12 s = one natural unit of
+ *time for argon, by length*sqrt(mass/energy)
+ ***************************************************************************************/
+
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // Edit these factors to be computed in terms of basic properties in
+ // natural units of the gas being simulated
+
+ printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\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(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+ scanf("%s", atype);
+
+ if (strcmp(atype, "He") == 0) {
+
+ VolFac = 1.8399744000000005e-29;
+ PressFac = 8152287.336171632;
+ TempFac = 10.864459551225972;
+ timefac = 1.7572698825166272e-12;
+ } else if (strcmp(atype, "Ne") == 0) {
+
+ VolFac = 2.0570823999999997e-29;
+ PressFac = 27223022.27659913;
+ TempFac = 40.560648991243625;
+ timefac = 2.1192341945685407e-12;
+ } else if (strcmp(atype, "Ar") == 0) {
+
+ VolFac = 3.7949992920124995e-29;
+ PressFac = 51695201.06691862;
+ TempFac = 142.0950000000000;
+ timefac = 2.09618e-12;
+ // strcpy(atype,"Ar");
+ } else if (strcmp(atype, "Kr") == 0) {
+
+ VolFac = 4.5882712000000004e-29;
+ PressFac = 59935428.40275003;
+ TempFac = 199.1817584391428;
+ timefac = 8.051563913585078e-13;
+ } else if (strcmp(atype, "Xe") == 0) {
+
+ VolFac = 5.4872e-29;
+ PressFac = 70527773.72794868;
+ TempFac = 280.30305642163006;
+ timefac = 9.018957925790732e-13;
+ } else {
+
+ VolFac = 3.7949992920124995e-29;
+ PressFac = 51695201.06691862;
+ TempFac = 142.0950000000000;
+ timefac = 2.09618e-12;
+ strcpy(atype, "Ar");
+ }
+ printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+ printf("\n YOU ARE SIMULATING %s GAS! \n", atype);
+ printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+
+ printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+ printf("\n YOU WILL NOW ENTER A FEW SIMULATION PARAMETERS\n");
+ printf(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
+ printf("\n\n ENTER THE INTIAL TEMPERATURE OF YOUR GAS IN KELVIN\n");
+ scanf("%lf", &Tinit);
+ // Make sure temperature is a positive number!
+ if (Tinit < 0.) {
+ printf("\n !!!!! ABSOLUTE TEMPERATURE MUST BE A POSITIVE "
+ "NUMBER! PLEASE "
+ "TRY AGAIN WITH A POSITIVE TEMPERATURE!!!\n");
+ exit(0);
+ }
+ // Convert initial temperature from kelvin to natural units
+ Tinit /= TempFac;
+
+ printf("\n\n ENTER THE NUMBER DENSITY IN moles/m^3\n");
+ printf(" FOR REFERENCE, NUMBER DENSITY OF AN IDEAL GAS AT STP IS ABOUT 40 "
+ "moles/m^3\n");
+ printf(" NUMBER DENSITY OF LIQUID ARGON AT 1 ATM AND 87 K IS ABOUT 35000 "
+ "moles/m^3\n");
+
+ scanf("%lf", &rho);
+
+ N = 10 * 500;
+ Vol = N / (rho * NA);
+
+ Vol /= VolFac;
+
+ // Limiting N to MAXPART for practical reasons
+ if (N >= MAXPART) {
+
+ printf("\n\n\n MAXIMUM NUMBER OF PARTICLES IS %i\n\n PLEASE "
+ "ADJUST YOUR "
+ "INPUT FILE ACCORDINGLY \n\n",
+ MAXPART);
+ exit(0);
+ }
+ // Check to see if the volume makes sense - is it too small?
+ // Remember VDW radius of the particles is 1 natural unit of length
+ // and volume = L*L*L, so if V = N*L*L*L = N, then all the particles
+ // will be initialized with an interparticle separation equal to 2xVDW
+ // radius
+ if (Vol < N) {
+
+ printf("\n\n\n YOUR DENSITY IS VERY HIGH!\n\n");
+ printf(" THE NUMBER OF PARTICLES IS %i AND THE AVAILABLE VOLUME "
+ "IS %f "
+ "NATURAL UNITS\n",
+ N, Vol);
+ printf(" SIMULATIONS WITH DENSITY GREATER THAN 1 PARTCICLE/(1 "
+ "Natural "
+ "Unit of Volume) MAY DIVERGE\n");
+ printf(" PLEASE ADJUST YOUR INPUT FILE ACCORDINGLY AND RETRY\n\n");
+ exit(0);
+ }
+ // Vol = L*L*L;
+ // Length of the box in natural units:
+ L = pow(Vol, (1. / 3));
+
+ // Files that we can write different quantities to
+ tfp = fopen(tfn, "w"); // The MD trajectory, coordinates of every
+ // particle at each timestep
+ ofp = fopen(ofn,
+ "w"); // Output of other quantities (T, P, gc, etc) at every timestep
+ afp = fopen(afn, "w"); // Average T, P, gc, etc from the simulation
+
+ int NumTime;
+ if (strcmp(atype, "He") == 0) {
+
+ // dt in natural units of time s.t. in SI it is 5 f.s. for all
+ // other gasses
+ dt = 0.2e-14 / timefac;
+ // We will run the simulation for NumTime timesteps.
+ // The total time will be NumTime*dt in natural units
+ // And NumTime*dt multiplied by the appropriate conversion factor
+ // for time in seconds
+ NumTime = 50000;
+ } else {
+ dt = 0.5e-14 / timefac;
+ NumTime = 200;
+ }
+
+ // Put all the atoms in simple crystal lattice and give them random
+ // velocities that corresponds to the initial temperature we have
+ // specified
+ initialize();
+
+ // Based on their positions, calculate the ininial intermolecular forces
+ // The accellerations of each particle will be defined from the forces and
+ // their mass, and this will allow us to update their positions via
+ // Newton's law
+ computeAccelerations();
+
+ // Print number of particles to the trajectory file
+ fprintf(tfp, "%i\n", N);
+
+ // We want to calculate the average Temperature and Pressure for the
+ // simulation The variables need to be set to zero initially
+ Pavg = 0;
+ Tavg = 0;
+
+ int tenp = floor(NumTime / 10);
+ fprintf(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 [");
+ for (i = 0; i < NumTime + 1; i++) {
+
+ // This just prints updates on progress of the calculation for the
+ // users convenience
+ if (i == tenp)
+ printf(" 10 |");
+ else if (i == 2 * tenp)
+ printf(" 20 |");
+ else if (i == 3 * tenp)
+ printf(" 30 |");
+ else if (i == 4 * tenp)
+ printf(" 40 |");
+ else if (i == 5 * tenp)
+ printf(" 50 |");
+ else if (i == 6 * tenp)
+ printf(" 60 |");
+ else if (i == 7 * tenp)
+ printf(" 70 |");
+ else if (i == 8 * tenp)
+ printf(" 80 |");
+ else if (i == 9 * tenp)
+ printf(" 90 |");
+ else if (i == 10 * tenp)
+ printf(" 100 ]\n");
+ fflush(stdout);
+
+ // This updates the positions and velocities using Newton's Laws
+ // 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, &PE, tfp);
+ Press *= PressFac;
+
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // Now we would like to calculate somethings about the system:
+ // Instantaneous mean velocity squared, Temperature, Pressure
+ // Potential, and Kinetic Energy
+ // We would also like to use the IGL to try to see if we can
+ // extract the gas constant
+ mvs = MeanSquaredVelocity();
+ KE = Kinetic();
+
+ // Temperature from Kinetic Theory
+ Temp = m * mvs / (3 * kB) * TempFac;
+
+ // Instantaneous gas constant and compressibility - not well
+ // defined because pressure may be zero in some instances because
+ // there will be zero wall collisions, pressure may be very high in
+ // some instances because there will be a number of collisions
+ gc = NA * Press * (Vol * VolFac) / (N * Temp);
+ Z = Press * (Vol * VolFac) / (N * kBSI * Temp);
+
+ Tavg += Temp;
+ Pavg += Press;
+
+ fprintf(ofp, " %8.4e %20.12f %20.12f %20.12f %20.12f %20.12f \n", i * dt * timefac,
+ Temp, Press, KE, PE, KE + PE);
+ }
+
+ // Because we have calculated the instantaneous temperature and pressure,
+ // we can take the average over the whole simulation here
+ Pavg /= NumTime;
+ Tavg /= NumTime;
+ Z = Pavg * (Vol * VolFac) / (N * kBSI * Tavg);
+ gc = NA * Pavg * (Vol * VolFac) / (N * Tavg);
+ fprintf(afp, " Total Time (s) T (K) P (Pa) PV/nT "
+ "(J/(mol K)) Z V (m^3) N\n");
+ fprintf(afp, " -------------- ----------- --------------- "
+ "-------------- --------------- ------------ -----------\n");
+ fprintf(afp,
+ " %8.4e %15.5f %15.5f %10.5f %10.5f %10.5e "
+ " %i\n",
+ i * dt * timefac, Tavg, Pavg, gc, Z, Vol * VolFac, N);
+
+ printf("\n TO ANIMATE YOUR SIMULATION, OPEN THE FILE \n '%s' WITH VMD "
+ "AFTER THE SIMULATION COMPLETES\n",
+ tfn);
+ printf("\n TO ANALYZE INSTANTANEOUS DATA ABOUT YOUR MOLECULE, OPEN THE FILE "
+ "\n "
+ " '%s' WITH YOUR FAVORITE TEXT EDITOR OR IMPORT THE DATA INTO EXCEL\n",
+ ofn);
+ printf("\n THE FOLLOWING THERMODYNAMIC AVERAGES WILL BE COMPUTED AND "
+ "WRITTEN TO THE FILE \n '%s':\n",
+ afn);
+ printf("\n AVERAGE TEMPERATURE (K): %15.5f\n", Tavg);
+ printf("\n AVERAGE PRESSURE (Pa): %15.5f\n", Pavg);
+ printf("\n PV/nT (J * mol^-1 K^-1): %15.5f\n", gc);
+ printf("\n PERCENT ERROR of pV/nT AND GAS CONSTANT: %15.5f\n",
+ 100 * fabs(gc - 8.3144598) / 8.3144598);
+ printf("\n THE COMPRESSIBILITY (unitless): %15.5f \n", Z);
+ printf("\n TOTAL VOLUME (m^3): %10.5e \n", Vol * VolFac);
+ printf("\n NUMBER OF PARTICLES (unitless): %i \n", N);
+
+ fclose(tfp);
+ fclose(ofp);
+ fclose(afp);
+
+ return 0;
+}
+
+void initialize() {
+ int n, p, i, j, k;
+ double pos;
+
+ // Number of atoms in each direction
+ n = int(ceil(pow(N, 1.0 / 3)));
+
+ // spacing between atoms along a given direction
+ pos = L / n;
+
+ // index for number of particles assigned positions
+ p = 0;
+ // initialize positions
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < n; j++) {
+ for (k = 0; k < n; k++) {
+ if (p < N) {
+
+ r[p][0] = (i + 0.5) * pos;
+ r[p][1] = (j + 0.5) * pos;
+ r[p][2] = (k + 0.5) * pos;
+ }
+ p++;
+ }
+ }
+ }
+
+ // Call function to initialize velocities
+ initializeVelocities();
+
+ /***********************************************
+ * Uncomment if you want to see what the initial positions and velocities
+ are printf(" Printing initial positions!\n"); for (i=0; i= N) {
+// break;
+// }
+// double quot, rnorm, term1, term2;
+//
+// fourD rj = load_4d(r[j + k]);
+// fourD rij = sub_4d(ri, rj);
+// fourD distTmp = mul_4d(rij, rij);
+// rij4[k] = rij;
+//
+// fourD temp = hadd_4d(distTmp, distTmp);
+// twoD lo = get_lower_2d(temp);
+// twoD hi = get_upper_2d(temp);
+// twoD sum = add_2d(lo, hi);
+//
+// // Convert the scalar result to a double
+// double finalSum;
+// store_lower_1d(&dist[k], sum);
+// }
+// fourD dists = load_4d(dist);
+// fourD quot = div_4d(create_4d(sigma * sigma), dists);
+// fourD term2 = mul_4d(quot, mul_4d(quot, quot));
+// fourD Pots = mul_4d(create_4d(epsilon_8), mul_4d(term2, sub_4d(term2, create_4d(1.))));
+// // Perform horizontal addition and store the result in a scalar
+// fourD temp = hadd_4d(Pots, Pots);
+// twoD lo = get_lower_2d(temp);
+// twoD hi = get_upper_2d(temp);
+// twoD sum = add_2d(lo, hi);
+//
+// // Convert the scalar result to a double
+// double finalSum;
+// store_lower_1d(&finalSum, sum);
+// Pot += finalSum;
+//
+// fourD distSqd = mul_4d(dists, mul_4d(dists, dists));
+// fourD rSqd_inv7 = mul_4d(distSqd, mul_4d(distSqd, dists));
+// fourD f8 = div_4d(sub_4d(create_4d(48.), mul_4d(create_4d(24.), distSqd)), rSqd_inv7);
+//
+// // Go back to the original loop
+// for (int k = 0; k < numParts; k++) {
+// if (j + k >= N) {
+// break;
+// }
+// fourD tmp = mul_4d(rij4[k], create_4d(f8[k]));
+// fourD aI = load_4d(a[i]);
+// fourD aJ = load_4d(a[j + k]);
+// fourD sum = add_4d(aI, tmp);
+// fourD sub = sub_4d(aJ, tmp);
+// store_4d(a[i], sum);
+// store_4d(a[j + k], sub);
+// }
+// }
+// }
+// //}
+//
+// return Pot;
+// }
+
+double PotentialAndAcceleration(double dt) {
+ memset(a, 0, sizeof(a));
+ double Pot = 0.;
+ # pragma omp parallel for reduction(+:Pot,a)
+ for (int i = 0; i < N - 1; i++) {
+ double ai0 = 0., ai1 = 0., ai2 = 0.;
+ double rI[3] = {r[i][0], r[i][1], r[i][2]};
+ for (int j = i + 1; j < N; j++) {
+ double quot, term2;
+ // component-by-componenent position of i relative to j
+ double rij[3];
+ // sum of squares of the components
+ double rSqd = 0.;
+
+ for (int k = 0; k < 3; k++) {
+ double rijk = rI[k] - r[j][k];
+
+ rij[k] = rijk;
+
+ rSqd += rijk * rijk;
+ }
+
+ // Here we remove the pow function and simplify the calculation
+ double rSqd_3 = rSqd * rSqd * rSqd;
+ double rSqd_7 = rSqd_3 * rSqd_3 * rSqd;
+ double f = (48. - (24. * rSqd_3)) / rSqd_7;
+ // from F = ma, where m = 1 in natural units!
+ double tmp = rij[0] * f;
+ double tmp1 = rij[1] * f;
+ double tmp2 = rij[2] * f;
+ ai0 += tmp;
+ ai1 += tmp1;
+ ai2 += tmp2;
+ a[j][0] -= tmp;
+ a[j][1] -= tmp1;
+ a[j][2] -= tmp2;
+
+ quot = sigma * sigma / rSqd;
+ term2 = quot * quot * quot;
+ Pot += term2 * (term2 - 1.);
+ }
+ a[i][0] += ai0;
+ a[i][1] += ai1;
+ a[i][2] += ai2;
+ }
+ return Pot*epsilon_8;
+}
+
+// Function to calculate the potential energy of the system
+double Potential() {
+ double quot, rSqd, rnorm, term1, term2, Pot;
+ int i, j;
+
+ Pot = 0.;
+
+ for (i = 0; i < N - 1; i++) {
+ for (j = i + 1; j < N; j++) {
+ rSqd = 0.;
+ double rI[3] = {r[i][0], r[i][1], r[i][2]};
+ for (int k = 0; k < 3; k++) {
+ double rijk = rI[k] - r[j][k];
+ rSqd += rijk * rijk;
+ }
+
+ quot = sigma * sigma / rSqd;
+ term2 = quot * quot * quot;
+ Pot += epsilon_8 * term2 * (term2 - 1.);
+ }
+ }
+ return Pot;
+}
+
+// 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 computeAccelerations() {
+ double f, rSqd, tmp = 0.;
+ memset(a, 0, sizeof(a));
+
+ for (int i = 0; i < N - 1; i++) { // loop over all distinct pairs i,j
+ double rI[3] = {r[i][0], r[i][1], r[i][2]};
+ for (int j = i + 1; j < N; j++) {
+ // component-by-componenent position of i relative to j
+ double rij[3];
+ // sum of squares of the components
+ double rSqd = 0.;
+
+ for (int k = 0; k < 3; k++) {
+ double rijk = rI[k] - r[j][k];
+
+ rij[k] = rijk;
+
+ rSqd += rijk * rijk;
+ }
+
+ // Here we remove the pow function and simplify the calculation
+ double rSqd_3 = rSqd * rSqd * rSqd;
+ double rSqd_7 = rSqd_3 * rSqd_3 * rSqd;
+ double f = (48. - (24. * rSqd_3)) / rSqd_7;
+ // from F = ma, where m = 1 in natural units!
+ for (int k = 0; k < 3; k++) {
+ double tmp = rij[k] * f;
+ a[i][k] += tmp;
+ a[j][k] -= tmp;
+ }
+ }
+ }
+}
+
+// returns sum of dv/dt*m/A (aka Pressure) from elastic collisions with walls
+double VelocityVerlet(double dt, int iter, double *PE, FILE *fp) {
+ int i, j, k;
+
+ double psum = 0.;
+
+ // Compute accelerations from forces at current position
+ // this call was removed (commented) for predagogical reasons
+ // computeAccelerations();
+ // Update positions and velocity with current velocity and acceleration
+ // printf(" Updated Positions!\n");
+ for (i = 0; i < N; i++) {
+ for (j = 0; j < 3; j++) {
+ double tmp = 0.5 * a[i][j] * dt;
+ r[i][j] += v[i][j] * dt + tmp * 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(dt);
+ // Update velocity with updated acceleration
+ for (i = 0; i < N; i++) {
+ for (j = 0; j < 3; j++) {
+ v[i][j] += 0.5 * a[i][j] * dt;
+ }
+ }
+
+ // Elastic walls
+ for (i = 0; i < N; i++) {
+ for (j = 0; j < 3; j++) {
+ if (r[i][j] < 0.) {
+ v[i][j] *= -1.; //- elastic walls
+ psum += 2 * m * fabs(v[i][j]) / 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]) / dt; // contribution to pressure
+ // from "right" walls
+ }
+ }
+ }
+
+ /* removed, uncomment to save atoms positions */
+ /*for (i=0; i= 1.0 || rsq == 0.0);
+
+ fac = sqrt(-2.0 * log(rsq) / rsq);
+ gset = v1 * fac;
+ available = true;
+
+ return v2 * fac;
+ } else {
+
+ available = false;
+ return gset;
+ }
+}