A2/MDS
2
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
MDS/src/MD-original.cpp

679 lines
21 KiB
C++
Raw Blame History

/*
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 <http://www.gnu.org/licenses/>.
Electronic Contact: foleyj10@wpunj.edu
Mail Contact: Prof. Jonathan Foley
Department of Chemistry, William Paterson University
300 Pompton Road
Wayne NJ 07470
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// Number of particles
int N;
// Lennard-Jones parameters in natural units!
double sigma = 1.;
double epsilon = 1.;
double m = 1.;
double kB = 1.;
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, 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 * 216;
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, 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();
PE = Potential();
// 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; i++) {
printf(" %6.3e %6.3e %6.3e\n",r[i][0],r[i][1],r[i][2]);
}
printf(" Printing initial velocities!\n");
for (i=0; i<N; i++) {
printf(" %6.3e %6.3e %6.3e\n",v[i][0],v[i][1],v[i][2]);
}
*/
}
// Function to calculate the averaged velocity squared
double MeanSquaredVelocity() {
double vx2 = 0;
double vy2 = 0;
double vz2 = 0;
double v2;
for (int i = 0; i < N; i++) {
vx2 = vx2 + v[i][0] * v[i][0];
vy2 = vy2 + v[i][1] * v[i][1];
vz2 = vz2 + v[i][2] * v[i][2];
}
v2 = (vx2 + vy2 + vz2) / N;
// printf(" Average of x-component of velocity squared is %f\n",v2);
return v2;
}
// Function to calculate the kinetic energy of the system
double Kinetic() { // Write Function here!
double v2, kin;
kin = 0.;
for (int i = 0; i < N; i++) {
v2 = 0.;
for (int j = 0; j < 3; j++) {
v2 += v[i][j] * v[i][j];
}
kin += m * v2 / 2.;
}
// printf(" Total Kinetic Energy is %f\n",N*mvs*m/2.);
return kin;
}
// Function to calculate the potential energy of the system
double Potential() {
double quot, r2, rnorm, term1, term2, Pot;
int i, j, k;
Pot = 0.;
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
if (j != i) {
r2 = 0.;
for (k = 0; k < 3; k++) {
r2 += (r[i][k] - r[j][k]) * (r[i][k] - r[j][k]);
}
rnorm = sqrt(r2);
quot = sigma / rnorm;
term1 = pow(quot, 12.);
term2 = pow(quot, 6.);
Pot += 4 * epsilon * (term1 - term2);
}
}
}
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() {
int i, j, k;
double f, rSqd;
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;
}
}
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;
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];
}
// 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));
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;
}
}
}
}
// 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++) {
r[i][j] += v[i][j] * dt + 0.5 * a[i][j] * dt * dt;
v[i][j] += 0.5 * a[i][j] * dt;
}
// 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();
// 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<N; i++) {
fprintf(fp,"%s",atype);
for (j=0; j<3; j++) {
fprintf(fp," %12.10e ",r[i][j]);
}
fprintf(fp,"\n");
}*/
// fprintf(fp,"\n \n");
return psum / (6 * L * L);
}
void initializeVelocities() {
int i, j;
for (i = 0; i < N; i++) {
for (j = 0; j < 3; j++) {
// Pull a number from a Gaussian Distribution
v[i][j] = gaussdist();
}
}
// Vcm = sum_i^N m*v_i/ sum_i^N M
// Compute center-of-mas velocity according to the formula above
double vCM[3] = {0, 0, 0};
for (i = 0; i < N; i++) {
for (j = 0; j < 3; j++) {
vCM[j] += m * v[i][j];
}
}
for (i = 0; i < 3; i++)
vCM[i] /= N * m;
// Subtract out the center-of-mass velocity from the
// velocity of each particle... effectively set the
// center of mass velocity to zero so that the system does
// not drift in space!
for (i = 0; i < N; i++) {
for (j = 0; j < 3; j++) {
v[i][j] -= vCM[j];
}
}
// Now we want to scale the average velocity of the system
// by a factor which is consistent with our initial temperature, Tinit
double vSqdSum, lambda;
vSqdSum = 0.;
for (i = 0; i < N; i++) {
for (j = 0; j < 3; j++) {
vSqdSum += v[i][j] * v[i][j];
}
}
lambda = sqrt(3 * (N - 1) * Tinit / vSqdSum);
for (i = 0; i < N; i++) {
for (j = 0; j < 3; j++) {
v[i][j] *= lambda;
}
}
}
// Numerical recipes Gaussian distribution number generator
double gaussdist() {
static bool available = false;
static double gset;
double fac, rsq, v1, v2;
if (!available) {
do {
v1 = 2.0 * rand() / double(RAND_MAX) - 1.0;
v2 = 2.0 * rand() / double(RAND_MAX) - 1.0;
rsq = v1 * v1 + v2 * v2;
} while (rsq >= 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;
}
}
Reference in a new issue View git blame Copy permalink