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

Fase 1 final

Browse source
This commit is contained in:
Afonso Franco 2023-10-22 16:06:13 +01:00
parent fbfb16c650
commit 896eb08a5c
Signed by: afonso
SSH key fingerprint: SHA256:JiuxZNdA5bRWXPMUJChI0AQ75yC+cXY4xM0IaVwEVys
80 changed files with 55 additions and 9341 deletions
Download patch file Download diff file Expand all files Collapse all files

100
src/MD.cpp
View file

@ -74,8 +74,11 @@ double Tinit; // 2;
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
@ -410,7 +413,6 @@ void initialize() {
int n, p, i, j, k;
double pos;
// Number of atoms in each direction
n = int(ceil(pow(N, 1.0 / 3)));
@ -589,34 +591,35 @@ double PotentialAndAcceleration(double dt) {
memset(a, 0, sizeof(a));
double Pot = 0.;
for (int i = 0; i < N-1; i++) {
double ri[3] = {r[i][0], r[i][1], r[i][2]};
for (int i = 0; i < N - 1; i++) {
double rI[3] = {r[i][0], r[i][1], r[i][2]};
for (int j = i + 1; j < N; j++) {
double quot, rnorm, term1, term2;
// CALCULATE POTENTIAL ENERGY
double dist = 0.;
double posItoJ[3]; // position of i relative to 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++) {
// POT
double distTmp = ri[k] - r[j][k];
dist += distTmp * distTmp;
// ACCEL
posItoJ[k] = distTmp;
double rijk = rI[k] - r[j][k];
rij[k] = rijk;
rSqd += rijk * rijk;
}
// From derivative of Lennard-Jones with sigma and epsilon
// set equal to 1 in natural units!
double distSqd = dist * dist * dist;
double rSqd_inv7 = distSqd * distSqd * dist;
double f = (48. - (24. * distSqd)) / rSqd_inv7;
// 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 = posItoJ[k] * f;
double tmp = rij[k] * f;
a[i][k] += tmp;
a[j][k] -= tmp;
}
quot = sigma * sigma / dist;
quot = sigma * sigma / rSqd;
term2 = quot * quot * quot;
Pot += epsilon_8 * term2 * (term2 - 1.);
}
@ -626,26 +629,23 @@ double PotentialAndAcceleration(double dt) {
// Function to calculate the potential energy of the system
double Potential() {
double quot, r2, rnorm, term1, term2, Pot;
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++) {
r2 = 0.;
rSqd = 0.;
double rI[3] = {r[i][0], r[i][1], r[i][2]};
for (int k = 0; k < 3; k++) {
double tmp = r[i][k] - r[j][k];
r2 += tmp * tmp;
double rijk = rI[k] - r[j][k];
rSqd += rijk * rijk;
}
quot = sigma / sqrt(r2);
term2 = quot * quot;
term2 = term2 * term2 * term2;
term1 = term2 * term2;
Pot += 8 * epsilon * (term1 - term2);
quot = sigma * sigma / rSqd;
term2 = quot * quot * quot;
Pot += epsilon_8 * term2 * (term2 - 1.);
}
}
return Pot;
@ -655,32 +655,32 @@ double Potential() {
// 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, tmp = 0.;
double rij[3]; // position of i relative to j
memset(a, 0, sizeof(a));
memset(a, 0, sizeof(a));
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 (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.;
// component-by-componenent position of i relative
// to j
for (k = 0; k < 3; k++) {
rij[k] = r[i][k] - r[j][k];
tmp = rij[k];
rSqd += tmp * tmp;
for (int k = 0; k < 3; k++) {
double rijk = rI[k] - r[j][k];
rij[k] = rijk;
rSqd += rijk * rijk;
}
// From derivative of Lennard-Jones with sigma and epsilon
// set equal to 1 in natural units!
double rSqd_inv7 = 1.0 / (rSqd * rSqd * rSqd * rSqd * rSqd * rSqd * rSqd);
double rSqd_inv4 = rSqd_inv7 * rSqd * rSqd * rSqd;
f = 24 * (2 * rSqd_inv7 - rSqd_inv4);
// from F = ma, where m = 1 in natural units!
for (k = 0; k < 3; k++) {
tmp = rij[k] * f;
// 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;
}