gcc flags changed

This commit is contained in:
Afonso Franco 2023-10-16 19:51:28 +01:00
parent 4c27b3c774
commit 763c84739f
Signed by: afonso
SSH key fingerprint: SHA256:JiuxZNdA5bRWXPMUJChI0AQ75yC+cXY4xM0IaVwEVys
3 changed files with 526 additions and 538 deletions

View file

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

BIN
out/MD

Binary file not shown.

View file

@ -65,28 +65,27 @@ char atype[10];
// Function prototypes // Function prototypes
// initialize positions on simple cubic lattice, also calls function to // initialize positions on simple cubic lattice, also calls function to
// initialize velocities // initialize velocities
void initialize (); void initialize();
// update positions and velocities using Velocity Verlet algorithm // update positions and velocities using Velocity Verlet algorithm
// print particle coordinates to file for rendering via VMD or other animation // print particle coordinates to file for rendering via VMD or other animation
// software return 'instantaneous pressure' // software return 'instantaneous pressure'
double VelocityVerlet (double dt, int iter, double *PE, FILE *fp); double VelocityVerlet(double dt, int iter, double *PE, FILE *fp);
// Compute Force using F = -dV/dr // Compute Force using F = -dV/dr
// solve F = ma for use in Velocity Verlet // solve F = ma for use in Velocity Verlet
void computeAccelerations (); void computeAccelerations();
// Numerical Recipes function for generation gaussian distribution // Numerical Recipes function for generation gaussian distribution
double gaussdist (); double gaussdist();
// Initialize velocities according to user-supplied initial Temperature // Initialize velocities according to user-supplied initial Temperature
// (Tinit) // (Tinit)
void initializeVelocities (); void initializeVelocities();
// Compute total potential energy from particle coordinates // Compute total potential energy from particle coordinates
double Potential (); double Potential();
// Compute mean squared velocity from particle velocities // Compute mean squared velocity from particle velocities
double MeanSquaredVelocity (); double MeanSquaredVelocity();
// Compute total kinetic energy from particle mass and velocities // Compute total kinetic energy from particle mass and velocities
double Kinetic (); double Kinetic();
int int main() {
main () {
// variable delcarations // variable delcarations
int i; int i;
@ -96,21 +95,21 @@ main () {
char trash[10000], prefix[1000], tfn[1000], ofn[1000], afn[1000]; char trash[10000], prefix[1000], tfn[1000], ofn[1000], afn[1000];
FILE *infp, *tfp, *ofp, *afp; FILE *infp, *tfp, *ofp, *afp;
printf ("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
printf (" WELCOME TO WILLY P CHEM MD!\n"); printf(" WELCOME TO WILLY P CHEM MD!\n");
printf (" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
printf ("\n ENTER A TITLE FOR YOUR CALCULATION!\n"); printf("\n ENTER A TITLE FOR YOUR CALCULATION!\n");
scanf ("%s", prefix); scanf("%s", prefix);
strcpy (tfn, prefix); strcpy(tfn, prefix);
strcat (tfn, "_traj.xyz"); strcat(tfn, "_traj.xyz");
strcpy (ofn, prefix); strcpy(ofn, prefix);
strcat (ofn, "_output.txt"); strcat(ofn, "_output.txt");
strcpy (afn, prefix); strcpy(afn, prefix);
strcat (afn, "_average.txt"); strcat(afn, "_average.txt");
printf ("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
printf (" TITLE ENTERED AS '%s'\n", prefix); printf(" TITLE ENTERED AS '%s'\n", prefix);
printf (" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
/* Table of values for Argon relating natural units to SI units: /* Table of values for Argon relating natural units to SI units:
* These are derived from Lennard-Jones parameters from the article * These are derived from Lennard-Jones parameters from the article
@ -131,42 +130,42 @@ main () {
// Edit these factors to be computed in terms of basic properties in // Edit these factors to be computed in terms of basic properties in
// natural units of the gas being simulated // natural units of the gas being simulated
printf ("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf("\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("\n FOR HELIUM, TYPE 'He' THEN PRESS 'return' TO CONTINUE\n");
printf (" FOR NEON, TYPE 'Ne' 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 ARGON, TYPE 'Ar' THEN PRESS 'return' TO CONTINUE\n");
printf (" FOR KRYPTON, TYPE 'Kr' 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(" FOR XENON, TYPE 'Xe' THEN PRESS 'return' TO CONTINUE\n");
printf (" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
scanf ("%s", atype); scanf("%s", atype);
if (strcmp (atype, "He") == 0) { if (strcmp(atype, "He") == 0) {
VolFac = 1.8399744000000005e-29; VolFac = 1.8399744000000005e-29;
PressFac = 8152287.336171632; PressFac = 8152287.336171632;
TempFac = 10.864459551225972; TempFac = 10.864459551225972;
timefac = 1.7572698825166272e-12; timefac = 1.7572698825166272e-12;
} else if (strcmp (atype, "Ne") == 0) { } else if (strcmp(atype, "Ne") == 0) {
VolFac = 2.0570823999999997e-29; VolFac = 2.0570823999999997e-29;
PressFac = 27223022.27659913; PressFac = 27223022.27659913;
TempFac = 40.560648991243625; TempFac = 40.560648991243625;
timefac = 2.1192341945685407e-12; timefac = 2.1192341945685407e-12;
} else if (strcmp (atype, "Ar") == 0) { } else if (strcmp(atype, "Ar") == 0) {
VolFac = 3.7949992920124995e-29; VolFac = 3.7949992920124995e-29;
PressFac = 51695201.06691862; PressFac = 51695201.06691862;
TempFac = 142.0950000000000; TempFac = 142.0950000000000;
timefac = 2.09618e-12; timefac = 2.09618e-12;
// strcpy(atype,"Ar"); // strcpy(atype,"Ar");
} else if (strcmp (atype, "Kr") == 0) { } else if (strcmp(atype, "Kr") == 0) {
VolFac = 4.5882712000000004e-29; VolFac = 4.5882712000000004e-29;
PressFac = 59935428.40275003; PressFac = 59935428.40275003;
TempFac = 199.1817584391428; TempFac = 199.1817584391428;
timefac = 8.051563913585078e-13; timefac = 8.051563913585078e-13;
} else if (strcmp (atype, "Xe") == 0) { } else if (strcmp(atype, "Xe") == 0) {
VolFac = 5.4872e-29; VolFac = 5.4872e-29;
PressFac = 70527773.72794868; PressFac = 70527773.72794868;
@ -178,34 +177,34 @@ main () {
PressFac = 51695201.06691862; PressFac = 51695201.06691862;
TempFac = 142.0950000000000; TempFac = 142.0950000000000;
timefac = 2.09618e-12; timefac = 2.09618e-12;
strcpy (atype, "Ar"); strcpy(atype, "Ar");
} }
printf ("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
printf ("\n YOU ARE SIMULATING %s GAS! \n", atype); printf("\n YOU ARE SIMULATING %s GAS! \n", atype);
printf ("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
printf ("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
printf ("\n YOU WILL NOW ENTER A FEW SIMULATION PARAMETERS\n"); printf("\n YOU WILL NOW ENTER A FEW SIMULATION PARAMETERS\n");
printf (" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); printf(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
printf ("\n\n ENTER THE INTIAL TEMPERATURE OF YOUR GAS IN KELVIN\n"); printf("\n\n ENTER THE INTIAL TEMPERATURE OF YOUR GAS IN KELVIN\n");
scanf ("%lf", &Tinit); scanf("%lf", &Tinit);
// Make sure temperature is a positive number! // Make sure temperature is a positive number!
if (Tinit < 0.) { if (Tinit < 0.) {
printf ("\n !!!!! ABSOLUTE TEMPERATURE MUST BE A POSITIVE " printf("\n !!!!! ABSOLUTE TEMPERATURE MUST BE A POSITIVE "
"NUMBER! PLEASE " "NUMBER! PLEASE "
"TRY AGAIN WITH A POSITIVE TEMPERATURE!!!\n"); "TRY AGAIN WITH A POSITIVE TEMPERATURE!!!\n");
exit (0); exit(0);
} }
// Convert initial temperature from kelvin to natural units // Convert initial temperature from kelvin to natural units
Tinit /= TempFac; Tinit /= TempFac;
printf ("\n\n ENTER THE NUMBER DENSITY IN moles/m^3\n"); 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 " printf(" FOR REFERENCE, NUMBER DENSITY OF AN IDEAL GAS AT STP IS ABOUT 40 "
"moles/m^3\n"); "moles/m^3\n");
printf (" NUMBER DENSITY OF LIQUID ARGON AT 1 ATM AND 87 K IS ABOUT 35000 " printf(" NUMBER DENSITY OF LIQUID ARGON AT 1 ATM AND 87 K IS ABOUT 35000 "
"moles/m^3\n"); "moles/m^3\n");
scanf ("%lf", &rho); scanf("%lf", &rho);
N = 10 * 216; N = 10 * 216;
Vol = N / (rho * NA); Vol = N / (rho * NA);
@ -215,11 +214,11 @@ main () {
// Limiting N to MAXPART for practical reasons // Limiting N to MAXPART for practical reasons
if (N >= MAXPART) { if (N >= MAXPART) {
printf ("\n\n\n MAXIMUM NUMBER OF PARTICLES IS %i\n\n PLEASE " printf("\n\n\n MAXIMUM NUMBER OF PARTICLES IS %i\n\n PLEASE "
"ADJUST YOUR " "ADJUST YOUR "
"INPUT FILE ACCORDINGLY \n\n", "INPUT FILE ACCORDINGLY \n\n",
MAXPART); MAXPART);
exit (0); exit(0);
} }
// Check to see if the volume makes sense - is it too small? // Check to see if the volume makes sense - is it too small?
// Remember VDW radius of the particles is 1 natural unit of length // Remember VDW radius of the particles is 1 natural unit of length
@ -228,30 +227,31 @@ main () {
// radius // radius
if (Vol < N) { if (Vol < N) {
printf ("\n\n\n YOUR DENSITY IS VERY HIGH!\n\n"); printf("\n\n\n YOUR DENSITY IS VERY HIGH!\n\n");
printf (" THE NUMBER OF PARTICLES IS %i AND THE AVAILABLE VOLUME " printf(" THE NUMBER OF PARTICLES IS %i AND THE AVAILABLE VOLUME "
"IS %f " "IS %f "
"NATURAL UNITS\n", "NATURAL UNITS\n",
N, Vol); N, Vol);
printf (" SIMULATIONS WITH DENSITY GREATER THAN 1 PARTCICLE/(1 " printf(" SIMULATIONS WITH DENSITY GREATER THAN 1 PARTCICLE/(1 "
"Natural " "Natural "
"Unit of Volume) MAY DIVERGE\n"); "Unit of Volume) MAY DIVERGE\n");
printf (" PLEASE ADJUST YOUR INPUT FILE ACCORDINGLY AND RETRY\n\n"); printf(" PLEASE ADJUST YOUR INPUT FILE ACCORDINGLY AND RETRY\n\n");
exit (0); exit(0);
} }
// Vol = L*L*L; // Vol = L*L*L;
// Length of the box in natural units: // Length of the box in natural units:
L = pow (Vol, (1. / 3)); L = pow(Vol, (1. / 3));
// Files that we can write different quantities to // Files that we can write different quantities to
tfp = fopen (tfn, "w"); // The MD trajectory, coordinates of every tfp = fopen(tfn, "w"); // The MD trajectory, coordinates of every
// particle at each timestep // particle at each timestep
ofp = fopen (ofn, ofp = fopen(
ofn,
"w"); // Output of other quantities (T, P, gc, etc) at every timestep "w"); // Output of other quantities (T, P, gc, etc) at every timestep
afp = fopen (afn, "w"); // Average T, P, gc, etc from the simulation afp = fopen(afn, "w"); // Average T, P, gc, etc from the simulation
int NumTime; int NumTime;
if (strcmp (atype, "He") == 0) { if (strcmp(atype, "He") == 0) {
// dt in natural units of time s.t. in SI it is 5 f.s. for all // dt in natural units of time s.t. in SI it is 5 f.s. for all
// other gasses // other gasses
@ -269,57 +269,58 @@ main () {
// Put all the atoms in simple crystal lattice and give them random // Put all the atoms in simple crystal lattice and give them random
// velocities that corresponds to the initial temperature we have // velocities that corresponds to the initial temperature we have
// specified // specified
initialize (); initialize();
// Based on their positions, calculate the ininial intermolecular forces // Based on their positions, calculate the ininial intermolecular forces
// The accellerations of each particle will be defined from the forces and // The accellerations of each particle will be defined from the forces and
// their mass, and this will allow us to update their positions via // their mass, and this will allow us to update their positions via
// Newton's law // Newton's law
computeAccelerations (); computeAccelerations();
// Print number of particles to the trajectory file // Print number of particles to the trajectory file
fprintf (tfp, "%i\n", N); fprintf(tfp, "%i\n", N);
// We want to calculate the average Temperature and Pressure for the // We want to calculate the average Temperature and Pressure for the
// simulation The variables need to be set to zero initially // simulation The variables need to be set to zero initially
Pavg = 0; Pavg = 0;
Tavg = 0; Tavg = 0;
int tenp = floor (NumTime / 10); int tenp = floor(NumTime / 10);
fprintf (ofp, " time (s) T(t) (K) P(t) (Pa) " fprintf(ofp,
" time (s) T(t) (K) P(t) (Pa) "
"Kinetic En. (n.u.) Potential En. (n.u.) Total En. (n.u.)\n"); "Kinetic En. (n.u.) Potential En. (n.u.) Total En. (n.u.)\n");
printf (" PERCENTAGE OF CALCULATION COMPLETE:\n ["); printf(" PERCENTAGE OF CALCULATION COMPLETE:\n [");
for (i = 0; i < NumTime + 1; i++) { for (i = 0; i < NumTime + 1; i++) {
// This just prints updates on progress of the calculation for the // This just prints updates on progress of the calculation for the
// users convenience // users convenience
if (i == tenp) if (i == tenp)
printf (" 10 |"); printf(" 10 |");
else if (i == 2 * tenp) else if (i == 2 * tenp)
printf (" 20 |"); printf(" 20 |");
else if (i == 3 * tenp) else if (i == 3 * tenp)
printf (" 30 |"); printf(" 30 |");
else if (i == 4 * tenp) else if (i == 4 * tenp)
printf (" 40 |"); printf(" 40 |");
else if (i == 5 * tenp) else if (i == 5 * tenp)
printf (" 50 |"); printf(" 50 |");
else if (i == 6 * tenp) else if (i == 6 * tenp)
printf (" 60 |"); printf(" 60 |");
else if (i == 7 * tenp) else if (i == 7 * tenp)
printf (" 70 |"); printf(" 70 |");
else if (i == 8 * tenp) else if (i == 8 * tenp)
printf (" 80 |"); printf(" 80 |");
else if (i == 9 * tenp) else if (i == 9 * tenp)
printf (" 90 |"); printf(" 90 |");
else if (i == 10 * tenp) else if (i == 10 * tenp)
printf (" 100 ]\n"); printf(" 100 ]\n");
fflush (stdout); fflush(stdout);
// This updates the positions and velocities using Newton's Laws // This updates the positions and velocities using Newton's Laws
// Also computes the Pressure as the sum of momentum changes from // Also computes the Pressure as the sum of momentum changes from
// wall collisions / timestep which is a Kinetic Theory of gasses // wall collisions / timestep which is a Kinetic Theory of gasses
// concept of Pressure // concept of Pressure
Press = VelocityVerlet (dt, i + 1, &PE, tfp); Press = VelocityVerlet(dt, i + 1, &PE, tfp);
Press *= PressFac; Press *= PressFac;
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@ -328,8 +329,8 @@ main () {
// Potential, and Kinetic Energy // Potential, and Kinetic Energy
// We would also like to use the IGL to try to see if we can // We would also like to use the IGL to try to see if we can
// extract the gas constant // extract the gas constant
mvs = MeanSquaredVelocity (); mvs = MeanSquaredVelocity();
KE = Kinetic (); KE = Kinetic();
// Temperature from Kinetic Theory // Temperature from Kinetic Theory
Temp = m * mvs / (3 * kB) * TempFac; Temp = m * mvs / (3 * kB) * TempFac;
@ -344,8 +345,8 @@ main () {
Tavg += Temp; Tavg += Temp;
Pavg += Press; Pavg += Press;
fprintf (ofp, " %8.4e %20.12f %20.12f %20.12f %20.12f %20.12f \n", i * dt * timefac, fprintf(ofp, " %8.4e %20.12f %20.12f %20.12f %20.12f %20.12f \n",
Temp, Press, KE, PE, KE + PE); i * dt * timefac, Temp, Press, KE, PE, KE + PE);
} }
// Because we have calculated the instantaneous temperature and pressure, // Because we have calculated the instantaneous temperature and pressure,
@ -354,48 +355,49 @@ main () {
Tavg /= NumTime; Tavg /= NumTime;
Z = Pavg * (Vol * VolFac) / (N * kBSI * Tavg); Z = Pavg * (Vol * VolFac) / (N * kBSI * Tavg);
gc = NA * Pavg * (Vol * VolFac) / (N * Tavg); gc = NA * Pavg * (Vol * VolFac) / (N * Tavg);
fprintf (afp, " Total Time (s) T (K) P (Pa) PV/nT " fprintf(afp, " Total Time (s) T (K) P (Pa) PV/nT "
"(J/(mol K)) Z V (m^3) N\n"); "(J/(mol K)) Z V (m^3) N\n");
fprintf (afp, " -------------- ----------- --------------- " fprintf(afp,
" -------------- ----------- --------------- "
"-------------- --------------- ------------ -----------\n"); "-------------- --------------- ------------ -----------\n");
fprintf (afp, fprintf(afp,
" %8.4e %15.5f %15.5f %10.5f %10.5f %10.5e " " %8.4e %15.5f %15.5f %10.5f %10.5f %10.5e "
" %i\n", " %i\n",
i * dt * timefac, Tavg, Pavg, gc, Z, Vol * VolFac, N); i * dt * timefac, Tavg, Pavg, gc, Z, Vol * VolFac, N);
printf ("\n TO ANIMATE YOUR SIMULATION, OPEN THE FILE \n '%s' WITH VMD " printf("\n TO ANIMATE YOUR SIMULATION, OPEN THE FILE \n '%s' WITH VMD "
"AFTER THE SIMULATION COMPLETES\n", "AFTER THE SIMULATION COMPLETES\n",
tfn); tfn);
printf ("\n TO ANALYZE INSTANTANEOUS DATA ABOUT YOUR MOLECULE, OPEN THE FILE " printf("\n TO ANALYZE INSTANTANEOUS DATA ABOUT YOUR MOLECULE, OPEN THE FILE "
"\n " "\n "
" '%s' WITH YOUR FAVORITE TEXT EDITOR OR IMPORT THE DATA INTO EXCEL\n", " '%s' WITH YOUR FAVORITE TEXT EDITOR OR IMPORT THE DATA INTO EXCEL\n",
ofn); ofn);
printf ("\n THE FOLLOWING THERMODYNAMIC AVERAGES WILL BE COMPUTED AND " printf("\n THE FOLLOWING THERMODYNAMIC AVERAGES WILL BE COMPUTED AND "
"WRITTEN TO THE FILE \n '%s':\n", "WRITTEN TO THE FILE \n '%s':\n",
afn); afn);
printf ("\n AVERAGE TEMPERATURE (K): %15.5f\n", Tavg); printf("\n AVERAGE TEMPERATURE (K): %15.5f\n", Tavg);
printf ("\n AVERAGE PRESSURE (Pa): %15.5f\n", Pavg); printf("\n AVERAGE PRESSURE (Pa): %15.5f\n", Pavg);
printf ("\n PV/nT (J * mol^-1 K^-1): %15.5f\n", gc); 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", printf("\n PERCENT ERROR of pV/nT AND GAS CONSTANT: %15.5f\n",
100 * fabs (gc - 8.3144598) / 8.3144598); 100 * fabs(gc - 8.3144598) / 8.3144598);
printf ("\n THE COMPRESSIBILITY (unitless): %15.5f \n", Z); printf("\n THE COMPRESSIBILITY (unitless): %15.5f \n", Z);
printf ("\n TOTAL VOLUME (m^3): %10.5e \n", Vol * VolFac); printf("\n TOTAL VOLUME (m^3): %10.5e \n",
printf ("\n NUMBER OF PARTICLES (unitless): %i \n", N); Vol * VolFac);
printf("\n NUMBER OF PARTICLES (unitless): %i \n", N);
fclose (tfp); fclose(tfp);
fclose (ofp); fclose(ofp);
fclose (afp); fclose(afp);
return 0; return 0;
} }
void void initialize() {
initialize () {
int n, p, i, j, k; int n, p, i, j, k;
double pos; double pos;
// Number of atoms in each direction // Number of atoms in each direction
n = int (ceil (pow (N, 1.0 / 3))); n = int(ceil(pow(N, 1.0 / 3)));
// spacing between atoms along a given direction // spacing between atoms along a given direction
pos = L / n; pos = L / n;
@ -418,7 +420,7 @@ initialize () {
} }
// Call function to initialize velocities // Call function to initialize velocities
initializeVelocities (); initializeVelocities();
/*********************************************** /***********************************************
* Uncomment if you want to see what the initial positions and velocities * Uncomment if you want to see what the initial positions and velocities
@ -434,8 +436,7 @@ initialize () {
} }
// Function to calculate the averaged velocity squared // Function to calculate the averaged velocity squared
double double MeanSquaredVelocity() {
MeanSquaredVelocity () {
double vx2 = 0; double vx2 = 0;
double vy2 = 0; double vy2 = 0;
@ -455,8 +456,7 @@ MeanSquaredVelocity () {
} }
// Function to calculate the kinetic energy of the system // Function to calculate the kinetic energy of the system
double double Kinetic() { // Write Function here!
Kinetic () { // Write Function here!
double v2, kin; double v2, kin;
@ -475,12 +475,6 @@ Kinetic () { // Write Function here!
return kin; return kin;
} }
double
getF (double dist) {}
double
getLocalPot (double dist) {}
// 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++) {
@ -490,9 +484,8 @@ getLocalPot (double dist) {}
// } // }
// } // }
double double PotentialAndAcceleration() {
PotentialAndAcceleration () { memset(a, 0, sizeof(a));
memset (a, 0, sizeof (a));
double Pot = 0.; double Pot = 0.;
// double rT[3][MAXPART]; // double rT[3][MAXPART];
// transposeMatrix (r, rT); // transposeMatrix (r, rT);
@ -521,7 +514,7 @@ PotentialAndAcceleration () {
double distSqd = dist * dist; double distSqd = dist * dist;
double rSqd_inv4 = 1.0 / (distSqd * distSqd); double rSqd_inv4 = 1.0 / (distSqd * distSqd);
double rSqd_inv7 = rSqd_inv4 / (distSqd * dist); double rSqd_inv7 = rSqd_inv4 / (distSqd * dist);
double f = 24.0 * (2.0 * rSqd_inv7 - rSqd_inv4); double f = 48.0 * rSqd_inv7 - 24 * rSqd_inv4;
// from F = ma, where m = 1 in natural units! // from F = ma, where m = 1 in natural units!
for (int k = 0; k < 3; k++) { for (int k = 0; k < 3; k++) {
double tmp = posItoJ[k] * f; double tmp = posItoJ[k] * f;
@ -534,8 +527,7 @@ PotentialAndAcceleration () {
} }
// Function to calculate the potential energy of the system // Function to calculate the potential energy of the system
double double Potential() {
Potential () {
double quot, r2, rnorm, term1, term2, Pot; double quot, r2, rnorm, term1, term2, Pot;
int i, j; int i, j;
@ -550,7 +542,7 @@ Potential () {
r2 += tmp * tmp; r2 += tmp * tmp;
} }
quot = sigma / sqrt (r2); quot = sigma / sqrt(r2);
term2 = quot * quot; term2 = quot * quot;
term2 = term2 * term2 * term2; term2 = term2 * term2 * term2;
term1 = term2 * term2; term1 = term2 * term2;
@ -564,13 +556,12 @@ Potential () {
// Uses the derivative of the Lennard-Jones potential to calculate // Uses the derivative of the Lennard-Jones potential to calculate
// the forces on each atom. Then uses a = F/m to calculate the // the forces on each atom. Then uses a = F/m to calculate the
// accelleration of each atom. // accelleration of each atom.
void void computeAccelerations() {
computeAccelerations () {
int i, j, k; int i, j, k;
double f, rSqd, tmp = 0.; double f, rSqd, tmp = 0.;
double rij[3]; // position of i relative to j double rij[3]; // position of i relative to j
memset (a, 0, sizeof (a)); memset(a, 0, sizeof(a));
for (i = 0; i < N - 1; i++) { // loop over all distinct pairs i,j for (i = 0; i < N - 1; i++) { // loop over all distinct pairs i,j
for (j = i + 1; j < N; j++) { for (j = i + 1; j < N; j++) {
// initialize r^2 to zero // initialize r^2 to zero
@ -600,8 +591,7 @@ computeAccelerations () {
} }
// returns sum of dv/dt*m/A (aka Pressure) from elastic collisions with walls // 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; int i, j, k;
double psum = 0.; double psum = 0.;
@ -621,7 +611,7 @@ VelocityVerlet (double dt, int iter, double *PE, FILE *fp) {
} }
// Update accellerations from updated positions // Update accellerations from updated positions
// computeAccelerations (); // computeAccelerations ();
*PE = PotentialAndAcceleration (); *PE = PotentialAndAcceleration();
// Update velocity with updated acceleration // 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++) {
@ -634,12 +624,12 @@ VelocityVerlet (double dt, int iter, double *PE, FILE *fp) {
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
if (r[i][j] < 0.) { if (r[i][j] < 0.) {
v[i][j] *= -1.; //- elastic walls v[i][j] *= -1.; //- elastic walls
psum += 2 * m * fabs (v[i][j]) / dt; // contribution to pressure psum += 2 * m * fabs(v[i][j]) / dt; // contribution to pressure
// from "left" walls // from "left" walls
} }
if (r[i][j] >= L) { if (r[i][j] >= L) {
v[i][j] *= -1.; //- elastic walls v[i][j] *= -1.; //- elastic walls
psum += 2 * m * fabs (v[i][j]) / dt; // contribution to pressure psum += 2 * m * fabs(v[i][j]) / dt; // contribution to pressure
// from "right" walls // from "right" walls
} }
} }
@ -658,8 +648,7 @@ VelocityVerlet (double dt, int iter, double *PE, FILE *fp) {
return psum / (6 * L * L); return psum / (6 * L * L);
} }
void void initializeVelocities() {
initializeVelocities () {
int i, j; int i, j;
@ -667,13 +656,13 @@ initializeVelocities () {
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
// Pull a number from a Gaussian Distribution // Pull a number from a Gaussian Distribution
v[i][j] = gaussdist (); v[i][j] = gaussdist();
} }
} }
// Vcm = sum_i^N m*v_i/ sum_i^N M // Vcm = sum_i^N m*v_i/ sum_i^N M
// Compute center-of-mas velocity according to the formula above // Compute center-of-mas velocity according to the formula above
double vCM[3] = { 0, 0, 0 }; double vCM[3] = {0, 0, 0};
for (i = 0; i < N; i++) { for (i = 0; i < N; i++) {
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
@ -707,7 +696,7 @@ initializeVelocities () {
} }
} }
lambda = sqrt (3 * (N - 1) * Tinit / vSqdSum); lambda = sqrt(3 * (N - 1) * Tinit / vSqdSum);
for (i = 0; i < N; i++) { for (i = 0; i < N; i++) {
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
@ -718,19 +707,18 @@ initializeVelocities () {
} }
// Numerical recipes Gaussian distribution number generator // Numerical recipes Gaussian distribution number generator
double double gaussdist() {
gaussdist () {
static bool available = false; static bool available = false;
static double gset; static double gset;
double fac, rsq, v1, v2; double fac, rsq, v1, v2;
if (!available) { if (!available) {
do { do {
v1 = 2.0 * rand () / double (RAND_MAX) - 1.0; v1 = 2.0 * rand() / double(RAND_MAX) - 1.0;
v2 = 2.0 * rand () / double (RAND_MAX) - 1.0; v2 = 2.0 * rand() / double(RAND_MAX) - 1.0;
rsq = v1 * v1 + v2 * v2; rsq = v1 * v1 + v2 * v2;
} while (rsq >= 1.0 || rsq == 0.0); } while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt (-2.0 * log (rsq) / rsq); fac = sqrt(-2.0 * log(rsq) / rsq);
gset = v1 * fac; gset = v1 * fac;
available = true; available = true;