A2/MDS
2
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

FINAL

Browse source
This commit is contained in:
Afonso Franco 2024-01-10 02:47:17 +00:00
parent b84aa107de
commit a8b0b317a4
Signed by: afonso
SSH key fingerprint: SHA256:JiuxZNdA5bRWXPMUJChI0AQ75yC+cXY4xM0IaVwEVys
10 changed files with 5638 additions and 70 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

33
Makefile
View file

@ -1,16 +1,33 @@
CFLAGS = -march=native -mtune=native -O2 -ftree-vectorize -mavx
CC = gcc
SRC = src/
CFLAGS = -march=native -mtune=native -mavx -O2 -ftree-vectorize -fopenmp
.DEFAULT_GOAL = MD.exe
.DEFAULT_GOAL = all
MD.exe: $(SRC)/MD2.cpp
$(CC) $(CFLAGS) $(SRC)MD2.cpp -lm -o MD.exe
all: MDseq.exe MDpar.exe
MDseq.exe: $(SRC)/MDseq.cpp
module load gcc/11.2.0;
$(CC) $(CFLAGS) $(SRC)MDseq.cpp -lm -o MDseq.exe
MDpar.exe: $(SRC)/MDpar.cpp
module load gcc/11.2.0;
$(CC) $(CFLAGS) $(SRC)MDpar.cpp -lm -fopenmp -o MDpar.exe
#FASE 3
MDvec.exe: $(SRC)/MDvec.cpp
module load gcc/11.2.0;
$(CC) $(CFLAGS) $(SRC)MDvec.cpp -lm -fopenmp -o MDvec.exe
clean:
rm ./MD.exe
rm ./MD*.exe
runseq:
./MDseq.exe < inputdata.txt
runpar:
OMP_NUM_THREADS=24 ./MDpar.exe < inputdata.txt
run:
./MD.exe < inputdata.txt
runfull: clean MD run
sbatch run.sh

0
relatorio.pdf → relatorio-fase1.pdf
View file

BIN
relatorio-fase2.pdf Normal file
View file

Binary file not shown.

3908
relatorio-fase3.pdf Normal file
View file

File diff suppressed because one or more lines are too long

9
run.sh Normal file
View file

@ -0,0 +1,9 @@
#!/bin/bash
#SBATCH --cpus-per-task=40
#SBATCH --time=00:10:00
#SBATCH --partition=cpar
#SBATCH --exclusive
export OMP_NUM_THREADS=28;
perf stat ./MDvec.exe < inputdata.txt

BIN
src/.DS_Store vendored
View file

Binary file not shown.

142
src/MD2.cpp → src/MDCUDA.cu
View file

@ -25,30 +25,11 @@
*/
#include <immintrin.h>
#include <math.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// 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)
#include <cuda.h>
#define ALIGNMENT 32
// Number of particles
@ -586,54 +567,93 @@ void transposeMatrix2(double matT[MAXPART][3], double mat[3][MAXPART]) {
// 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.;
__global__ void computeAccelerationsAndPotentialKernel(double *a_dev, double *r_dev, double *Pot_dev, double sigma, double epsilon_8) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
for (int k = 0; k < 3; k++) {
double rijk = rI[k] - r[j][k];
if (i < N - 1) {
double ai0 = 0.0, ai1 = 0.0, ai2 = 0.0;
double rI[3] = {r_dev[i * 3], r_dev[i * 3 + 1], r_dev[i * 3 + 2]};
double Pot_local = 0.0;
for (int j = i + 1; j < N; j++) {
double quot, term2;
double rij[3];
double rSqd = 0.0;
rij[k] = rijk;
for (int k = 0; k < 3; k++) {
double rijk = rI[k] - r_dev[j * 3 + k];
rij[k] = rijk;
rSqd += rijk * rijk;
}
rSqd += rijk * rijk;
}
double rSqd_3 = rSqd * rSqd * rSqd;
double rSqd_7 = rSqd_3 * rSqd_3 * rSqd;
double f = (48.0 - (24.0 * rSqd_3)) / rSqd_7;
// 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;
double tmp = rij[0] * f;
double tmp1 = rij[1] * f;
double tmp2 = rij[2] * f;
ai0 += tmp;
ai1 += tmp1;
ai2 += tmp2;
// Have to fix this data race
a_dev[j * 3] -= tmp;
a_dev[j * 3 + 1] -= tmp1;
a_dev[j * 3 + 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;
quot = sigma * sigma / rSqd;
term2 = quot * quot * quot;
Pot_local += term2 * (term2 - 1.0);
}
Pot_dev[0] = Pot_local;
// Have to fix this data race
a_dev[i * 3] += ai0;
a_dev[i * 3 + 1] += ai1;
a_dev[i * 3 + 2] += ai2;
}
}
double PotentialAndAcceleration(double dt) {
double *a_dev, *r_dev, *Pot_dev;
double Pot = 0.0;
// Allocate device memory
cudaMalloc((void**)&a_dev, N * 3 * sizeof(double));
cudaMalloc((void**)&r_dev, N * 3 * sizeof(double));
cudaMalloc((void**)&Pot_dev, N * sizeof(double));
// Copy data from host to device
cudaMemcpy(r_dev, r, N * 3 * sizeof(double), cudaMemcpyHostToDevice);
cudaMemset(Pot_dev, 0, N * sizeof(double), cudaMemcpyHostToDevice);
// Set a_dev to 0
cudaMemset(a_dev, 0, N * 3 * sizeof(double));
// Set grid and block sizes
int threadsPerBlock = 256;
int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;
// Launch the kernel
computeAccelerationsAndPotentialKernel<<<blocksPerGrid, threadsPerBlock>>>(a_dev, r_dev, Pot_dev, sigma, epsilon_8);
// Copy results back to host
cudaMemcpy(a, a_dev, N * 3 * sizeof(double), cudaMemcpyDeviceToHost);
cudaMemcpy(&Pot, Pot_dev, sizeof(double), cudaMemcpyDeviceToHost);
// Free device memory
cudaFree(a_dev);
cudaFree(r_dev);
// Accumulate the potential energy
for (int i = 0; i < N; i++) {
Pot += Pot_dev[i];
}
cudaFree(Pot_dev);
return Pot * epsilon_8;
}
// Function to calculate the potential energy of the system
double Potential() {
double quot, rSqd, rnorm, term1, term2, Pot;

741
src/MDpar.cpp Normal file
View file

@ -0,0 +1,741 @@
/*
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 <immintrin.h>
#include <math.h>
#include <omp.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 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; 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;
}
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];
}
}
}
void transposeMatrix2(double matT[MAXPART][3], double mat[3][MAXPART]) {
for (int i = 0; i < MAXPART; i++) {
for (int j = 0; j < 3; j++) {
matT[i][j] = mat[j][i];
}
}
}
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<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;
}
}

2
src/MD.cpp → src/MDseq.cpp
View file

@ -227,7 +227,7 @@ int main() {
scanf("%lf", &rho);
N = 10 * 216;
N = 10 * 500;
Vol = N / (rho * NA);
Vol /= VolFac;

873
src/MDvec.cpp Normal file
View file

@ -0,0 +1,873 @@
/*
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 <immintrin.h>
#include <math.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// Readability
typedef __m256d fourDoubles;
#define load_4_doubles(ptr) _mm256_loadu_pd(ptr)
#define store_4_doubles(ptr, val) _mm256_storeu_pd(ptr, val)
#define add_4_doubles(a, b) _mm256_add_pd(a, b)
#define subtract_4_doubles(a, b) _mm256_sub_pd(a, b)
#define multiply_4_doubles(a, b) _mm256_mul_pd(a, b)
#define divide_4_doubles(a, b) _mm256_div_pd(a, b)
#define create_4_doubles(a) _mm256_set1_pd(a)
#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 rx[MAXPART];
double ry[MAXPART];
double rz[MAXPART];
// Velocity
double v[MAXPART][3];
// Acceleration
double ax[MAXPART];
double ay[MAXPART];
double az[MAXPART];
// 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) {
rx[p] = (i + 0.5) * pos;
ry[p] = (j + 0.5) * pos;
rz[p] = (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;
}
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];
}
}
}
void transposeMatrix2(double matT[MAXPART][3], double mat[3][MAXPART]) {
for (int i = 0; i < MAXPART; i++) {
for (int j = 0; j < 3; j++) {
matT[i][j] = mat[j][i];
}
}
}
double PotentialAndAccelerationSIMD(double dt) {
// SAME AS OTHER FUCTION BUT WITH 4 POINTS AT A TIME
memset(ax, 0, sizeof(ax));
memset(ay, 0, sizeof(ay));
memset(az, 0, sizeof(az));
fourDoubles sigmaSqd = create_4_doubles(sigma * sigma);
double Pot = 0.;
#pragma omp parallel for reduction(+ : Pot, ax, ay, az)
for (int i = 0; i < N; i++) {
double aix = 0., aiy = 0., aiz = 0.;
fourDoubles rIx = create_4_doubles(rx[i]);
fourDoubles rIy = create_4_doubles(ry[i]);
fourDoubles rIz = create_4_doubles(rz[i]);
for (int j = i + 1; j < N; j += 4) {
// component-by-componenent position of the i point relative to the j
// point
fourDoubles ri[3];
// sum of squares of the components
fourDoubles rSqd = {0., 0., 0., 0.};
fourDoubles rJx = load_4_doubles(rx + j);
fourDoubles rIJx = subtract_4_doubles(rIx, rJx);
ri[0] = rIJx;
fourDoubles rJy = load_4_doubles(ry + j);
fourDoubles rIJy = subtract_4_doubles(rIy, rJy);
ri[1] = rIJy;
fourDoubles rJz = load_4_doubles(rz + j);
fourDoubles rIJz = subtract_4_doubles(rIz, rJz);
ri[2] = rIJz;
// Calculate the rSqd for each point
rSqd = add_4_doubles(multiply_4_doubles(rIJx, rIJx),
add_4_doubles(multiply_4_doubles(rIJy, rIJy),
multiply_4_doubles(rIJz, rIJz)));
// Here we remove the pow function and simplify the calculation
fourDoubles rSqd_3 =
multiply_4_doubles(rSqd, multiply_4_doubles(rSqd, rSqd));
fourDoubles rSqd_7 =
multiply_4_doubles(rSqd_3, multiply_4_doubles(rSqd_3, rSqd));
fourDoubles f = divide_4_doubles(
subtract_4_doubles(create_4_doubles(48.),
multiply_4_doubles(create_4_doubles(24.), rSqd_3)),
rSqd_7);
fourDoubles rIJFx = multiply_4_doubles(ri[0], f);
fourDoubles rIJFy = multiply_4_doubles(ri[1], f);
fourDoubles rIJFz = multiply_4_doubles(ri[2], f);
// Update the acceleration of the i particle
aix = aix + rIJFx[0] + rIJFx[1] + rIJFx[2] + rIJFx[3];
aiy = aiy + rIJFy[0] + rIJFy[1] + rIJFy[2] + rIJFy[3];
aiz = aiz + rIJFz[0] + rIJFz[1] + rIJFz[2] + rIJFz[3];
// Update the acceleration of the j particle
store_4_doubles(ax + j,
subtract_4_doubles(load_4_doubles(ax + j), rIJFx));
store_4_doubles(ay + j,
subtract_4_doubles(load_4_doubles(ay + j), rIJFy));
store_4_doubles(az + j,
subtract_4_doubles(load_4_doubles(az + j), rIJFz));
// Update the potential energy
fourDoubles quot = divide_4_doubles(sigmaSqd, rSqd);
fourDoubles quot3 =
multiply_4_doubles(quot, multiply_4_doubles(quot, quot));
fourDoubles Pots = multiply_4_doubles(
quot3, subtract_4_doubles(quot3, create_4_doubles(1.)));
Pot += Pots[0] + Pots[1] + Pots[2] + Pots[3];
}
// Store the acceleration of the i particle
ax[i] += aix;
ay[i] += aiy;
az[i] += aiz;
}
return Pot * epsilon_8;
}
// Function to calculate the potential energy of the system
double Potential() {
fourDoubles sigmaSqd = create_4_doubles(sigma * sigma);
double Pot = 0.;
fourDoubles Pots = {0., 0., 0., 0.};
#pragma omp parallel for reduction(+ : Pot)
for (int i = 0; i < N; i++) {
double aix = 0., aiy = 0., aiz = 0.;
fourDoubles rIx = create_4_doubles(rx[i]);
fourDoubles rIy = create_4_doubles(ry[i]);
fourDoubles rIz = create_4_doubles(rz[i]);
for (int j = i + 1; j < N; j += 4) {
// component-by-componenent position of the i point relative to the j
// point
fourDoubles ri[3];
// sum of squares of the components
fourDoubles rSqd = {0., 0., 0., 0.};
fourDoubles rJx = load_4_doubles(rx + j);
fourDoubles rIJx = subtract_4_doubles(rIx, rJx);
ri[0] = rIJx;
fourDoubles rJy = load_4_doubles(ry + j);
fourDoubles rIJy = subtract_4_doubles(rIy, rJy);
ri[1] = rIJy;
fourDoubles rJz = load_4_doubles(rz + j);
fourDoubles rIJz = subtract_4_doubles(rIz, rJz);
ri[2] = rIJz;
// Calculate the rSqd for each point
rSqd = add_4_doubles(multiply_4_doubles(rIJx, rIJx),
add_4_doubles(multiply_4_doubles(rIJy, rIJy),
multiply_4_doubles(rIJz, rIJz)));
// Update the potential energy
fourDoubles quot = divide_4_doubles(sigmaSqd, rSqd);
fourDoubles quot3 =
multiply_4_doubles(quot, multiply_4_doubles(quot, quot));
Pots = multiply_4_doubles(
quot3, subtract_4_doubles(quot3, create_4_doubles(1.)));
Pot += Pots[0] + Pots[1] + Pots[2] + Pots[3];
}
}
return Pot * epsilon_8;
}
// 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() {
// SAME AS OTHER FUCTION BUT WITH 4 POINTS AT A TIME
memset(ax, 0, sizeof(ax));
memset(ay, 0, sizeof(ay));
memset(az, 0, sizeof(az));
#pragma omp parallel for reduction(+ : ax, ay, az)
for (int i = 0; i < N; i++) {
double aix = 0., aiy = 0., aiz = 0.;
fourDoubles rIx = create_4_doubles(rx[i]);
fourDoubles rIy = create_4_doubles(ry[i]);
fourDoubles rIz = create_4_doubles(rz[i]);
for (int j = i + 1; j < N; j += 4) {
// component-by-componenent position of the i point relative to the j
// point
fourDoubles ri[3];
// sum of squares of the components
fourDoubles rSqd = {0., 0., 0., 0.};
fourDoubles rJx = load_4_doubles(rx + j);
fourDoubles rIJx = subtract_4_doubles(rIx, rJx);
ri[0] = rIJx;
fourDoubles rJy = load_4_doubles(ry + j);
fourDoubles rIJy = subtract_4_doubles(rIy, rJy);
ri[1] = rIJy;
fourDoubles rJz = load_4_doubles(rz + j);
fourDoubles rIJz = subtract_4_doubles(rIz, rJz);
ri[2] = rIJz;
// Calculate the rSqd for each point
rSqd = add_4_doubles(multiply_4_doubles(rIJx, rIJx),
add_4_doubles(multiply_4_doubles(rIJy, rIJy),
multiply_4_doubles(rIJz, rIJz)));
// Here we remove the pow function and simplify the calculation
fourDoubles rSqd_3 =
multiply_4_doubles(rSqd, multiply_4_doubles(rSqd, rSqd));
fourDoubles rSqd_7 =
multiply_4_doubles(rSqd_3, multiply_4_doubles(rSqd_3, rSqd));
fourDoubles f = divide_4_doubles(
subtract_4_doubles(create_4_doubles(48.),
multiply_4_doubles(create_4_doubles(24.), rSqd_3)),
rSqd_7);
fourDoubles rIJFx = multiply_4_doubles(ri[0], f);
fourDoubles rIJFy = multiply_4_doubles(ri[1], f);
fourDoubles rIJFz = multiply_4_doubles(ri[2], f);
// Update the acceleration of the i particle
aix = aix + rIJFx[0] + rIJFx[1] + rIJFx[2] + rIJFx[3];
aiy = aiy + rIJFy[0] + rIJFy[1] + rIJFy[2] + rIJFy[3];
aiz = aiz + rIJFz[0] + rIJFz[1] + rIJFz[2] + rIJFz[3];
// Update the acceleration of the j particle
store_4_doubles(ax + j,
subtract_4_doubles(load_4_doubles(ax + j), rIJFx));
store_4_doubles(ay + j,
subtract_4_doubles(load_4_doubles(ay + j), rIJFy));
store_4_doubles(az + j,
subtract_4_doubles(load_4_doubles(az + j), rIJFz));
}
// Store the acceleration of the i particle
ax[i] += aix;
ay[i] += aiy;
az[i] += aiz;
}
}
// 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++) {
double tmp = 0.5 * ax[i] * dt;
rx[i] += v[i][0] * dt + tmp * dt;
v[i][0] += tmp;
tmp = 0.5 * ay[i] * dt;
ry[i] += v[i][1] * dt + tmp * dt;
v[i][1] += tmp;
tmp = 0.5 * az[i] * dt;
rz[i] += v[i][2] * dt + tmp * dt;
v[i][2] += 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 = PotentialAndAccelerationSIMD(dt);
// Update velocity with updated acceleration
for (i = 0; i < N; i++) {
v[i][0] += 0.5 * ax[i] * dt;
v[i][1] += 0.5 * ay[i] * dt;
v[i][2] += 0.5 * az[i] * dt;
}
// Elastic walls
for (i = 0; i < N; i++) {
for (j = 0; j < 3; j++) {
double *r;
if (j == 0)
r = rx;
else if (j == 1)
r = ry;
else if (j == 2)
r = rz;
if (r[i] < 0.) {
v[i][j] *= -1.; //- elastic walls
psum += 2 * m * fabs(v[i][j]) / dt; // contribution to pressure
// from "left" walls
}
if (r[i] >= 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;
}
}