1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165

#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>
using namespace std;
vector<vector <double>> readin_pos(){
vector <vector <double>> pos;
ifstream inFile;
inFile.open("positions.txt");
if(inFile.fail()){ // Error test
cerr<< "Error opening file"<<endl;
exit(1);
}
vector <double> pos_xy;
double pos_x;
double pos_y;
while(!inFile.eof()){
inFile>> pos_x;
inFile>> pos_y;
pos_xy={pos_x, pos_y};
pos.push_back(pos_xy);
}
inFile.close();
pos.pop_back();
return pos;
}
vector<vector <double>> readin_vel(){
vector <vector <double>> vel;
ifstream inFile;
inFile.open("velocities.txt");
if(inFile.fail()){ // Error test
cerr<< "Error opening file"<<endl;
exit(1);
}
vector <double> vel_xy;
double vel_x;
double vel_y;
while(!inFile.eof()){
inFile>> vel_x;
inFile>> vel_y;
vel_xy={vel_x, vel_y};
vel.push_back(vel_xy);
}
inFile.close();
vel.pop_back();
return vel;
}
void calculate_forces(vector <vector<double>> F, vector <vector<double>> pos, double R, double L, double Epot){ //also calculate E_pot
double r;
double rx;
double ry;
double f_abs;
double fx;
double fy;
double total_fx;
double total_fy;
Epot=0;
//copy of posvector, used to consider images of particles in neighbouring boxes if necessary
vector <vector <double>> image_pos=pos;
//matrix storing the 2d forces between all particles
vector <double> **F_ij=new vector <double> *[pos.size()];
for(int i=0;i<pos.size();i++){
F_ij[i]=new vector <double>[pos.size()];
for(int j=0;j<pos.size();j++){ //initializing
F_ij[i][j]={0,0};
}
}
//filling the matrix
for(int i=0;i<pos.size();i++){
F_ij[i][i]={0,0}; //no force of particle on itself
for(int j=0; j<i;j++){
F_ij[i][j]={F_ij[j][i].at(0),F_ij[j][i].at(1)}; //Newton's 3rd law
}
for(int j=i+1; j<pos.size();j++){
rx=image_pos.at(i).at(0)image_pos.at(j).at(0); //distances between particles
ry=image_pos.at(i).at(1)image_pos.at(j).at(1);
if(rx>=L/2){
rx=L;
image_pos.at(j).at(0)+=L; //imaging particle to neighbouring box
}
else if(rx<=L/2){
rx+=L;
image_pos.at(j).at(0)=L;
}
if(ry>=L/2){
ry=L;
image_pos.at(j).at(1)+=L;
}
else if(ry<=L/2){
ry+=L;
image_pos.at(j).at(1)=L;
}
r=sqrt(pow(rx,2)+pow(ry,2));
if(r<=R){ //calculate forces
f_abs=48*pow(r,14)24*pow(r,8);
fx=f_abs*(image_pos.at(i).at(0)image_pos.at(j).at(0));
fy=f_abs*(image_pos.at(i).at(1)image_pos.at(j).at(1));
F_ij[i][j]={fx, fy};
Epot+=4*(pow(r,12)pow(r,6))+1;
}
}
}
//calculate total force vector F, holding the total forces for all particles
for(int i=0;i<pos.size();i++){
total_fx=0;
total_fy=0;
for(int j=0; j<pos.size(); j++){
total_fx+=F_ij[i][j].at(0);
total_fy+=F_ij[i][j].at(1);
}
F.at(i)={total_fx, total_fy};
}
delete[] F_ij;
}
void update_pos(vector <vector<double>> pos, vector <vector<double>> vel, vector <vector<double>> F, double dt){
for (int i=0;i<pos.size();i++){
pos.at(i).at(0)+=dt*vel.at(i).at(0)+(1./2)*pow(dt,2)*F.at(i).at(0);
pos.at(i).at(1)+=dt*vel.at(i).at(1)+(1./2)*pow(dt,2)*F.at(i).at(1);
}
}
void update_vel(vector <vector<double>> vel, vector <vector<double>> F, vector <vector <double>> F_old, double dt){
for (int i=0; i<vel.size();i++){
vel.at(i).at(0)=vel.at(i).at(0)+(1./2)*dt*(F.at(i).at(0)+F_old.at(i).at(0));
vel.at(i).at(1)=vel.at(i).at(1)+(1./2)*dt*(F.at(i).at(1)+F_old.at(i).at(1));
}
}
void calculate_Ekin(vector <vector <double>> vel, double Ekin){
Ekin=0;
for(int i=0;i<vel.size();i++){
Ekin+=(1./2)*(pow(vel.at(i).at(0),2)+pow(vel.at(i).at(1),2));
}
}
int main(){
double dt=0.0005;
double R=pow(2,1/6);
double L=14;
double Ekin;
double Epot;
vector <vector <double>> F_old; //storage for old forces
vector <vector <double>> pos=readin_pos();
vector <vector <double>> vel=readin_vel();
vector <vector <double>> F(pos.size(), vector <double>(2,0));
calculate_forces(F, pos, R, L, Epot);
calculate_Ekin(vel, Ekin);
//simulation step  velocityVerlet
for(int i=1; i<=200000;i++){
update_pos(pos,vel,F,dt);
F_old=F;
calculate_forces(F,pos,R,L, Epot);
update_vel(vel,F,F_old,dt);
calculate_Ekin(vel, Ekin);
cout <<i<<endl;
}
}
 