I'm trying to catch values that are less than 0 or greater than 8
double L = 8;
if (rx[i]<0.0) { rx[i]+=L; b[i][0]--;}
else if (rx[i]>L) { rx[i]-=L; b[i][0]++;}
the problem is that I noticed that when rx[i] is 7.9996768475 it counts it as if its greater than 8 and it adds 1 to b[i][0].
I'm running an MD simulation and Im trying to save the unwrapped (unfolded) coordinates and so my box length (boundaries) are 0.0 and 8.0 but when Im printing my unfolded coordinates, I see that when values are close to 0 or close to 8 it counts them as if it crossed the boundary.
here is my code, go to the function integrate to see how i apply the boundary conditions:
//
// main.c
// MDsimulation
//
// Created by Fernanda Vargas on 3/23/21.
//
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
double v[512][3];
double x[512][3];
double xm[512][3];
double f[512][3];
double p[16667];
int b[512][3];
double realrand(){
int intRand = rand ();
return ((double) intRand / (double) RAND_MAX);
}
double halfrand(){
return (realrand() - 0.5);
}
void initialize(double * rx, double * ry, double * rz, double T, double dt){
int npart = 512;
double sumvx=0, sumvy=0, sumvz=0;
double sumv2=0;
double fs =0;
FILE* fp = fopen("MCsimulationLastConfiguration.txt","r");
for (int i=0; i<npart; i++){
/*get latest particle configuration*/
/*assign x y z values to p*/
fscanf(fp, "%lf %lf %lf",&rx[i],&ry[i],&rz[i]);
/*give random values to velocities*/
v[i][0] = halfrand();
v[i][1] = halfrand();
v[i][2] = halfrand();
/*velocity center of mass*/
sumvx += v[i][0];
sumvy += v[i][1];
sumvz += v[i][2];
/*kinetic energy*/
sumv2 +=v[i][0]*v[i][0]+v[i][1]*v[i][1]+v[i][2]*v[i][2];
}
fclose(fp);
/*velocity center of mass*/
sumvx = sumvx/npart;
sumvy = sumvy/npart;
sumvz = sumvz/npart;
/*mean-squared velocity*/
sumv2= sumv2/npart;
/*scale factor of the velocities*/
fs = sqrt(3*T/sumv2);
for (int i=0; i<npart; i++){
v[i][0] = (v[i][0]-sumvx)*fs;
v[i][1] = (v[i][1]-sumvy)*fs;
v[i][2] = (v[i][2]-sumvz)*fs;
xm[i][0] = rx[i] - v[i][0]*dt;
xm[i][1] = ry[i] - v[i][1]*dt;
xm[i][2] = rz[i] - v[i][2]*dt;
}
}
/*calculate force and energy*/
double forceenergy(double * rx, double * ry, double * rz, double Lhalf, double L, double sigma, double V, double rc){
double ene = 0;
double npart = 512;
double drmag, rstar, dx, dy, dz;
for (int i=0; i<npart; i++){
f[i][0] = 0;
f[i][1] = 0;
f[i][2] = 0;
}
for (int i=0;i<512;i++){
for(int j = i+1; j<512; j++){
/*choose 2 particles interacting*/
/*find dr calculating differences between x,y,z of the particles*/
dx = (rx[i]-rx[j]);
dy = (ry[i]-ry[j]);
dz = (rz[i]-rz[j]);
if (dx>Lhalf) dx-=L;
else if (dx<-Lhalf) dx+=L;
if (dy>Lhalf) dy-=L;
else if (dy<-Lhalf) dy+=L;
if (dz>Lhalf) dz-=L;
else if (dz<-Lhalf) dz+=L;
drmag=sqrt(dx*dx + dy*dy + dz*dz);
rstar = drmag/sigma;
if (rstar<2.5*sigma){
double rstar2 = (1/(rstar*rstar));
double rstar6 = rstar2*rstar2*rstar2;
double rstar12 = rstar6*rstar6;
ene += 4 *(rstar12-rstar6) - rc;
double force = (48/drmag)*(rstar12-0.5*rstar6);
f[i][0] += dx*force/rstar;
f[j][0] -= dx*force/rstar;
f[i][1] += dy*force/rstar;
f[j][1] -= dy*force/rstar;
f[i][2] += dz*force/rstar;
f[j][2] -= dz*force/rstar;
}
}
}
return ene;
}
double pressure(double * rx, double * ry, double * rz, double Lhalf, double L, double sigma, double V){
double drmag, rstar, dx, dy, dz;
double innerproduct = 0.0;
for (int i=0;i<512;i++){
for(int j = i+1; j<512; j++){
/*choose 2 particles interacting*/
/*find dr calculating differences between x,y,z of the particles*/
dx = (rx[i]-rx[j]);
dy = (ry[i]-ry[j]);
dz = (rz[i]-rz[j]);
if (dx>Lhalf) dx-=L;
else if (dx<-Lhalf) dx+=L;
if (dy>Lhalf) dy-=L;
else if (dy<-Lhalf) dy+=L;
if (dz>Lhalf) dz-=L;
else if (dz<-Lhalf) dz+=L;
drmag=sqrt(dx*dx + dy*dy + dz*dz);
rstar = drmag/sigma;
if (rstar<2.5*sigma){
double rstar2 = (1/(rstar*rstar));
double rstar6 = rstar2*rstar2*rstar2;
double rstar12 = rstar6*rstar6;
double force = (48/drmag)*(rstar12-0.5*rstar6);
innerproduct += (force*rstar);
}
else {
innerproduct += 0;
}
}
}
double fraction = (double) 1/3;
double pre = 1 + (fraction*innerproduct)/V;
return pre;
}
/*integrating the equations of motion*/
double integrate(double * rx, double * ry, double * rz, double dt, double L){
int npart = 512;
double sumvx=0, sumvy=0, sumvz=0;
double xx, xy, xz;
double sumv2=0;
for(int i=0; i <npart; i++){
double dxmn[3] = {rx[i] - xm[i][0], ry[i]-xm[i][1], rz[i]-xm[i][2]};
for (int d=0;d<3;d++){
if (dxmn[d]<-L/2) {dxmn[d]+=L;}
else if (dxmn[d]>L/2) {dxmn[d]-=L;}
}
xx = dxmn[0] + rx[i] + dt*dt*f[i][0];
xy = dxmn[1] + ry[i] + dt*dt*f[i][1];
xz = dxmn[2] + rz[i] + dt*dt*f[i][2];
v[i][0] = (xx - rx[i] + dxmn[0])/(2*dt);
v[i][1] = (xy - ry[i] + dxmn[1])/(2*dt);
v[i][2] = (xz - rz[i] + dxmn[2])/(2*dt);
/* Apply periodic boundary conditions and save a counter of how the coordinate has displaced (+ or -) and whats the accumulation*/
if (rx[i]<0.0) { rx[i]+=L; b[i][0]--;}
else if (rx[i]>L) { rx[i]-=L; b[i][0]++;}
if (ry[i]<0.0) { ry[i]+=L; b[i][1]--;}
else if (ry[i]>L) { ry[i]-=L; b[i][1]++;}
if (rz[i]<0.0) { rz[i]+=L; b[i][2]--;}
else if (rz[i]>L) { rz[i]-=L; b[i][2]++;}
/*velocity center of mass*/
sumvx += v[i][0];
sumvy += v[i][1];
sumvz += v[i][2];
/*kinetic energy*/
sumv2 +=v[i][0]*v[i][0]+v[i][1]*v[i][1]+v[i][2]*v[i][2]; /*check this*/
/*update positions previous time*/
xm[i][0] = rx[i];
xm[i][1] = ry[i];
xm[i][2] = rz[i];
//update positions current time
rx[i] = xx;
ry[i] = xy;
rz[i] = xz;
}
return sumv2;
}
int main(){
/*define variables (reduced units)*/
double * rx, * ry, * rz;
double T = 1; /*temperature*/
double L = 8.0;
double Lhalf = L/2;
double sigma = 1;
double V = L*L*L;
double rc = -0.0163168911;
int updatefrequency = 1;
double ene, pre,etot;
double dt = 0.001;
double npart = 512;
srand(100); /*strain of random numbers*/
/*initialize b*/
for (int i=0; i<512; i++){
b[i][0]=0;
b[i][1]=0;
b[i][2]=0;
}
/* Allocate the position arrays */
rx = (double*)malloc(512*sizeof(double));
ry = (double*)malloc(512*sizeof(double));
rz = (double*)malloc(512*sizeof(double));
/*initialize configuration*/
initialize(rx,ry,rz,T,dt);
int t=0;
FILE* fe = fopen("MDsimulationOutputEnergy.txt","w");
FILE* fepp = fopen("MDsimulationOutputEnergyperpart.txt","w");
FILE* fp = fopen("MDsimulationOutputPressure.txt","w");
FILE* fcon = fopen("MDsimulationConfigurations.txt","w");
FILE* fval = fopen("MDsimulationVelocities.txt","w");
FILE* funfol = fopen("MDsimulationUnfoldedCoordinates.txt","w");
FILE* fke = fopen("MDsimulationKE.txt","w");
/*50000*/
for (int i=0;i<50000;i++){
ene = forceenergy(rx,ry,rz,Lhalf,L,sigma,V,rc);
pre = pressure(rx,ry,rz,Lhalf,L,sigma,V);
double sumv2 = integrate(rx,ry,rz,dt,L);
double etot = (ene + 0.5*sumv2)/npart;
t += dt;
/*sample = sample(); what average am i calculating?*/
/*save energy and pressure in files*/
fprintf(fe,"%10.10f\n",ene);
fprintf(fp,"%10.10f\n",pre);
fprintf(fepp,"%10.10f\n",etot);
fprintf(fke,"%10.10f\n",sumv2);
if (i%updatefrequency==0){
/*save configurations to calculate RDF and velocity*/
for (int m=0;m<512;m++) {
fprintf(fcon,"%10.10f %10.10f %10.10f\n",rx[m], ry[m], rz[m]);
fprintf(fval,"%10.10f %10.10f %10.10f\n",v[m][0], v[m][1], v[m][2]);
fprintf(funfol,"%d %d %d %d %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f\n",m,b[m][0],b[m][1],b[m][2],xm[m][0],xm[m][1],xm[m][2],xm[m][0]+b[m][0]*L, xm[m][1]+b[m][1]*L, xm[m][2]+b[m][2]*L);
}
}
}
fclose(fe);
fclose(fp);
fclose(fcon);
fclose(fval);
fclose(fepp);
fclose(funfol);
printf("Ended Here"); /* An Indicator to make sure that the program has ended here
*/
return 0;
}
Aucun commentaire:
Enregistrer un commentaire