interact.c (Source)

#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "atoms.h"
#include "header.h"
#include "types.h"
#define getclosest(d,bo2,b)  ((d > bo2)?(d-b):((d < -bo2)?(d+b):d))
#define OUTPUTEVERY 1000
#define MAXNEIGHBORS 500
const double kB = 1.380658E-23;
Atom* allatoms;
int numatoms;
int** neighbors;
double bigX,bigY,bigZ;
double bXo2,bYo2,bZo2;
double mincutoff;
double maxdr;
double dt,hdt;
double desiredT;
double T;
double velMul;
int iterations;
void correctT()
{
int i;
T = 0;
for(i=0; i<numatoms; i++)
{
T += 0.5 * elements[allatoms[i].type].mass * (allatoms[i].vx*allatoms[i].vx + allatoms[i].vy*allatoms[i].vy + allatoms[i].vz*allatoms[i].vz);
}
T /= (double)numatoms;
T /= 1.5 * kB;
}
void makeneighbors()
{
int i,j,k;
double dx,dy,dz;
Atom* thisatom;
Atom* thatatom;
maxdr = 0.0;
printf("calculating neighbors at i=%d\n\n", iterations);
for(i=0; i<numatoms; i++)
{
thisatom = allatoms + i;
k = 0;
thisatom->ndx = thisatom->x;
thisatom->ndy = thisatom->y;
thisatom->ndz = thisatom->z;
for(j=i+1; j<numatoms && k<MAXNEIGHBORS; j++)
{
thatatom = allatoms + j;
dx = thisatom->x - thatatom->x;
dy = thisatom->y - thatatom->y;
dz = thisatom->z - thatatom->z;
dx = getclosest(dx,bXo2,bigX);
dy = getclosest(dy,bYo2,bigY);
dz = getclosest(dz,bZo2,bigZ);
if(dx<mincutoff&&dx>-mincutoff&&dy<mincutoff&&dy>-mincutoff&&dz<mincutoff&&dz>-mincutoff)
{
neighbors[i][k++] = j;
}
}
neighbors[i][k] = -1;
}
}
void interactions()
{
int i,j,k,l,m,n;
Atom* thisatom;
Atom* atom2;
Atom* atom3;
Atom* thatatom;
int atom2off;
int atom3off;
int atom4off;
double epsilon,sigma,thismass,thatmass,k0,x0,z;
double rmsq;
double dx,dy,dz;
double rsq,rmor2,rmor6,rmor12,r,rmx0;
double pot,mag;
double fx,fy,fz;
int bonded,doffset;
for(i=0; i<numatoms; i++)
{
thisatom = allatoms + i;
thismass = elements[thisatom->type].mass;
/*order 1 bonded interactions*/
for(j=0; j<NUMBONDS; j++)
{
if(thisatom->bonds[j] <= 0)
continue;
thatatom = thisatom + thisatom->bonds[j];
thatmass = elements[thatatom->type].mass;
switch(thisatom->type)
{
case H:
k0 = 483;
x0 = 1.12E-10;
break;
default:
switch(thatatom->type)
{
case H:
k0 = 483;
x0 = 1.12E-10;
break;
default:
k0 = 450;
x0 = 1.53E-10;
}
break;
}
dx = thisatom->x - thatatom->x;
dy = thisatom->y - thatatom->y;
dz = thisatom->z - thatatom->z;
dx = getclosest(dx,bXo2,bigX);
dy = getclosest(dy,bYo2,bigY);
dz = getclosest(dz,bZo2,bigZ);
rsq = dx*dx + dy*dy + dz*dz;
r = sqrt(rsq);
rmx0 = r-x0;
pot = 0.5 * k0 * rmx0 * rmx0;
mag = -k0 * rmx0 / r;
fx = dx*mag;
fy = dy*mag;
fz = dz*mag;
thisatom->newpotential += pot;
thisatom->ax += fx/thismass;
thisatom->ay += fy/thismass;
thisatom->az += fz/thismass;
thatatom->newpotential += pot;
thatatom->ax -= fx/thatmass;
thatatom->ay -= fy/thatmass;
thatatom->az -= fz/thatmass;
}
/*order 2 bonded interactions*/
for(j=0; j<NUMBONDS; j++)
{
atom2off = thisatom->bonds[j];
if(atom2off == 0)
continue;
atom2 = thisatom+atom2off;
for(k=0; k<NUMBONDS; k++)
{
atom3off = atom2off+atom2->bonds[k];
if(atom2->bonds[k] == 0 || atom3off <= 0)
continue;
thatatom = thisatom + atom3off;
thatmass = elements[thatatom->type].mass;
switch(thisatom->type)
{
case H:
switch(thatatom->type)
{
case H:
x0 = 1.83E-10;
break;
default:
x0 = 2.18E-10;
break;
}
break;
default:
switch(thatatom->type)
{
case H:
x0 = 2.18E-10;
break;
default:
x0 = 2.50E-10;
break;
}
break;
}
k0 = 10;
dx = thisatom->x - thatatom->x;
dy = thisatom->y - thatatom->y;
dz = thisatom->z - thatatom->z;
dx = getclosest(dx,bXo2,bigX);
dy = getclosest(dy,bYo2,bigY);
dz = getclosest(dz,bZo2,bigZ);
rsq = dx*dx + dy*dy + dz*dz;
r = sqrt(rsq);
rmx0 = r-x0;
pot = 0.5 * k0 * rmx0 * rmx0;
mag = -k0 * rmx0 / r;
fx = dx*mag;
fy = dy*mag;
fz = dz*mag;
thisatom->newpotential += pot;
thisatom->ax += fx/thismass;
thisatom->ay += fy/thismass;
thisatom->az += fz/thismass;
thatatom->newpotential += pot;
thatatom->ax -= fx/thatmass;
thatatom->ay -= fy/thatmass;
thatatom->az -= fz/thatmass;
}
}
/*order 3 bonded interactions*/
for(j=0; j<NUMBONDS; j++)
{
atom2off = thisatom->bonds[j];
if(atom2off == 0)
continue;
atom2 = thisatom + atom2off;
for(k=0; k<NUMBONDS; k++)
{
atom3off = atom2off + atom2->bonds[k];
if(atom2->bonds[k] == 0 || atom3off <= 0)
continue;
atom3 = thisatom + atom3off;
for(l=0; l<NUMBONDS; l++)
{
atom4off = atom3off + atom3->bonds[l];
if(atom3->bonds[k] == 0 || atom4off <= 0)
continue;
thatatom = thisatom + atom4off;
thatmass = elements[thatatom->type].mass;
switch(thisatom->type)
{
case H:
switch(thatatom->type)
{
case H:
x0 = 2.278E-10;
z = 2.392E-19;
break;
default:
x0 = 2.446E-10;
z = 2.234E-19;
break;
}
break;
default:
switch(thatatom->type)
{
case H:
x0 = 2.446E-10;
z = 2.234E-19;
break;
default:
x0 = 2.551E-10;
z = 2.660E-19;
break;
}
break;
}
dx = thisatom->x - thatatom->x;
dy = thisatom->y - thatatom->y;
dz = thisatom->z - thatatom->z;
dx = getclosest(dx,bXo2,bigX);
dy = getclosest(dy,bYo2,bigY);
dz = getclosest(dz,bZo2,bigZ);
rsq = dx*dx + dy*dy + dz*dz;
pot = z*x0*x0/rsq;
mag = 2.0*pot/rsq;
fx = dx*mag;
fy = dy*mag;
fz = dz*mag;
thisatom->newpotential += pot;
thisatom->ax += fx/thismass;
thisatom->ay += fy/thismass;
thisatom->az += fz/thismass;
thatatom->newpotential += pot;
thatatom->ax -= fx/thatmass;
thatatom->ay -= fy/thatmass;
thatatom->az -= fz/thatmass;
}
}
}
/*order 0, non-bonded (Van Der Waal's) interactions*/
for(j=0; j<MAXNEIGHBORS; j++)
{
if(neighbors[i][j] == -1)
break;
thatatom = allatoms + neighbors[i][j];
doffset = neighbors[i][j] - i;
bonded = 0;
for(n=0; n<NUMBONDS; n++)
{
if(doffset == thisatom->bonds[n])
bonded = 1;
for(m=0; m<NUMBONDS; m++)
{
if(doffset == (thisatom+thisatom->bonds[n])->bonds[m]+thisatom->bonds[n])
bonded = 1;
}
}
if(bonded)
continue;
thatmass = elements[thatatom->type].mass;
fx = fy = fz = 0.0;
if(thisatom->type == thatatom->type)
{
epsilon = elements[thisatom->type].epsilon;
sigma = 2.0 * elements[thisatom->type].radius;
}
else
{
epsilon = sqrt(elements[thisatom->type].epsilon*elements[thatatom->type].epsilon);
sigma = elements[thisatom->type].radius + elements[thatatom->type].radius;
}
rmsq = sigma*sigma;
dx = thisatom->x - thatatom->x;
dy = thisatom->y - thatatom->y;
dz = thisatom->z - thatatom->z;
dx = getclosest(dx,bXo2,bigX);
dy = getclosest(dy,bYo2,bigY);
dz = getclosest(dz,bZo2,bigZ);
rsq = dx*dx + dy*dy + dz*dz;
rmor2 = rmsq/rsq;
rmor6 = rmor2*rmor2*rmor2;
rmor12 = rmor6*rmor6;
pot = -4.0*epsilon*(rmor6-rmor12);
mag = -24.0*(epsilon/rsq)*(rmor6-2.0*rmor12);
fx = dx*mag;
fy = dy*mag;
fz = dz*mag;
thisatom->newpotential += pot;
thisatom->ax += fx/thismass;
thisatom->ay += fy/thismass;
thisatom->az += fz/thismass;
thatatom->newpotential += pot;
thatatom->ax -= fx/thatmass;
thatatom->ay -= fy/thatmass;
thatatom->az -= fz/thatmass;
}
}
}
void integrate()
{
int i;
double axhdt, ayhdt, azhdt;
double dx,dy,dz;
Atom* thisatom;
for(i=0; i<numatoms; i++)
{
thisatom = allatoms + i;
thisatom->potential = thisatom->newpotential;
/*temperature correction*/
thisatom->vx *= velMul;
thisatom->vy *= velMul;
thisatom->vz *= velMul;
/* velocity verlet integration */
axhdt = thisatom->ax * hdt;
ayhdt = thisatom->ay * hdt;
azhdt = thisatom->az * hdt;
thisatom->vx += axhdt;
thisatom->vy += ayhdt;
thisatom->vz += azhdt;
thisatom->x += thisatom->vx * dt + axhdt*dt;
thisatom->y += thisatom->vy * dt + ayhdt*dt;
thisatom->z += thisatom->vz * dt + azhdt*dt;
thisatom->vx += axhdt;
thisatom->vy += ayhdt;
thisatom->vz += azhdt;
if(thisatom->x > bigX) thisatom->x -= bigX;
else if(thisatom->x < 0.0) thisatom->x += bigX;
if(thisatom->y > bigY) thisatom->y -= bigY;
else if(thisatom->y < 0.0) thisatom->y += bigY;
if(thisatom->z > bigZ) thisatom->z -= bigZ;
else if(thisatom->z < 0.0) thisatom->z += bigZ;
dx = thisatom->x - (thisatom->ndx);
dy = thisatom->y - (thisatom->ndy);
dz = thisatom->z - (thisatom->ndz);
dx = getclosest(dx,bXo2,bigX);
dy = getclosest(dy,bYo2,bigY);
dz = getclosest(dz,bZo2,bigZ);
maxdr = max(max(maxdr,dx),max(dy,dz));;
thisatom->newpotential = 0.0;
thisatom->ax = 0.0;
thisatom->ay = 0.0;
thisatom->az = 0.0;
}
}
int main(int argc, char** argv)
{
int i,j;
int outevery = OUTPUTEVERY;
int keepgoing;
Headerinfo tempheader;
FILE* file;
int outnum;
char buffer[15];
int time0;
time0 = (int)time(NULL);
outnum = 0;
if(argc < 2)
{
printf("no infile specified\n\n");
return 0;
}
file = fopen(argv[1], "rb");
if(file == NULL)
{
printf("failed to open infile %s for reading\n\n", argv[1]);
return errno;
}
printf("using infile %s\n\n", argv[1]);
i = (int)fread(&tempheader, sizeof(Headerinfo), 1, file);
if(i < 1)
{
printf("error reading header from infile\n\n");
return 0;
}
iterations = tempheader.iterations;
if(iterations > 0)
keepgoing = 1;
else
keepgoing = 0;
mincutoff = tempheader.mincutoff;
dt = tempheader.dt;
bigX = tempheader.bigX;
bXo2 = bigX/2.0;
bigY = tempheader.bigY;
bYo2 = bigY/2.0;
bigZ = tempheader.bigZ;
bZo2 = bigZ/2.0;
desiredT = tempheader.desiredT;
T = desiredT;
numatoms = tempheader.numatoms;
allatoms = (Atom*)malloc(numatoms*sizeof(Atom));
if(allatoms == NULL)
{
printf("unable to allocate memory for allatoms\n\n");
return 0;
}
i = (int)fread(allatoms, sizeof(Atom), numatoms, file);
if(i < numatoms)
{
printf("error reading atoms from infile\n\n");
return 0;
}
fclose(file);
neighbors = (int**)malloc(numatoms*sizeof(int*));
if(neighbors == NULL)
{
printf("unable to allocate memory for neighbors\n\n");
free(allatoms);
return 0;
}
for(i=0; i<numatoms; i++)
{
neighbors[i] = (int*)malloc(MAXNEIGHBORS*sizeof(int));
if(neighbors[i] == NULL)
{
printf("unable to allocate memory for neighbors[%d]\n\n", i);
for(j=0; j<i; j++)
{
free(neighbors[i]);
}
free(neighbors);
free(allatoms);
return 0;
}
}
for(i=0; i<numatoms; i++)
{
for(j=0; j<MAXNEIGHBORS; j++)
{
neighbors[i][j] = -1;
}
}
makeneighbors();
hdt = 0.5 * dt;
while(keepgoing)
{
if(maxdr>0.5*mincutoff)
{
makeneighbors();
}
interactions();
correctT();
velMul = sqrt(desiredT/T);
integrate();
iterations--;
if(T>10000)
{
iterations = min(10, iterations);
outevery = 1;
}
if(iterations == outevery*(iterations/outevery))
{
file = NULL;
do
{
if(file != NULL) fclose(file);
outnum++;
sprintf(buffer,"o%07d.atm", outnum);
file = fopen(buffer, "r");
} while(file != NULL);
printf("%s\n", buffer);
printf("T=%f\n", T);
printf("i=%d\n", iterations);
printf("t=%d\n", (int)time(NULL) - time0);
printf("\n");
file = fopen(buffer, "wb");
tempheader.iterations = iterations;
tempheader.mincutoff = mincutoff;
tempheader.dt = dt;
tempheader.bigX = bigX;
tempheader.bigY = bigY;
tempheader.bigZ = bigZ;
tempheader.numatoms = numatoms;
tempheader.desiredT = desiredT;
tempheader.T = T;
fwrite(&tempheader, sizeof(Headerinfo), 1, file);
fwrite(allatoms, sizeof(Atom), numatoms, file);
fclose(file);
}
if(iterations > 0)
keepgoing = 1;
else
keepgoing = 0;
}
free(allatoms);
return 1;
}