/*
This code runs Monte Carlo simulations for Lennard Jones particles of a single type, in the canonical (NVT) ensemble.
The simulation cell is assumed to be cubic.
R. K. Lindsey (2023)
*/
// Project 1 - Lindsey Shereda
#include<iostream>
#include<iomanip>
#include<fstream>
#include<string>
#include<vector>
#include<cmath>
#include<random>
using namespace std;
double nav = 6.02*pow(10.0,23); // Avagadro's number
double cm2m = 1.0/100.0; // x cm * cm2m = m
double A2m = 1.0/pow(10.0,10); // x Angstrom * A to m = m
double kB = 1.380649*pow(10,-23.0); // Units: J⋅K^−1
struct xyz
{
double x;
double y;
double z;
// Added this operator so I can identify unequal vectors (to make identifying when atom positions are equal easier)
bool operator != (const xyz&other) const{
return x != other.x || y != other.y || z != other.z;
}
};
double get_dist(const xyz & a1, const xyz & a2, const xyz & boxdim, xyz & rij_vec)
{
double scalar_dist;
// Updates rij_vec with a1 & a2 information
rij_vec.x = a2.x-a1.x;
rij_vec.y = a2.y-a1.y;
rij_vec.z = a2.z-a1.z;
// Calculates the scalar distance
scalar_dist = sqrt(pow(rij_vec.x,2)+pow(rij_vec.y,2)+pow(rij_vec.z,2));
return scalar_dist;
}
double update_max_displacement(double fraction_accepted, double boxdim, double max_displacement)
{
// Calculates new max_displacement based on acceptance criteria
max_displacement -= round((max_displacement/boxdim)*boxdim);
// Returns the new max displacement
return max_displacement;
}
void write_frame(ofstream & trajstream, vector<xyz> & coords, string atmtyp, xyz & boxdim, int mcstep)
{
trajstream << "ITEM: TIMESTEP" << endl;
trajstream << mcstep << endl;
trajstream << "ITEM: NUMBER OF ATOMS" << endl;
trajstream << coords.size() << endl;
trajstream << "ITEM: BOX BOUNDS pp pp pp" << endl;
trajstream << "0 " << boxdim.x << endl;
trajstream << "0 " << boxdim.y << endl;
trajstream << "0 " << boxdim.z << endl;
trajstream << "ITEM: ATOMS id type element xu yu zu" << endl;
for (int i=0; i<coords.size(); i++)
trajstream << i << " 1 " << atmtyp << " " << coords[i].x << " "<< coords[i].y << " " << coords[i].z << endl;
}
double get_LJ_eij(double sigma, double epsilon, double rcut, double rij)
{
double LJ_energy;
if(rij <= rcut) {
LJ_energy = 4*epsilon*(pow((sigma / rij),12)-pow((sigma / rij),6));
}
return LJ_energy;
}
void get_LJ_fij(double sigma, double epsilon, double rcut, double rij, const xyz & rij_vec, xyz & fij)
{
double LJ_force;
if(rij <= rcut) {
fij.x = 48*epsilon*(((pow(sigma,12))/(pow(rij_vec.x,13)))-(0.5*((pow(sigma,6))/(pow(rij_vec.x,7)))));
fij.y = 48*epsilon*(((pow(sigma,12))/(pow(rij_vec.y,13)))-(0.5*((pow(sigma,6))/(pow(rij_vec.y,7)))));
fij.z = 48*epsilon*(((pow(sigma,12))/(pow(rij_vec.z,13)))-(0.5*((pow(sigma,6))/(pow(rij_vec.z,7)))));
}
}
void get_single_particle_LJ_contributions(double sigma, double epsilon, double rcut, const vector<xyz> & coords, int selected_atom, const xyz & selected_atom_coords, const xyz & boxdim, double & energy_selected, xyz & stress_selected)
{
energy_selected = 0;
stress_selected.x = 0;
stress_selected.y = 0;
stress_selected.z = 0;
xyz force;
force.x = 0;
force.y = 0;
force.z = 0;
static double rij;
static xyz rij_vec;
/* Write code to determine the contributions to the total system energy, forces, and stresses due to the selected atom.
Self interactions should not be included. */
for(int i=0; i<coords.size(); i++){
// If statement filters out self interactions
if (coords[i] != selected_atom_coords){
// Gets scalar distance and updates distance vector
double scalar_distance = get_dist(coords[i],selected_atom_coords,boxdim,rij_vec);
// Gets LJ energy and accumulates into energy selected
energy_selected += get_LJ_eij(sigma,epsilon,rcut,rij);
// Updates LJ force vector and adds to stress selected
get_LJ_fij(sigma,epsilon,rcut,rij,rij_vec,force);
stress_selected.x = stress_selected.x + force.x;
stress_selected.y = stress_selected.y + force.y;
stress_selected.z = stress_selected.z + force.z;
} else {
return;
}
}
}
int main(int argc, char* argv[])
{
///////////////////////////////////////////////////////////////////////////////////////////////
// Set user-defined variables (Read in from input file at some point)
///////////////////////////////////////////////////////////////////////////////////////////////
int seed = stoi(argv[1]);// Seed for random number generator - read from commandline
double redden = stod(argv[2]);// Reduced density - read from commandline
//////// Static variables
mt19937 mt(time(nullptr)); // Seed the (psuedo) random number
double sigma = 3.4; // Units: Angstrom
double epsilon = 120; // Units: K/k_B
double rcut = 4*sigma; // Model outer cutoff
int natoms = 500; // Number of atoms in the system
string atmtyp = "Ar"; // Atom type (atomic symbol)
double molmass = 39.948; // Mass of an particle of atmtyp (Atomic mass in the case of an argon atom)
double density = redden / pow(A2m,3.0) * pow(cm2m,3.0) / nav * molmass / pow(sigma,3.0); // System density; Units: g/cm^3
double numden; // System number density; Units: atoms/Ang^3 - will be calculated later on
double temp = 1.2*epsilon; // Temperature, in K
double nsteps = 5e6; // Number of MC steps
int iofrq = 2e4; // Frequency to output statistics and trajectory
int nequil = 2e6; // Equilibration period (chemical potential and heat capacity only collected after this many steps)
xyz boxdim; // Simulation box x, y, and z lengths - will be determined later on
vector<xyz> coords; // Coordinates for all atoms in our system; coords[0].x gives the x coordinate of the 0th atom - will be generated later on
//////// Variables used for file i/o
ofstream trajstream; // Trajectory file - uses the LAMMPS lammpstrj format
//////// Print information for user
cout << "# Number of atoms: " << natoms << endl;
cout << "# Atom type: " << atmtyp << endl;
cout << "# Molar Mass (g/mol): " << molmass << endl;
cout << "# Density (g/cm^3): " << density << endl;
cout << "# LJ sigma (Angstrom): " << sigma << endl;
cout << "# LJ epsilon/kB (K): " << epsilon << endl;
cout << "# LJ cutoff (Angstrom): " << rcut << endl;
cout << "# LJ cutoff (sigma): " << rcut/sigma << endl;
cout << "# Temperature (K): " << temp << endl;
cout << "# Number MC steps: " << nsteps << endl;
cout << "# Output frequency: " << iofrq << endl;
///////////////////////////////////////////////////////////////////////////////////////////////
// Generate the system coordinates
// Generates an initial configuation of atom on a 3D lattice with a user-specified number of
// atoms, at the user-specified density.
///////////////////////////////////////////////////////////////////////////////////////////////
// Determine the box length that will yield the user-specified density. Start by computing the target number density.
numden = density / molmass * nav / pow(cm2m,3.0) * pow(A2m,3.0); // Density in atoms per Ang^3
cout << "# Num. den. (atoms/Ang): " << numden << endl;
boxdim.x = pow(natoms/numden, 1.0/3.0);
boxdim.y = boxdim.x;
boxdim.z = boxdim.x;
cout << "# Box length (x): " << boxdim.x << endl;
cout << "# Box length (y): " << boxdim.y << endl;
cout << "# Box length (z): " << boxdim.z << endl;
// Compute the number of gridpoints to use in each direction (ceiling of the cube root of number of atoms).
// Set the spacing in each dimension based on the box length and the number of gridpoints+1 (to prevent overlap
// across the periodic boundary.
int ngridpoints = ceil(pow(natoms,1.0/3.0));
double spacing = boxdim.x / (ngridpoints+1);
cout << "# Init. spacing (Ang): " << spacing << endl;
xyz tmp_atom;
int added_atoms = 0;
for (int x=0; x<ngridpoints; x++)
{
for(int y=0; y<ngridpoints; y++)
{
for(int z=0; z<ngridpoints; z++)
{
if(added_atoms >= natoms)
continue;
tmp_atom.x = x * spacing;
tmp_atom.y = y * spacing;
tmp_atom.z = z * spacing;
coords.push_back(tmp_atom);
added_atoms++;
}
}
}
///////////////////////////////////////////////////////////////////////////////////////////////
// Print information on the initial configuration
///////////////////////////////////////////////////////////////////////////////////////////////
// Print this initial configuration to a file
trajstream.open("MC_traj.lammpstrj");
write_frame(trajstream, coords, atmtyp, boxdim, -1);
// Print some information for the user
cout << "# reduc. density (rho*): " << redden << endl;
cout << "# reduc. temp (T*): " << kB*temp / epsilon << endl;
///////////////////////////////////////////////////////////////////////////////////////////////
// Determine initial system energy and pressure
///////////////////////////////////////////////////////////////////////////////////////////////
double rij; // Scalar distance between atoms
xyz rij_vec; // Distance vector between atoms
double energy = 0; // Energy for the configuration
double widom_factor = 0; // The exponential term in the Widom chemical potential expression
double widom_trials = 0; // Number of Widom insertion attempts
double widom_avg = 0; // Average value of the wWidom factor
double Eavg = 0; // Average energy
double Esqavg = 0; // Square of average energy
xyz force; // Force on a given atom
xyz stensor; // System stress tensor (diagonal components only, i.e., xx, yy, zz)
stensor.x = 0;
stensor.y = 0;
stensor.z = 0;
for (int i=0; i<natoms; i++)
{
for (int j=i+1; j<natoms; j++)
{
// Determine atom pair's contirbution to total system energy
energy = get_LJ_eij(sigma,epsilon,rcut,rij);
return energy;
// Determines the atom pair's contribution to the total system pressure, updates stensor vector
get_LJ_fij(sigma,epsilon,rcut,rij,rij_vec,force);
stensor.x = force.x;
stensor.y = force.y;
stensor.z = force.z;
}
}
///////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////
// Begin the simulation
///////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////
double max_displacement = 0.5*boxdim.x; // Start trial displacements at one angstrom
int selected_atom;
double eold_selected; // Pre-trial-move energy of the selected atom
double enew_selected; // Post-trial-move energy of the selected atom
double delta_energy; // Difference between old and new (trial) energy
xyz sold_selected; // Pre-trial-move stress tensor diagonal of the selected atom
xyz snew_selected; // Post-trial-move stress tensor diagonal of the selected atom
xyz trial_displacement;
xyz trial_position;
int naccepted_moves = 0; // Number of accepted moves
double fraction_accepted; // Fraction of attempted moves that have been accepted
int nrunningav_moves = 0; // Move index for running averages (doesn't start until equilibration period has ended)
double pressure;
double Cv;
double stat_avgE = 0;
double stat_avgEsq = 0;
double stat_avgP = 0;
for (int i=0; i<nsteps; i++)
{
// Select a random particle. The syntax below shows how to use the random number generator. This generate a random integer between 0 and natoms-1
selected_atom = int(mt()*natoms);
// Determine contributions to the system's energy, force, and stress due to the selected atom
get_single_particle_LJ_contributions(sigma, epsilon, rcut, coords, selected_atom, coords[selected_atom], boxdim, eold_selected, sold_selected);
// Attempt to randomly move a particle - this is a multistep process
// 1. Generate the trial **displacement** in x, y, and z - the particle should be able to move in positive
// and negative directions, i.e., +/- the maximum displacement
trial_displacement.x = (1-2*rand())*max_displacement;
trial_displacement.y = (1-2*rand())*max_displacement;
trial_displacement.z = (1-2*rand())*max_displacement;
// 3. Apply PBC if the particle has moved outside the box
if(trial_displacement.x > max_displacement){
trial_displacement.x -= round((trial_displacement.x / boxdim.x)*boxdim.x);
}
if(trial_displacement.y > max_displacement){
trial_displacement.y -= round((trial_displacement.y / boxdim.y)*boxdim.y);
}
if(trial_displacement.z > max_displacement){
trial_displacement.z -= round((trial_displacement.z / boxdim.z)*boxdim.z);
}
// 2. Generate the trial **position** in x, y, and z based on the *corrected* displacement
trial_position.x = coords[selected_atom].x + trial_displacement.x;
trial_position.y = coords[selected_atom].y + trial_displacement.y;
trial_position.z = coords[selected_atom].z + trial_displacement.z;
// 4. Determine the energy contribution of that particle with the system **in it's trial position**
get_single_particle_LJ_contributions(sigma, epsilon, rcut, coords, selected_atom, trial_position, boxdim, enew_selected, snew_selected);
if (i >= nequil) // Only do Widom tests for the equilibrated portion of the simulation
{
// 5. Do a widom insertion test
double ewidom;
xyz swidom;
xyz widom_position;
// 5.a Generate another position for a new ghost atom
widom_position.x = (1-2*rand())*max_displacement;
widom_position.y = (1-2*rand())*max_displacement;
widom_position.z = (1-2*rand())*max_displacement;
// 5.b Calculate change in system energy due to insertion of new ghost atom
get_single_particle_LJ_contributions(sigma, epsilon, rcut, coords, -1, widom_position, boxdim, ewidom, swidom);
// 5.c Update the Widom factor
widom_factor = exp(-ewidom/temp);
widom_avg += widom_factor;
widom_trials++;
}
else
{
// Needed to avoid a divide-by-zero during equilibration phase when these values aren't collected
widom_avg = 1;
widom_trials = 1;
}
// 6. Accept or reject the move
delta_energy = enew_selected - eold_selected;
if (mt() < exp(-delta_energy/temp)) // Then accept
{
// Then the system energy has decreased **or** our random number is less than our probability to accept
// Update the system position
coords[selected_atom].x = trial_position.x;
coords[selected_atom].y = trial_position.y;
coords[selected_atom].z = trial_position.z;
// Update the system energy
energy += delta_energy;
// Update the system stress tensors
stensor.x += snew_selected.x - sold_selected.x;
stensor.y += snew_selected.y - sold_selected.y;
stensor.z += snew_selected.z - sold_selected.z;
// Update the number of accepted moves
naccepted_moves += naccepted_moves + 1;
}
// Update maximum diplacement to target a 50% acceptance rate
fraction_accepted = float(naccepted_moves)/float(i+1);
max_displacement = update_max_displacement(fraction_accepted, boxdim.x, max_displacement);
// print statistics if ineeded - don't forget to convert the stress tensor to pressure
// Compute instantaneous properties
// Gets Volume through density variable and converts cm^3 to m^3 in the equation
pressure = ((density / molmass*nav)/(1*10^-6))*kB*temp + 1/3*((density / molmass*nav)/(10*10^-6)) * (stensor.x + stensor.y + stensor.z);
Cv = 0;
if (i >= nequil) // Compute values for running averages, only using the equilibrated portion of the trajectory
{
stat_avgE += energy;
stat_avgEsq += energy*energy;
stat_avgP += pressure/epsilon*pow(sigma,3.0);
nrunningav_moves++;
double avgE = stat_avgE /float(nrunningav_moves);
double avgEsq = stat_avgEsq / float(nrunningav_moves);
Cv = (stat_avgEsq - pow(stat_avgE,2)) / (kB * pow(temp,2));
}
if ( (i+1) % iofrq == 0)
{
write_frame(trajstream, coords, atmtyp, boxdim, i);
cout << "Step: " << setw(10) << left << i;
cout << " NAcc: " << setw(10) << left << setprecision(3) << naccepted_moves;
cout << " fAcc: " << setw(10) << left << fixed << setprecision(3) << fraction_accepted;
cout << " Maxd: " << setw(10) << left << fixed << setprecision(5) << max_displacement;
cout << " E*/N: " << setw(10) << left << fixed << setprecision(5) << energy/natoms/epsilon;
cout << " P*: " << setw(10) << left << fixed << setprecision(5) << pressure/epsilon;
cout << " P*cold: " << setw(10) << left << fixed << setprecision(5) << ((density/molmass*nav)/(1*10^-6))*kB*temp + 1/3*((density/molmass*nav)/(10*10^-6)) * (stensor.x + stensor.y + stensor.z);
cout << " Mu_xs: " << setw(10) << left << fixed << setprecision(5) << log(widom_factor) / (1/kB*temp);
cout << " Cv*/N_xs: " << setw(15) << left << fixed << setprecision(5) << (stat_avgEsq - pow(stat_avgE,2)) / (kB * pow(temp,2));
cout << " E(kJ/mol): " << setw(10) << left << fixed << setprecision(3) << energy * 0.003814; // KJ/mol per K
cout << " P(bar): " << setw(10) << left << fixed << setprecision(3) << pressure * 0.003814 * 10.0e30 * 1000/(6.02*10.0e23)*1.0e-5; // KJ/mol/A^3 to bar
cout << endl;
}
}
trajstream.close();
stat_avgE /= float(nrunningav_moves);
stat_avgEsq /= float(nrunningav_moves);
Cv = (stat_avgEsq - pow(stat_avgE,2)) / (kB * pow(temp,2));
stat_avgE *= 1.0/natoms/epsilon;
cout << "# Computed average properties: " << endl;
cout << " # E*/N: " << setw(10) << left << fixed << setprecision(5) << stat_avgE << endl;
cout << " # P*: " << setw(10) << left << fixed << setprecision(5) << stat_avgP / float(nrunningav_moves) << endl;
cout << " # Cv*/N_xs: " << setw(15) << left << fixed << setprecision(5) << Cv << endl;
cout << " # Mu_xs: " << setw(10) << left << fixed << setprecision(5) << log(widom_factor) / (1/kB*temp) << endl;
} Write, Run & Share C++ code online using OneCompiler's C++ online compiler for free. It's one of the robust, feature-rich online compilers for C++ language, running on the latest version 17. Getting started with the OneCompiler's C++ compiler is simple and pretty fast. The editor shows sample boilerplate code when you choose language as C++ and start coding!
OneCompiler's C++ online compiler supports stdin and users can give inputs to programs using the STDIN textbox under the I/O tab. Following is a sample program which takes name as input and print your name with hello.
#include <iostream>
#include <string>
using namespace std;
int main()
{
string name;
cout << "Enter name:";
getline (cin, name);
cout << "Hello " << name;
return 0;
}
C++ is a widely used middle-level programming language.
When ever you want to perform a set of operations based on a condition If-Else is used.
if(conditional-expression) {
//code
}
else {
//code
}
You can also use if-else for nested Ifs and If-Else-If ladder when multiple conditions are to be performed on a single variable.
Switch is an alternative to If-Else-If ladder.
switch(conditional-expression){
case value1:
// code
break; // optional
case value2:
// code
break; // optional
......
default:
code to be executed when all the above cases are not matched;
}
For loop is used to iterate a set of statements based on a condition.
for(Initialization; Condition; Increment/decrement){
//code
}
While is also used to iterate a set of statements based on a condition. Usually while is preferred when number of iterations are not known in advance.
while (condition) {
// code
}
Do-while is also used to iterate a set of statements based on a condition. It is mostly used when you need to execute the statements atleast once.
do {
// code
} while (condition);
Function is a sub-routine which contains set of statements. Usually functions are written when multiple calls are required to same set of statements which increases re-usuability and modularity. Function gets run only when it is called.
return_type function_name(parameters);
function_name (parameters)
return_type function_name(parameters) {
// code
}