Removed Sqrt function and optimized multiplications
This commit is contained in:
parent
44e386f674
commit
4c27b3c774
SSH key fingerprint:
SHA256:JiuxZNdA5bRWXPMUJChI0AQ75yC+cXY4xM0IaVwEVys
9 changed files with 626 additions and 268 deletions
167
src/MD.cpp
167
src/MD.cpp
|
|
@ -38,6 +38,7 @@ 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)
|
||||
|
|
@ -95,11 +96,9 @@ main () {
|
|||
char trash[10000], prefix[1000], tfn[1000], ofn[1000], afn[1000];
|
||||
FILE *infp, *tfp, *ofp, *afp;
|
||||
|
||||
printf (
|
||||
"\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
printf ("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
printf (" WELCOME TO WILLY P CHEM MD!\n");
|
||||
printf (
|
||||
" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
printf (" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
printf ("\n ENTER A TITLE FOR YOUR CALCULATION!\n");
|
||||
scanf ("%s", prefix);
|
||||
strcpy (tfn, prefix);
|
||||
|
|
@ -109,11 +108,9 @@ main () {
|
|||
strcpy (afn, prefix);
|
||||
strcat (afn, "_average.txt");
|
||||
|
||||
printf (
|
||||
"\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
printf ("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
printf (" TITLE ENTERED AS '%s'\n", prefix);
|
||||
printf (
|
||||
" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
printf (" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
|
||||
/* Table of values for Argon relating natural units to SI units:
|
||||
* These are derived from Lennard-Jones parameters from the article
|
||||
|
|
@ -134,17 +131,14 @@ main () {
|
|||
// 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\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");
|
||||
printf (" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
scanf ("%s", atype);
|
||||
|
||||
if (strcmp (atype, "He") == 0) {
|
||||
|
|
@ -186,17 +180,13 @@ main () {
|
|||
timefac = 2.09618e-12;
|
||||
strcpy (atype, "Ar");
|
||||
}
|
||||
printf (
|
||||
"\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
printf ("\n !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
|
||||
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");
|
||||
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!
|
||||
|
|
@ -210,12 +200,10 @@ main () {
|
|||
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");
|
||||
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);
|
||||
|
||||
|
|
@ -258,9 +246,8 @@ main () {
|
|||
// 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
|
||||
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;
|
||||
|
|
@ -299,10 +286,8 @@ main () {
|
|||
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");
|
||||
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++) {
|
||||
|
||||
|
|
@ -359,8 +344,8 @@ main () {
|
|||
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);
|
||||
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,
|
||||
|
|
@ -369,27 +354,22 @@ main () {
|
|||
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,
|
||||
" 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);
|
||||
" %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 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);
|
||||
|
|
@ -399,8 +379,7 @@ main () {
|
|||
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 TOTAL VOLUME (m^3): %10.5e \n", Vol * VolFac);
|
||||
printf ("\n NUMBER OF PARTICLES (unitless): %i \n", N);
|
||||
|
||||
fclose (tfp);
|
||||
|
|
@ -495,50 +474,57 @@ Kinetic () { // Write Function here!
|
|||
// printf(" Total Kinetic Energy is %f\n",N*mvs*m/2.);
|
||||
return kin;
|
||||
}
|
||||
|
||||
double
|
||||
getF (double dist) {}
|
||||
|
||||
double
|
||||
getLocalPot (double dist) {}
|
||||
|
||||
// void
|
||||
// transposeMatrix (double mat[MAXPART][3], double matT[3][MAXPART]) {
|
||||
// for (int i = 0; i < 3; i++) {
|
||||
// for (int j = 0; j < MAXPART; j++) {
|
||||
// matT[i][j] = mat[j][i];
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
|
||||
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));
|
||||
double Pot = 0.;
|
||||
// double rT[3][MAXPART];
|
||||
// transposeMatrix (r, rT);
|
||||
|
||||
Pot = 0.;
|
||||
|
||||
for (i = 0; i < N - 1; i++) {
|
||||
for (j = i + 1; j < N; j++) {
|
||||
for (int i = 0; i < N - 1; i++) {
|
||||
for (int j = i + 1; j < N; j++) {
|
||||
double quot, rnorm, term1, term2;
|
||||
// CALCULATE POTENTIAL ENERGY
|
||||
dist = 0.;
|
||||
distTmp = 0;
|
||||
double dist = 0.;
|
||||
double posItoJ[3]; // position of i relative to j
|
||||
|
||||
for (int k = 0; k < 3; k++) {
|
||||
// POT
|
||||
distTmp = r[i][k] - r[j][k];
|
||||
double 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);
|
||||
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 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);
|
||||
|
||||
double distSqd = dist * dist;
|
||||
double rSqd_inv4 = 1.0 / (distSqd * distSqd);
|
||||
double rSqd_inv7 = rSqd_inv4 / (distSqd * dist);
|
||||
double f = 24.0 * (2.0 * rSqd_inv7 - rSqd_inv4);
|
||||
// from F = ma, where m = 1 in natural units!
|
||||
for (k = 0; k < 3; k++) {
|
||||
tmp = posItoJ[k] * f;
|
||||
for (int k = 0; k < 3; k++) {
|
||||
double tmp = posItoJ[k] * f;
|
||||
a[i][k] += tmp;
|
||||
a[j][k] -= tmp;
|
||||
}
|
||||
|
|
@ -599,8 +585,7 @@ computeAccelerations () {
|
|||
}
|
||||
// From derivative of Lennard-Jones with sigma and epsilon
|
||||
// set equal to 1 in natural units!
|
||||
double rSqd_inv7
|
||||
= 1.0 / (rSqd * rSqd * rSqd * rSqd * rSqd * rSqd * rSqd);
|
||||
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);
|
||||
|
||||
|
|
@ -648,16 +633,14 @@ VelocityVerlet (double dt, int iter, double *PE, FILE *fp) {
|
|||
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
|
||||
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
|
||||
v[i][j] *= -1.; //- elastic walls
|
||||
psum += 2 * m * fabs (v[i][j]) / dt; // contribution to pressure
|
||||
// from "right" walls
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue