Skip to content

Instantly share code, notes, and snippets.

@carbolymer
Created October 17, 2012 21:57
Show Gist options
  • Save carbolymer/3908546 to your computer and use it in GitHub Desktop.
Save carbolymer/3908546 to your computer and use it in GitHub Desktop.
#include <iostream>
#include <math.h>
#include <fstream>
#include <vector>
#include <stdlib.h>
#include <time.h>
#include "argon.h"
using namespace std;
#define FILE_WRITE_STEP 200
/*
polozenia i pedy do jednego pliku:
x y z px py pz
kazdy krok czasowy oddzielic \n\n
w drugim pliku zapisywac czas, potencjal i energie kinetyczna
http://omega.if.pw.edu.pl/~mars/labKMS/
*/
int main()
{
state systemState;
double
a,
b0[3],
b1[3],
b2[3],
T,
m,
f,
L,
tau,
kb = 8.31e-3,
R = 0.38e-9,
epsilon = 1,
momentum,
*pairPotentials,
avTemperature = 0,
avHamiltonian = 0;
int
n,
n3,
n6,
w = 0, i,j,k,nparticle,So,Sd;
vect *pairDistances, *pairForces, totalMomentum, tMomentum;
fstream parametersFile("argon.in", fstream::in);
fstream crystalFile("avs.dat", fstream::out);
fstream stateFile("state.dat", fstream::out);
parametersFile >> n; // 1
parametersFile >> m; // 2
parametersFile >> epsilon; // 3
parametersFile >> R; // 4
parametersFile >> f; // 5
parametersFile >> L; // 6
parametersFile >> a; // 7
parametersFile >> T; // 8
parametersFile >> tau; // 9
parametersFile >> So; // 10
parametersFile >> Sd; // 11
cout << "Initial parameters:" << endl;
cout << "n:\t\t" << n << endl;
cout << "m:\t\t" << m << endl;
cout << "epsilon:\t" << epsilon << endl;
cout << "R:\t\t" << R << endl;
cout << "f:\t\t" << f << endl;
cout << "L:\t\t" << L << endl;
cout << "a:\t\t" << a << endl;
cout << "T_0:\t\t" << T << endl;
cout << "S_o:\t\t" << So << endl;
cout << "S_d:\t\t" << Sd << endl << endl;
srand(time(NULL));
n3 = n*n*n;
n6 = n3*n3;
momentum = sqrt(3.0*T*kb*m);
b0[0]=a;
b0[1]=0.0;
b0[2]=0.0;
b1[0]=a*0.5;
b1[1]=a*sqrt(3.)/2.;
b1[2]=0.0;
b2[0]=a*0.5;
b2[1]=a*sqrt(3.)/6.;
b2[2]=a*sqrt(2./3.);
pairDistances = new vect[n6]; // i+j*n3
pairForces = new vect[n6];
pairPotentials = new double[n6];
systemState.position = new vect[n3];
systemState.momentum = new vect[n3];
systemState.force = new vect[n3];
totalMomentum.x = 0;
totalMomentum.y = 0;
totalMomentum.z = 0;
// initial conditions
for(i=0; i<n; ++i)
for(j=0; j<n; ++j)
for(k=0; k<n; ++k)
{
systemState.position[w].x = (i-(n-1)*0.5)*b0[0]+(j-(n-1)*0.5)*b1[0]+(k-(n-1)*0.5)*b2[0];
systemState.position[w].y = (i-(n-1)*0.5)*b0[1]+(j-(n-1)*0.5)*b1[1]+(k-(n-1)*0.5)*b2[1];
systemState.position[w].z = (i-(n-1)*0.5)*b0[2]+(j-(n-1)*0.5)*b1[2]+(k-(n-1)*0.5)*b2[2];
// totalMomentum.x += systemState.momentum[w].x = momentum*(2*drand48()-1);
// totalMomentum.y += systemState.momentum[w].y = momentum*(2*drand48()-1);
// totalMomentum.z += systemState.momentum[w].z = momentum*(2*drand48()-1);
totalMomentum.x += systemState.momentum[w].x = momentum*(2*drand48()-1);
totalMomentum.y += systemState.momentum[w].y = momentum*(2*drand48()-1);
totalMomentum.z += systemState.momentum[w].z = momentum*(2*drand48()-1);
++w;
}
totalMomentum.x /= n3;
totalMomentum.y /= n3;
totalMomentum.z /= n3;
w = 0;
for(i=0; i<n; ++i)
for(j=0; j<n; ++j)
for(k=0; k<n; ++k)
{
// at the beginning system is not moving
systemState.momentum[w].x -= totalMomentum.x;
systemState.momentum[w].y -= totalMomentum.y;
systemState.momentum[w].z -= totalMomentum.z;
// crystalFile << systemState.position[w].x <<" "<< systemState.position[w].y <<" "<< systemState.position[w].z << " ";
// crystalFile << systemState.momentum[w].x <<" "<< systemState.momentum[w].y <<" "<< systemState.momentum[w].z << endl;
++w;
}
crystalFile << endl;
//
// **
//
for(int step = 0; step < So + Sd ; ++step)
{
systemState.pressure = 0;
systemState.potential = 0;
systemState.temperature = 0;
systemState.hamiltonian = 0;
// calculations for pairs
// i=j
for(i = 0; i < n3; ++i)
{
nparticle = i+i*n3;
pairDistances[nparticle].x = 0;
pairDistances[nparticle].y = 0;
pairDistances[nparticle].z = 0;
pairForces[nparticle].x = 0;
pairForces[nparticle].y = 0;
pairForces[nparticle].z = 0;
systemState.potential += pairPotentials[nparticle] = 0;
systemState.force[i].x = 0;
systemState.force[i].y = 0;
systemState.force[i].z = 0;
}
// i > j
for(i = 0; i < n3; ++i)
for(j = 0; j < i; ++j)
{
int npair = i + j*n3;
pairDistances[npair].x = systemState.position[i].x - systemState.position[j].x;
pairDistances[npair].y = systemState.position[i].y - systemState.position[j].y;
pairDistances[npair].z = systemState.position[i].z - systemState.position[j].z;
double distance = len(pairDistances[npair]);
double coefficient = 12.*epsilon*(pow(R/distance,12)-pow(R/distance,6))/distance/distance;
pairForces[npair].x = coefficient*pairDistances[npair].x;
pairForces[npair].y = coefficient*pairDistances[npair].y;
pairForces[npair].z = coefficient*pairDistances[npair].z;
pairPotentials[npair] = epsilon*(pow(R/distance,12)-2.*pow(R/distance,6));
systemState.force[i].x += pairForces[npair].x;
systemState.force[i].y += pairForces[npair].y;
systemState.force[i].z += pairForces[npair].z;
}
// i < j
for(i = 0; i < n3; ++i)
for(j = i+1; j < n3; ++j)
{
int npair = i + j*n3;
int npair2 = j+i*n3;
pairDistances[npair].x = -pairDistances[npair2].x;
pairDistances[npair].y = -pairDistances[npair2].y;
pairDistances[npair].z = -pairDistances[npair2].z;
pairForces[npair].x = -pairForces[npair2].x;
pairForces[npair].y = -pairForces[npair2].y;
pairForces[npair].z = -pairForces[npair2].z;
systemState.potential += pairPotentials[npair] = pairPotentials[npair2];
systemState.force[i].x += pairForces[npair].x;
systemState.force[i].y += pairForces[npair].y;
systemState.force[i].z += pairForces[npair].z;
}
// calculation for particles
for(i = 0; i < n3; ++i)
{
double distance = len(systemState.position[i]);
double coefficient;
if(distance >= L)
{
coefficient = f*(L-distance)/distance;
systemState.force[i].x += coefficient*systemState.position[i].x;
systemState.force[i].y += coefficient*systemState.position[i].y;
systemState.force[i].z += coefficient*systemState.position[i].z;
systemState.potential += f*pow(distance-L,2)/2.;
systemState.pressure += len(systemState.force[i]);
}
}
systemState.pressure /= 4.*M_PI*L*L;
// motion equations
for(i = 0; i < n3; ++i)
{
double energy = 0;
if(step > 0)
{
// 18c
systemState.momentum[i].x = tMomentum.x + systemState.force[i].x*tau/2;
systemState.momentum[i].y = tMomentum.y + systemState.force[i].y*tau/2;
systemState.momentum[i].z = tMomentum.z + systemState.force[i].z*tau/2;
}
energy = pow(len(systemState.momentum[i]),2)/2./m;
systemState.temperature += 2./3./n3/kb * energy;
systemState.hamiltonian += energy;
// 18a
tMomentum.x = systemState.momentum[i].x + systemState.force[i].x*tau/2;
tMomentum.y = systemState.momentum[i].y + systemState.force[i].y*tau/2;
tMomentum.z = systemState.momentum[i].z + systemState.force[i].z*tau/2;
// 18b
systemState.position[i].x = systemState.position[i].x + systemState.momentum[i].x*tau/m;
systemState.position[i].y = systemState.position[i].y + systemState.momentum[i].y*tau/m;
systemState.position[i].z = systemState.position[i].z + systemState.momentum[i].z*tau/m;
if(step % FILE_WRITE_STEP == 0 && step > So)
{
crystalFile << systemState.position[i].x <<" "<< systemState.position[i].y <<" "<< systemState.position[i].z << " ";
crystalFile << systemState.momentum[i].x <<" "<< systemState.momentum[i].y <<" "<< systemState.momentum[i].z << endl;
}
}
systemState.hamiltonian += systemState.potential;
if(step > So)
{
avTemperature += systemState.temperature/Sd;
avHamiltonian += systemState.hamiltonian/Sd;
}
if(step % FILE_WRITE_STEP == 0 && step > So)
{
stateFile << tau*step << " " << systemState.hamiltonian << " " << systemState.potential << " " << systemState.pressure << " " << systemState.temperature << endl;
crystalFile << endl << endl;
}
if(step % 100 == 0 && step != 0)
cout << "Progress: \t" << step << "\t/ " << So + Sd << endl;
}
cout << "Simulation complete. " << endl;
cout << "average H:\t" << avHamiltonian << endl;
cout << "average T:\t" << avTemperature << endl;
crystalFile.close();
parametersFile.close();
return 0;
}
#ifndef _ARGON_H_
#define _ARGON_H_
typedef struct
{
double x, y, z;
} vect;
typedef struct
{
vect
*position,
*momentum,
*force;
double
potential,
pressure,
hamiltonian,
temperature;
} state;
double len(vect position)
{
return sqrt(position.x*position.x+position.y*position.y+position.z*position.z);
}
#endif
5
39.948
1
0.38
1e4
2.3
0.38
100000
5e-4
100
30000
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment