]> hackdaworld.org Git - physik/posic.git/commitdiff
implemented link cell method, segfaulting by now!
authorhackbard <hackbard>
Fri, 28 Apr 2006 15:02:29 +0000 (15:02 +0000)
committerhackbard <hackbard>
Fri, 28 Apr 2006 15:02:29 +0000 (15:02 +0000)
moldyn.c
moldyn.h
posic.c

index 17720675c4f5c0c76dea880ac40eb5787f772eb0..4b580890f3eaab417494fa497130f53570130ce4 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -36,6 +36,17 @@ int moldyn_usage(char **argv) {
        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;
@@ -44,6 +55,9 @@ int moldyn_usage(char **argv) {
 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));
 
@@ -51,6 +65,8 @@ int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
        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++) {
@@ -85,6 +101,44 @@ int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
                                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);
@@ -216,7 +270,7 @@ int destroy_lattice(t_atom *atom) {
        return 0;
 }
 
-int thermal_init(t_moldyn *moldyn,t_random *random,int count) {
+int thermal_init(t_moldyn *moldyn,t_random *random) {
 
        /*
         * - gaussian distribution of velocities
@@ -233,7 +287,7 @@ int thermal_init(t_moldyn *moldyn,t_random *random,int count) {
 
        /* gaussian distribution of velocities */
        v3_zero(&p_total);
-       for(i=0;i<count;i++) {
+       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);
@@ -250,19 +304,19 @@ int thermal_init(t_moldyn *moldyn,t_random *random,int count) {
        }
 
        /* 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(moldyn,count);
+       scale_velocity(moldyn);
 
        return 0;
 }
 
-int scale_velocity(t_moldyn *moldyn,int count) {
+int scale_velocity(t_moldyn *moldyn) {
 
        int i;
        double e,c;
@@ -274,10 +328,10 @@ int scale_velocity(t_moldyn *moldyn,int count) {
         * - 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*moldyn->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;
@@ -299,7 +353,7 @@ double get_e_kin(t_atom *atom,int count) {
 
 double get_e_pot(t_moldyn *moldyn) {
 
-       return(moldyn->potential(moldyn));
+       return moldyn->energy;
 }
 
 double get_total_energy(t_moldyn *moldyn) {
@@ -343,31 +397,7 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
  * numerical tricks
  */
 
-/* verlet list */
-
-int verlet_list_init(t_moldyn *moldyn) {
-
-       int i,fd;
-
-       fd=open("/dev/null",O_WRONLY);
-
-       for(i=0;i<moldyn->count;i++)
-               list_init(&(moldyn->atom[i].verlet),fd);
-
-       moldyn->r_verlet=1.1*moldyn->cutoff;
-       /* +moldyn->tau*\
-               sqrt(3.0*K_BOLTZMANN*moldyn->t/moldyn->atom[0].mass); */
-
-       printf("debug: r verlet = %.15f\n",moldyn->r_verlet);
-       printf("       r cutoff = %.15f\n",moldyn->cutoff);
-       printf("       dim      = %.15f\n",moldyn->dim.x);
-
-       /* make sure to update the list in the beginning */
-       moldyn->dr_max1=moldyn->r_verlet;
-       moldyn->dr_max2=moldyn->r_verlet;
-
-       return 0;
-}
+/* linked list / cell method */
 
 int link_cell_init(t_moldyn *moldyn) {
 
@@ -390,64 +420,27 @@ int link_cell_init(t_moldyn *moldyn) {
        return 0;
 }
 
-int verlet_list_update(t_moldyn *moldyn) {
-
-       int i,j;
-       t_3dvec d;
-       t_atom *atom;
-
-       atom=moldyn->atom;
-
-       puts("debug: list update start");
-
-       for(i=0;i<moldyn->count;i++) {
-               list_destroy(&(atom[i].verlet));
-               for(j=0;j<moldyn->count;j++) {
-                       if(i!=j) {
-                               v3_sub(&d,&(atom[i].r),&(atom[j].r));
-                               v3_per_bound(&d,&(moldyn->dim));
-                               if(v3_norm(&d)<=moldyn->r_verlet)
-                                       list_add_immediate_ptr(&(atom[i].verlet),&(atom[j]));
-                       }
-               }
-       }
-
-       moldyn->dr_max1=0.0;
-       moldyn->dr_max2=0.0;
-
-       puts("debug: list update end");
-
-       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;
-       nx=moldyn->lc.nx; ny=moldyn->lc.ny; nz=moldyn->lc.nz;
+       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++) {
-               for(i=0;i<nx;i++) {
-                       if((atom[count].r.x>=i*moldyn->dim.x) && \
-                          (atom[count].r.x<(i+1)*moldyn->dim.x))
-                               break;
-               }
-               for(j=0;j<ny;j++) {
-                       if((atom[count].r.y>=j*moldyn->dim.y) && \
-                          (atom[count].r.y<(j+1)*moldyn->dim.y))
-                               break;
-               }
-               for(k=0;k<nz;k++) {
-                       if((atom[count].r.z>=k*moldyn->dim.z) && \
-                          (atom[count].r.z<(k+1)*moldyn->dim.z))
-                               break;
-               }
+               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]));
        }
@@ -455,14 +448,58 @@ int link_cell_update(t_moldyn *moldyn) {
        return 0;
 }
 
-int verlet_list_shutdown(t_moldyn *moldyn) {
+int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
 
-       int i;
+       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;
 
-       for(i=0;i<moldyn->count;i++)
-               list_shutdown(&(moldyn->atom[i].verlet));
+       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 0;
+       return count2;
 }
 
 int link_cell_shutdown(t_moldyn *moldyn) {
@@ -491,7 +528,6 @@ int moldyn_integrate(t_moldyn *moldyn) {
        int i;
        unsigned int e,m,s,d,v;
        t_3dvec p;
-       double rlc;
 
        int fd;
        char fb[128];
@@ -503,20 +539,20 @@ int moldyn_integrate(t_moldyn *moldyn) {
        d=moldyn->dwrite;
        v=moldyn->vwrite;
 
-       /* verlet list */       
-       rlc=moldyn->r_verlet-moldyn->cutoff;
-
        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;
+
        /* create the neighbour list */
-       //verlet_list_update(moldyn);
        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 */
@@ -524,15 +560,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
                /* neighbour list update */
                link_cell_update(moldyn);
-               //if(moldyn->dr_max1+moldyn->dr_max2>rlc) {
-               //      printf("\n");
-               //      verlet_list_update(moldyn);
-               //}
 
                /* integration step */
                moldyn->integrate(moldyn);
 
-               /* check for log & visualiziation */
+               /* check for log & visualization */
                if(e) {
                        if(!(i%e))
                                dprintf(moldyn->efd,
@@ -581,41 +613,30 @@ int moldyn_integrate(t_moldyn *moldyn) {
 int velocity_verlet(t_moldyn *moldyn) {
 
        int i,count;
-       double tau,tau_square,dr;
+       double tau,tau_square;
        t_3dvec delta;
        t_atom *atom;
 
        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 */
                v3_scale(&delta,&(atom[i].v),tau);
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
-               v3_add(&(atom[i].dr),&(atom[i].dr),&delta);
                v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
-               v3_add(&(atom[i].dr),&(atom[i].dr),&delta);
                v3_per_bound(&(atom[i].r),&(moldyn->dim));
 
-               /* set maximum dr (possible list update) [obsolete] */
-               //dr=v3_norm(&(atom[i].dr));
-               //if(dr>moldyn->dr_max1) {
-               //      moldyn->dr_max2=moldyn->dr_max1;
-               //      moldyn->dr_max1=dr;
-               //}
-               //else if(dr>moldyn->dr_max2) moldyn->dr_max2=dr;
-
                /* velocities */
                v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
                v3_add(&(atom[i].v),&(atom[i].v),&delta);
        }
 
        /* forces depending on chosen potential */
-       moldyn->force(moldyn);
+       moldyn->potential_force_function(moldyn);
 
        for(i=0;i<count;i++) {
                /* again velocities */
@@ -635,84 +656,112 @@ int velocity_verlet(t_moldyn *moldyn) {
 
 /* 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,*btom;
        t_linkcell *lc;
-       int i;
+       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 i,j,k;
        int ni,nj,nk;
 
        params=moldyn->pot_params;
@@ -725,80 +774,101 @@ double potential_lennard_jones(t_moldyn *moldyn) {
 
        u=0.0;
        for(i=0;i<count;i++) {
-               /* verlet list [obsolete] */
-               //list_reset(&(atom[i].verlet));
-               //if(atom[i].verlet.current==NULL) continue;
-
-               /* determine cell */
+               /* determine cell + neighbours */
                ni=atom[i].r.x/lc->x;
                nj=atom[i].r.y/lc->y;
                nk=atom[i].r.z/lc->z;
-
-               while(1) {
-                       btom=atom[i].verlet.current->data;
+               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));
-                       v3_per_bound(&distance,&(moldyn->dim));
-                       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);
-                       if(list_next(&(atom[i].verlet))==L_NO_NEXT_ELEMENT)
-                               break;
+                       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,count;
-       t_atom *atom,*btom;
-       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->epsilon4;
-       sig6=6*params->sigma6;
-       sig12=12*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++) {
-               list_reset(&(atom[i].verlet));
-               if(atom[i].verlet.current==NULL) continue;
-               while(1) {
-                       btom=atom[i].verlet.current->data;
-                       v3_sub(&distance,&(atom[i].r),&(btom->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;
-                               /* actually there would be a '-',       *
-                                * but f=-d/dr potential                */
-                               d=h1+h2;
-                               d*=eps;
-                               v3_scale(&force,&distance,d);
-                               v3_add(&(atom[i].f),&(atom[i].f),&force);
-                               //v3_sub(&(atom[j].f),&(atom[j].f),&force);
                        }
-                       if(list_next(&(atom[i].verlet))==L_NO_NEXT_ELEMENT)
-                               break;
                }
        }
 
+       moldyn->energy=u;
+       
        return 0;
 }
 
index 53a22fdd0d0d3e2e3eb09caa451ea63383c89173..e9f1495355df4dad8a44f3779f678649b328a99d 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -18,10 +18,8 @@ typedef struct s_atom {
        t_3dvec r;      /* positions */
        t_3dvec v;      /* velocities */
        t_3dvec f;      /* forces */
-       t_3dvec dr;     /* delta r for verlet list updates */
        int element;    /* number of element in pse */
        double mass;    /* atom mass */
-       t_list verlet;  /* verlet neighbour list */
 } t_atom;
 
 typedef struct s_linkcell {
@@ -38,15 +36,11 @@ typedef struct s_moldyn {
        t_atom *atom;
        t_3dvec dim;
        /* potential, force & parameters */
-       double (*potential)(struct s_moldyn *moldyn);
-       int (*force)(struct s_moldyn *moldyn);
+       int (*potential_force_function)(struct s_moldyn *moldyn);
        void *pot_params;
-       /* cut off radius, verlet list & co */
+       /* cut off radius */
        double cutoff;
        double cutoff_square;
-       double r_verlet;
-       double dr_max1;
-       double dr_max2;
        /* linked list / cell method */
        t_linkcell lc;
        /* temperature */
@@ -55,6 +49,9 @@ typedef struct s_moldyn {
        int (*integrate)(struct s_moldyn *moldyn);
        int time_steps;
        double tau;
+       double tau_square;
+       /* energy */
+       double energy;
        /* logging & visualization */
        unsigned char lvstat;
        unsigned int ewrite;
@@ -107,6 +104,13 @@ typedef struct s_lj_params {
 #define MOLDYN_TAU                     1.0e-15
 #define MOLDYN_RUNS                    1000000
 
+#define MOLDYN_INTEGRATE_VERLET                0x00
+#define MOLDYN_INTEGRATE_DEFAULT MOLDYN_INTEGRATE_VERLET
+
+#define MOLDYN_POTENTIAL_HO            0x00
+#define MOLDYN_POTENTIAL_LJ            0x01
+#define MOLDYN_POTENTIAL_DEFAULT MOLDYN_POTENTIAL_LJ
+
 /* phsical values */
 
 #define K_BOLTZMANN            1.3807e-27                      /* Nm/K */
@@ -134,8 +138,8 @@ int moldyn_shutdown(t_moldyn *moldyn);
 int create_lattice(unsigned char type,int element,double mass,double lc,
                    int a,int b,int c,t_atom **atom);
 int destroy_lattice(t_atom *atom);
-int thermal_init(t_moldyn *moldyn,t_random *random,int count);
-int scale_velocity(t_moldyn *moldyn,int count);
+int thermal_init(t_moldyn *moldyn,t_random *random);
+int scale_velocity(t_moldyn *moldyn);
 double get_e_kin(t_atom *atom,int count);
 double get_e_pot(t_moldyn *moldyn);
 double get_total_energy(t_moldyn *moldyn);
@@ -143,19 +147,15 @@ t_3dvec get_total_p(t_atom *atom,int count);
 
 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t);
 
-int verlet_list_init(t_moldyn *moldyn);
 int link_cell_init(t_moldyn *moldyn);
-int verlet_list_update(t_moldyn *moldyn);
 int link_cell_update(t_moldyn *moldyn);
-int verlet_list_shutdown(t_moldyn *moldyn);
+int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell);
 int link_cell_shutdown(t_moldyn *moldyn);
 
 int moldyn_integrate(t_moldyn *moldyn);
 int velocity_verlet(t_moldyn *moldyn);
 
-double potential_harmonic_oscillator(t_moldyn *moldyn);
-int force_harmonic_oscillator(t_moldyn *moldyn);
-double potential_lennard_jones(t_moldyn *moldyn);
-int force_lennard_jones(t_moldyn *moldyn);
+int harmonic_oscillator(t_moldyn *moldyn);
+int lennard_jones(t_moldyn *moldyn);
 
 #endif
diff --git a/posic.c b/posic.c
index 7cde6e2acc4966d277c78397754bd7b11865581a..236e208bcce0bea5516bbe1a70d5efe1decda416 100644 (file)
--- a/posic.c
+++ b/posic.c
@@ -25,7 +25,6 @@ int main(int argc,char **argv) {
        double e;
        double help;
        t_3dvec p;
-       int count;
 
        t_lj_params lj;
        t_ho_params ho;
@@ -53,12 +52,12 @@ int main(int argc,char **argv) {
        vis.dim.y=b*LC_SI;
        vis.dim.z=c*LC_SI;
 
-       /* init lattice */
+       /* init lattice
        printf("placing silicon atoms ... ");
-       count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si);
-       printf("(%d) ok!\n",count);
-       /* testing purpose
-       count=2;
+       md.count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si);
+       printf("(%d) ok!\n",md.count);
+       testing purpose */
+       md.count=2;
        si=malloc(2*sizeof(t_atom));
        si[0].r.x=0.13*sqrt(3.0)*LC_SI/2.0;
        si[0].r.y=0;
@@ -70,22 +69,15 @@ int main(int argc,char **argv) {
        si[1].r.z=0;
        si[1].element=SI;
        si[1].mass=M_SI;
-       */
+       /* */
 
        /* moldyn init (now si is a valid address) */
-       md.count=count;
        md.atom=si;
-       md.potential=potential_lennard_jones;
-       md.force=force_lennard_jones;
-       //md.potential=potential_harmonic_oscillator;
-       //md.force=force_harmonic_oscillator;
+       md.potential_force_function=lennard_jones;
+       //md.potential_force_function=harmonic_oscillator;
        md.cutoff=R_CUTOFF*LC_SI;
-       md.cutoff_square=md.cutoff*md.cutoff;
        md.pot_params=&lj;
        //md.pot_params=&ho;
-       md.integrate=velocity_verlet;
-       //md.time_steps=RUNS;
-       //md.tau=TAU;
        md.status=0;
        md.visual=&vis;
        /* dimensions of the simulation cell */
@@ -93,23 +85,19 @@ int main(int argc,char **argv) {
        md.dim.y=b*LC_SI;
        md.dim.z=c*LC_SI;
 
-       /* verlet list init */
-       // later integrated in moldyn_init function!
-       verlet_list_init(&md);
-
        printf("setting thermal fluctuations (T=%f K)\n",md.t);
-       thermal_init(&md,&random,count);
-       //for(a=0;a<count;a++) v3_zero(&(si[0].v));
+       // thermal_init(&md,&random);
+       for(a=0;a<md.count;a++) v3_zero(&(si[0].v));
 
        /* check kinetic energy */
 
-       e=get_e_kin(si,count);
+       e=get_e_kin(si,md.count);
        printf("kinetic energy: %.40f [J]\n",e);
        printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
-              1.5*count*K_BOLTZMANN*md.t,md.t);
+              1.5*md.count*K_BOLTZMANN*md.t,md.t);
 
        /* check total momentum */
-       p=get_total_p(si,count);
+       p=get_total_p(si,md.count);
        printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
 
        /* potential paramters */
@@ -125,7 +113,6 @@ int main(int argc,char **argv) {
        printf("estimated accurate time step: %.30f [s]\n",
               estimate_time_step(&md,3.0,md.t));
 
-
        /*
         * let's do the actual md algorithm now
         *
@@ -139,7 +126,7 @@ int main(int argc,char **argv) {
 
        /* close */
 
-       verlet_list_shutdown(&md);
+       link_cell_shutdown(&md);
 
        rand_close(&random);