*
*/
-#include "moldyn.h"
-
+#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
#include <math.h>
+#include "moldyn.h"
+
#include "math/math.h"
#include "init/init.h"
#include "random/random.h"
#include "visual/visual.h"
+#include "list/list.h"
+
+int moldyn_usage(char **argv) {
+
+ printf("\n%s usage:\n\n",argv[0]);
+ printf("--- general options ---\n");
+ printf("-E <steps> <file> (log total energy)\n");
+ printf("-M <steps> <file> (log total momentum)\n");
+ printf("-D <steps> <file> (dump total information)\n");
+ printf("-S <steps> <filebase> (single save file)\n");
+ printf("-V <steps> <filebase> (rasmol file)\n");
+ printf("--- physics options ---\n");
+ printf("-T <temperature> [K] (%f)\n",MOLDYN_TEMP);
+ printf("-t <timestep tau> [s] (%.15f)\n",MOLDYN_TAU);
+ printf("-R <runs> (%d)\n",MOLDYN_RUNS);
+ printf(" -- integration algo --\n");
+ printf(" -I <number> (%d)\n",MOLDYN_INTEGRATE_DEFAULT);
+ printf(" 0: velocity verlet\n");
+ printf(" -- potential --\n");
+ printf(" -P <number> <param1 param2 ...>\n");
+ printf(" 0: harmonic oscillator\n");
+ printf(" param1: spring constant\n");
+ printf(" param2: equilibrium distance\n");
+ printf(" 1: lennard jones\n");
+ printf(" param1: epsilon\n");
+ printf(" param2: sigma\n");
+ printf("\n");
+
+ return 0;
+}
+
+int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
+
+ int i;
+ t_ho_params hop;
+ t_lj_params ljp;
+ double s,e;
+
+ memset(moldyn,0,sizeof(t_moldyn));
+
+ /* default values */
+ moldyn->t=MOLDYN_TEMP;
+ moldyn->tau=MOLDYN_TAU;
+ moldyn->time_steps=MOLDYN_RUNS;
+ moldyn->integrate=velocity_verlet;
+ moldyn->potential_force_function=lennard_jones;
+
+ /* parse argv */
+ for(i=1;i<argc;i++) {
+ if(argv[i][0]=='-') {
+ switch(argv[i][1]){
+ case 'E':
+ moldyn->ewrite=atoi(argv[++i]);
+ strncpy(moldyn->efb,argv[++i],64);
+ break;
+ case 'M':
+ moldyn->mwrite=atoi(argv[++i]);
+ strncpy(moldyn->mfb,argv[++i],64);
+ break;
+ case 'D':
+ moldyn->dwrite=atoi(argv[++i]);
+ strncpy(moldyn->dfb,argv[++i],64);
+ break;
+ case 'S':
+ moldyn->swrite=atoi(argv[++i]);
+ strncpy(moldyn->sfb,argv[++i],64);
+ break;
+ case 'V':
+ moldyn->vwrite=atoi(argv[++i]);
+ strncpy(moldyn->vfb,argv[++i],64);
+ break;
+ case 'T':
+ moldyn->t=atof(argv[++i]);
+ break;
+ case 't':
+ moldyn->tau=atof(argv[++i]);
+ break;
+ case 'R':
+ moldyn->time_steps=atoi(argv[++i]);
+ break;
+ case 'I':
+ /* integration algorithm */
+ switch(atoi(argv[++i])) {
+ case MOLDYN_INTEGRATE_VERLET:
+ moldyn->integrate=velocity_verlet;
+ break;
+ default:
+ printf("unknown integration algo %s\n",argv[i]);
+ moldyn_usage(argv);
+ return -1;
+ }
+
+ case 'P':
+ /* potential + params */
+ switch(atoi(argv[++i])) {
+ case MOLDYN_POTENTIAL_HO:
+ hop.spring_constant=atof(argv[++i]);
+ hop.equilibrium_distance=atof(argv[++i]);
+ moldyn->pot_params=malloc(sizeof(t_ho_params));
+ memcpy(moldyn->pot_params,&hop,sizeof(t_ho_params));
+ moldyn->potential_force_function=harmonic_oscillator;
+ break;
+ case MOLDYN_POTENTIAL_LJ:
+ e=atof(argv[++i]);
+ s=atof(argv[++i]);
+ ljp.epsilon4=4*e;
+ ljp.sigma6=s*s*s*s*s*s;
+ ljp.sigma12=ljp.sigma6*ljp.sigma6;
+ moldyn->pot_params=malloc(sizeof(t_lj_params));
+ memcpy(moldyn->pot_params,&ljp,sizeof(t_lj_params));
+ moldyn->potential_force_function=lennard_jones;
+ break;
+ default:
+ printf("unknown potential %s\n",argv[i]);
+ moldyn_usage(argv);
+ return -1;
+ }
+
+ default:
+ printf("unknown option %s\n",argv[i]);
+ moldyn_usage(argv);
+ return -1;
+ }
+ } else {
+ moldyn_usage(argv);
+ return -1;
+ }
+ }
+
+ return 0;
+}
+
+int moldyn_log_init(t_moldyn *moldyn,void *v) {
+
+ moldyn->lvstat=0;
+ t_visual *vis;
+
+ vis=v;
+
+ if(moldyn->ewrite) {
+ moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC);
+ if(moldyn->efd<0) {
+ perror("[moldyn] efd open");
+ return moldyn->efd;
+ }
+ dprintf(moldyn->efd,"# moldyn total energy logfile\n");
+ moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E;
+ }
+
+ if(moldyn->mwrite) {
+ moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC);
+ if(moldyn->mfd<0) {
+ perror("[moldyn] mfd open");
+ return moldyn->mfd;
+ }
+ dprintf(moldyn->mfd,"# moldyn total momentum logfile\n");
+ moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M;
+ }
+
+ if(moldyn->swrite)
+ moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
+ if(moldyn->dwrite) {
+ moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
+ if(moldyn->dfd<0) {
+ perror("[moldyn] dfd open");
+ return moldyn->dfd;
+ }
+ write(moldyn->dfd,moldyn,sizeof(t_moldyn));
+ moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
+ }
+
+ if((moldyn->vwrite)&&(vis)) {
+ moldyn->visual=vis;
+ visual_init(vis,moldyn->vfb);
+ moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
+ }
+
+ moldyn->lvstat|=MOLDYN_LVSTAT_INITIALIZED;
+
+ return 0;
+}
+
+int moldyn_shutdown(t_moldyn *moldyn) {
+
+ if(moldyn->efd) close(moldyn->efd);
+ if(moldyn->mfd) close(moldyn->efd);
+ if(moldyn->dfd) close(moldyn->efd);
+ if(moldyn->visual) visual_tini(moldyn->visual);
+
+ return 0;
+}
int create_lattice(unsigned char type,int element,double mass,double lc,
int a,int b,int c,t_atom **atom) {
return 0;
}
-int thermal_init(t_atom *atom,t_random *random,int count,double t) {
+int thermal_init(t_moldyn *moldyn,t_random *random) {
/*
* - gaussian distribution of velocities
int i;
double v,sigma;
t_3dvec p_total,delta;
+ t_atom *atom;
+
+ atom=moldyn->atom;
/* gaussian distribution of velocities */
v3_zero(&p_total);
- for(i=0;i<count;i++) {
- sigma=sqrt(2.0*K_BOLTZMANN*t/atom[i].mass);
+ for(i=0;i<moldyn->count;i++) {
+ sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
/* x direction */
v=sigma*rand_get_gauss(random);
atom[i].v.x=v;
}
/* zero total momentum */
- v3_scale(&p_total,&p_total,1.0/count);
- for(i=0;i<count;i++) {
+ v3_scale(&p_total,&p_total,1.0/moldyn->count);
+ for(i=0;i<moldyn->count;i++) {
v3_scale(&delta,&p_total,1.0/atom[i].mass);
v3_sub(&(atom[i].v),&(atom[i].v),&delta);
}
/* velocity scaling */
- scale_velocity(atom,count,t);
+ scale_velocity(moldyn);
return 0;
}
-int scale_velocity(t_atom *atom,int count,double t) {
+int scale_velocity(t_moldyn *moldyn) {
int i;
double e,c;
+ t_atom *atom;
+
+ atom=moldyn->atom;
/*
* - velocity scaling (E = 3/2 N k T), E: kinetic energy
*/
e=0.0;
- for(i=0;i<count;i++)
+ for(i=0;i<moldyn->count;i++)
e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
- c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*t));
- for(i=0;i<count;i++)
+ c=sqrt((2.0*e)/(3.0*moldyn->count*K_BOLTZMANN*moldyn->t));
+ for(i=0;i<moldyn->count;i++)
v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
return 0;
double get_e_pot(t_moldyn *moldyn) {
- return(moldyn->potential(moldyn));
+ return moldyn->energy;
}
double get_total_energy(t_moldyn *moldyn) {
return tau;
}
+/*
+ * numerical tricks
+ */
+
+/* linked list / cell method */
+
+int link_cell_init(t_moldyn *moldyn) {
+
+ t_linkcell *lc;
+
+ lc=&(moldyn->lc);
+
+ /* partitioning the md cell */
+ lc->nx=moldyn->dim.x/moldyn->cutoff;
+ lc->x=moldyn->dim.x/lc->nx;
+ lc->ny=moldyn->dim.y/moldyn->cutoff;
+ lc->y=moldyn->dim.y/lc->ny;
+ lc->nz=moldyn->dim.z/moldyn->cutoff;
+ lc->z=moldyn->dim.z/lc->nz;
+
+ lc->subcell=malloc(lc->nx*lc->ny*lc->nz*sizeof(t_list));
+
+ link_cell_update(moldyn);
+
+ return 0;
+}
+
+int link_cell_update(t_moldyn *moldyn) {
+
+ int count,i,j,k;
+ int nx,ny,nz;
+ t_atom *atom;
+ t_linkcell *lc;
+
+ atom=moldyn->atom;
+ lc=&(moldyn->lc);
+
+ nx=lc->nx;
+ ny=lc->ny;
+ nz=lc->nz;
+
+ for(i=0;i<nx*ny*nz;i++)
+ list_destroy(&(moldyn->lc.subcell[i]));
+
+ for(count=0;count<moldyn->count;count++) {
+ i=atom[count].r.x/lc->x;
+ j=atom[count].r.y/lc->y;
+ k=atom[count].r.z/lc->z;
+ list_add_immediate_ptr(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]),
+ &(atom[count]));
+ }
+
+ return 0;
+}
+
+int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
+
+ t_linkcell *lc;
+ int a;
+ int count1,count2;
+ int ci,cj,ck;
+ int nx,ny,nz;
+ int x,y,z;
+ unsigned char bx,by,bz;
+
+ lc=&(moldyn->lc);
+ nx=lc->nx;
+ ny=lc->ny;
+ nx=lc->nz;
+ count1=1;
+ count2=27;
+ a=nx*ny;
+
+ cell[0]=lc->subcell[i+j*nx+k*a];
+ for(ci=-1;ci<=1;ci++) {
+ bx=0;
+ x=i+ci;
+ if((x<0)||(x>=nx)) {
+ x=(x+nx)%nx;
+ bx=1;
+ }
+ for(cj=-1;cj<=1;cj++) {
+ by=0;
+ y=j+cj;
+ if((y<0)||(y>=ny)) {
+ y=(y+ny)%ny;
+ by=1;
+ }
+ for(ck=-1;ck<=1;ck++) {
+ bz=0;
+ z=k+ck;
+ if((z<0)||(z>=nz)) {
+ z=(z+nz)%nz;
+ bz=1;
+ }
+ if(!(x|y|z)) continue;
+ if(bx|by|bz) {
+ cell[--count2]=lc->subcell[x+y*nx+z*a];
+ }
+ else {
+ cell[count1++]=lc->subcell[x+y*nx+z*a];
+ }
+ }
+ }
+ }
+
+ return count2;
+}
+
+int link_cell_shutdown(t_moldyn *moldyn) {
+
+ int i;
+ t_linkcell *lc;
+
+ lc=&(moldyn->lc);
+
+ for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
+ list_shutdown(&(moldyn->lc.subcell[i]));
+
+ return 0;
+}
/*
*
int moldyn_integrate(t_moldyn *moldyn) {
int i;
- int write;
+ unsigned int e,m,s,d,v;
+ t_3dvec p;
+
+ int fd;
+ char fb[128];
+
+ /* logging & visualization */
+ e=moldyn->ewrite;
+ m=moldyn->mwrite;
+ s=moldyn->swrite;
+ d=moldyn->dwrite;
+ v=moldyn->vwrite;
+
+ if(!(moldyn->lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
+ printf("[moldyn] warning, lv system not initialized\n");
+ return -1;
+ }
+
+ /* sqaure of some variables */
+ moldyn->tau_square=moldyn->tau*moldyn->tau;
+ moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
- write=moldyn->write;
+ /* create the neighbour list */
+ link_cell_update(moldyn);
/* calculate initial forces */
- moldyn->force(moldyn);
+ moldyn->potential_force_function(moldyn);
for(i=0;i<moldyn->time_steps;i++) {
+ /* show runs */
+ printf(".");
+
+ /* neighbour list update */
+ link_cell_update(moldyn);
+
/* integration step */
moldyn->integrate(moldyn);
- /* check for visualiziation */
- if(!(i%write)) {
- visual_atoms(moldyn->visual,i*moldyn->tau,
- moldyn->atom,moldyn->count);
- printf("finished %d / %d\n",i,moldyn->time_steps);
+ /* check for log & visualization */
+ if(e) {
+ if(!(i%e))
+ dprintf(moldyn->efd,
+ "%.15f %.45f\n",i*moldyn->tau,
+ get_total_energy(moldyn));
+ }
+ if(m) {
+ if(!(i%m)) {
+ p=get_total_p(moldyn->atom,moldyn->count);
+ dprintf(moldyn->mfd,
+ "%.15f %.45f\n",i*moldyn->tau,
+ v3_norm(&p));
+ }
+ }
+ if(s) {
+ if(!(i%s)) {
+ snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
+ moldyn->t,i*moldyn->tau);
+ fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
+ if(fd<0) perror("[moldyn] save fd open");
+ else {
+ write(fd,moldyn,sizeof(t_moldyn));
+ write(fd,moldyn->atom,
+ moldyn->count*sizeof(t_atom));
+ }
+ }
+ }
+ if(d) {
+ if(!(i%d))
+ write(moldyn->dfd,moldyn->atom,
+ moldyn->count*sizeof(t_atom));
+
+ }
+ if(v) {
+ if(!(i%v))
+ visual_atoms(moldyn->visual,i*moldyn->tau,
+ moldyn->atom,moldyn->count);
}
}
atom=moldyn->atom;
count=moldyn->count;
tau=moldyn->tau;
-
- tau_square=tau*tau;
+ tau_square=moldyn->tau_square;
for(i=0;i<count;i++) {
/* new positions */
}
/* forces depending on chosen potential */
- moldyn->force(moldyn);
+ moldyn->potential_force_function(moldyn);
for(i=0;i<count;i++) {
/* again velocities */
/* harmonic oscillator potential and force */
-double potential_harmonic_oscillator(t_moldyn *moldyn) {
+int harmonic_oscillator(t_moldyn *moldyn) {
t_ho_params *params;
- t_atom *atom;
- int i,j;
+ t_atom *atom,*btom;
+ t_linkcell *lc;
+ t_list *this,neighbour[27];
+ int i,j,c;
int count;
- t_3dvec distance;
+ t_3dvec force,distance;
double d,u;
double sc,equi_dist;
+ int ni,nj,nk;
params=moldyn->pot_params;
atom=moldyn->atom;
+ lc=&(moldyn->lc);
sc=params->spring_constant;
equi_dist=params->equilibrium_distance;
count=moldyn->count;
u=0.0;
for(i=0;i<count;i++) {
- for(j=0;j<i;j++) {
- v3_sub(&distance,&(atom[i].r),&(atom[j].r));
+ /* determine cell + neighbours */
+ ni=atom[i].r.x/lc->x;
+ nj=atom[i].r.y/lc->y;
+ nk=atom[i].r.z/lc->z;
+ c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
+
+ /* processing cell of atom i */
+ this=&(neighbour[0]);
+ list_reset(this); /* list has 1 element at minimum */
+ do {
+ btom=this->current->data;
+ if(btom==&(atom[i]))
+ continue;
+ v3_sub(&distance,&(atom[i].r),&(btom->r));
d=v3_norm(&distance);
u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
- }
- }
-
- return u;
-}
-
-int force_harmonic_oscillator(t_moldyn *moldyn) {
-
- t_ho_params *params;
- int i,j,count;
- t_atom *atom;
- t_3dvec distance;
- t_3dvec force;
- double d;
- double sc,equi_dist;
+ v3_scale(&force,&distance,-sc*(1.0-(equi_dist/d)));
+ v3_add(&(atom[i].f),&(atom[i].f),&force);
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+
+ /* neighbours not doing boundary condition overflow */
+ for(j=1;j<c;j++) {
+ this=&(neighbour[j]);
+ list_reset(this); /* there might not be a single atom */
+ if(this->start!=NULL) {
+
+ do {
+ btom=this->current->data;
+ v3_sub(&distance,&(atom[i].r),&(btom->r));
+ d=v3_norm(&distance);
+ if(d<=moldyn->cutoff) {
+ u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
+ v3_scale(&force,&distance,
+ -sc*(1.0-(equi_dist/d)));
+ v3_add(&(atom[i].f),&(atom[i].f),
+ &force);
+ }
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
- atom=moldyn->atom;
- count=moldyn->count;
- params=moldyn->pot_params;
- sc=params->spring_constant;
- equi_dist=params->equilibrium_distance;
+ }
+ }
- for(i=0;i<count;i++) v3_zero(&(atom[i].f));
+ /* neighbours due to boundary conditions */
+ for(j=c;j<27;j++) {
+ this=&(neighbour[j]);
+ list_reset(this); /* check boundary conditions */
+ if(this->start!=NULL) {
+
+ do {
+ btom=this->current->data;
+ v3_sub(&distance,&(atom[i].r),&(btom->r));
+ v3_per_bound(&distance,&(moldyn->dim));
+ d=v3_norm(&distance);
+ if(d<=moldyn->cutoff) {
+ u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
+ v3_scale(&force,&distance,
+ -sc*(1.0-(equi_dist/d)));
+ v3_add(&(atom[i].f),&(atom[i].f),
+ &force);
+ }
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
- for(i=0;i<count;i++) {
- for(j=0;j<i;j++) {
- v3_sub(&distance,&(atom[i].r),&(atom[j].r));
- v3_per_bound(&distance,&(moldyn->dim));
- d=v3_norm(&distance);
- if(d<=moldyn->cutoff) {
- v3_scale(&force,&distance,
- (-sc*(1.0-(equi_dist/d))));
- v3_add(&(atom[i].f),&(atom[i].f),&force);
- v3_sub(&(atom[j].f),&(atom[j].f),&force);
}
}
}
+ moldyn->energy=u;
+
return 0;
}
-
/* lennard jones potential & force for one sort of atoms */
-double potential_lennard_jones(t_moldyn *moldyn) {
+int lennard_jones(t_moldyn *moldyn) {
t_lj_params *params;
- t_atom *atom;
- int i,j;
+ t_atom *atom,*btom;
+ t_linkcell *lc;
+ t_list *this,neighbour[27];
+ int i,j,c;
int count;
- t_3dvec distance;
- double d,help;
- double u;
+ t_3dvec force,distance;
+ double d,h1,h2,u;
double eps,sig6,sig12;
+ int ni,nj,nk;
params=moldyn->pot_params;
atom=moldyn->atom;
+ lc=&(moldyn->lc);
count=moldyn->count;
- eps=params->epsilon;
+ eps=params->epsilon4;
sig6=params->sigma6;
sig12=params->sigma12;
u=0.0;
for(i=0;i<count;i++) {
- for(j=0;j<i;j++) {
- v3_sub(&distance,&(atom[j].r),&(atom[i].r));
- d=1.0/v3_absolute_square(&distance); /* 1/r^2 */
- help=d*d; /* 1/r^4 */
- help*=d; /* 1/r^6 */
- d=help*help; /* 1/r^12 */
- u+=eps*(sig12*d-sig6*help);
+ /* determine cell + neighbours */
+ ni=atom[i].r.x/lc->x;
+ nj=atom[i].r.y/lc->y;
+ nk=atom[i].r.z/lc->z;
+ c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
+
+ /* processing cell of atom i */
+ this=&(neighbour[0]);
+ list_reset(this); /* list has 1 element at minimum */
+ do {
+ btom=this->current->data;
+ if(btom==&(atom[i]))
+ continue;
+ v3_sub(&distance,&(atom[i].r),&(btom->r));
+ d=1.0/v3_absolute_square(&distance); /* 1/r^2 */
+ h1=d*d; /* 1/r^4 */
+ h2*=d; /* 1/r^6 */
+ h1=h2*h2; /* 1/r^12 */
+ u+=eps*(sig12*h1-sig6*h2);
+ h2*=d; /* 1/r^8 */
+ h1*=d; /* 1/r^14 */
+ h2*=6*sig6;
+ h1*=12*sig12;
+ d=-h1+h2;
+ d*=eps;
+ v3_scale(&force,&distance,d);
+ v3_add(&(atom[i].f),&(atom[i].f),&force);
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+
+ /* neighbours not doing boundary condition overflow */
+ for(j=1;j<c;j++) {
+ this=&(neighbour[j]);
+ list_reset(this); /* there might not be a single atom */
+ if(this->start!=NULL) {
+
+ do {
+ btom=this->current->data;
+ v3_sub(&distance,&(atom[i].r),&(btom->r));
+ d=v3_absolute_square(&distance); /* r^2 */
+ if(d<=moldyn->cutoff_square) {
+ d=1.0/d; /* 1/r^2 */
+ h1=d*d; /* 1/r^4 */
+ h2*=d; /* 1/r^6 */
+ h1=h2*h2; /* 1/r^12 */
+ u+=eps*(sig12*h1-sig6*h2);
+ h2*=d; /* 1/r^8 */
+ h1*=d; /* 1/r^14 */
+ h2*=6*sig6;
+ h1*=12*sig12;
+ d=-h1+h2;
+ d*=eps;
+ v3_scale(&force,&distance,d);
+ v3_add(&(atom[i].f),&(atom[i].f),
+ &force);
+ }
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+
+ }
}
- }
-
- return u;
-}
-
-int force_lennard_jones(t_moldyn *moldyn) {
-
- t_lj_params *params;
- int i,j,count;
- t_atom *atom;
- t_3dvec distance;
- t_3dvec force;
- double d,h1,h2;
- double eps,sig6,sig12;
- atom=moldyn->atom;
- count=moldyn->count;
- params=moldyn->pot_params;
- eps=params->epsilon;
- sig6=params->sigma6;
- sig12=params->sigma12;
-
- for(i=0;i<count;i++) v3_zero(&(atom[i].f));
+ /* neighbours due to boundary conditions */
+ for(j=c;j<27;j++) {
+ this=&(neighbour[j]);
+ list_reset(this); /* check boundary conditions */
+ if(this->start!=NULL) {
+
+ do {
+ btom=this->current->data;
+ v3_sub(&distance,&(atom[i].r),&(btom->r));
+ v3_per_bound(&distance,&(moldyn->dim));
+ d=v3_absolute_square(&distance); /* r^2 */
+ if(d<=moldyn->cutoff_square) {
+ d=1.0/d; /* 1/r^2 */
+ h1=d*d; /* 1/r^4 */
+ h2*=d; /* 1/r^6 */
+ h1=h2*h2; /* 1/r^12 */
+ u+=eps*(sig12*h1-sig6*h2);
+ h2*=d; /* 1/r^8 */
+ h1*=d; /* 1/r^14 */
+ h2*=6*sig6;
+ h1*=12*sig12;
+ d=-h1+h2;
+ d*=eps;
+ v3_scale(&force,&distance,d);
+ v3_add(&(atom[i].f),&(atom[i].f),
+ &force);
+ }
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
- for(i=0;i<count;i++) {
- for(j=0;j<i;j++) {
- v3_sub(&distance,&(atom[j].r),&(atom[i].r));
- v3_per_bound(&distance,&(moldyn->dim));
- d=v3_absolute_square(&distance);
- if(d<=moldyn->cutoff_square) {
- h1=1.0/d; /* 1/r^2 */
- d=h1*h1; /* 1/r^4 */
- h2=d*d; /* 1/r^8 */
- h1*=d; /* 1/r^6 */
- h1*=h2; /* 1/r^14 */
- h1*=sig12;
- h2*=sig6;
- d=12.0*h1-6.0*h2;
- d*=eps;
- v3_scale(&force,&distance,d);
- v3_add(&(atom[j].f),&(atom[j].f),&force);
- v3_sub(&(atom[i].f),&(atom[i].f),&force);
}
}
}
+ moldyn->energy=u;
+
return 0;
}