ready for tersoff 3bp debugging

diff --git a/moldyn.c b/moldyn.c
index 80f4fa1..af0132d 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -71,9 +71,18 @@ int set_cutoff(t_moldyn *moldyn,double cutoff) {
        return 0;
 }
 
-int set_temperature(t_moldyn *moldyn,double t) {
+int set_temperature(t_moldyn *moldyn,double t_ref) {
 
-       moldyn->t=t;
+       moldyn->t_ref=t_ref;
+
+       return 0;
+}
+
+int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
+
+       moldyn->pt_scale=(ptype|ttype);
+       moldyn->t_tc=ttc;
+       moldyn->p_tc=ptc;
 
        return 0;
 }
@@ -278,7 +287,7 @@ int destroy_atoms(t_moldyn *moldyn) {
        return 0;
 }
 
-int thermal_init(t_moldyn *moldyn) {
+int thermal_init(t_moldyn *moldyn,u8 equi_init) {
 
        /*
         * - gaussian distribution of velocities
@@ -298,7 +307,7 @@ int thermal_init(t_moldyn *moldyn) {
        /* gaussian distribution of velocities */
        v3_zero(&p_total);
        for(i=0;i<moldyn->count;i++) {
-               sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
+               sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
                /* x direction */
                v=sigma*rand_get_gauss(random);
                atom[i].v.x=v;
@@ -321,12 +330,12 @@ int thermal_init(t_moldyn *moldyn) {
        }
 
        /* velocity scaling */
-       scale_velocity(moldyn,VSCALE_INIT_EQUI);
+       scale_velocity(moldyn,equi_init);
 
        return 0;
 }
 
-int scale_velocity(t_moldyn *moldyn,u8 type) {
+int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 
        int i;
        double e,scale;
@@ -339,29 +348,40 @@ int scale_velocity(t_moldyn *moldyn,u8 type) {
         * - velocity scaling (E = 3/2 N k T), E: kinetic energy
         */
 
+       /* get kinetic energy / temperature & count involved atoms */
        e=0.0;
        count=0;
        for(i=0;i<moldyn->count;i++) {
-               if(atom[i].attr&ATOM_ATTR_HB) {
+               if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
                        e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
                        count+=1;
                }
        }
-
-       /* temporary hack for e,t = 0 */
+       if(count!=0) moldyn->t=(2.0*e)/(3.0*count*K_BOLTZMANN);
+       else return 0;  /* no atoms involved in scaling! */
+       
+       /* (temporary) hack for e,t = 0 */
        if(e==0.0) {
-               if(moldyn->t!=0.0)
-                       thermal_init(moldyn);
+       moldyn->t=0.0;
+               if(moldyn->t_ref!=0.0)
+                       thermal_init(moldyn,equi_init);
                else
-                       return 0;
+                       return 0; /* no scaling needed */
        }
 
-       /* direct scaling */
-       scale=(1.5*count*K_BOLTZMANN*moldyn->t)/e;
-       if(type&VSCALE_INIT_EQUI) scale*=2.0; /* equipartition theorem */
+
+       /* get scaling factor */
+       scale=moldyn->t_ref/moldyn->t;
+       if(equi_init&TRUE)
+               scale*=2.0;
+       else
+               if(moldyn->pt_scale&T_SCALE_BERENDSEN)
+                       scale=1.0+moldyn->tau*(scale-1.0)/moldyn->t_tc;
        scale=sqrt(scale);
+
+       /* velocity scaling */
        for(i=0;i<moldyn->count;i++)
-               if(atom[i].attr&ATOM_ATTR_HB)
+               if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
                        v3_scale(&(atom[i].v),&(atom[i].v),scale);
 
        return 0;
@@ -666,6 +686,10 @@ int moldyn_integrate(t_moldyn *moldyn) {
                /* integration step */
                moldyn->integrate(moldyn);
 
+               /* p/t scaling */
+               if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
+                       scale_velocity(moldyn,FALSE);
+
                /* increase absolute time */
                moldyn->time+=moldyn->tau;
 
@@ -777,7 +801,9 @@ int potential_force_calc(t_moldyn *moldyn) {
        int i,j,k,count;
        t_atom *itom,*jtom,*ktom;
        t_linkcell *lc;
-       t_list neighbour_i[27],neighbour_j[27];
+       t_list neighbour_i[27];
+       t_list neighbour_i2[27];
+       //t_list neighbour_j[27];
        t_list *this,*that;
        u8 bc_ij,bc_ijk;
        int countn,dnlc;
@@ -839,47 +865,61 @@ int potential_force_calc(t_moldyn *moldyn) {
                                           !(jtom->attr&ATOM_ATTR_3BP))
                                                continue;
 
-                                       link_cell_neighbour_index(moldyn,
-                                          (jtom->r.x+moldyn->dim.x/2)/lc->x,
-                                          (jtom->r.y+moldyn->dim.y/2)/lc->y,
-                                          (jtom->r.z+moldyn->dim.z/2)/lc->z,
-                                          neighbour_j);
-
-                                       /* neighbours of j */
-                                       for(k=0;k<lc->countn;k++) {
-
-                                               that=&(neighbour_j[k]);
-                                               list_reset(that);
-                                       
-                                               if(that->start==NULL)
-                                                       continue;
-
-                                               bc_ijk=(k<lc->dnlc)?0:1;
-
-                                               do {
-
-                       ktom=that->current->data;
-
-                       if(!(ktom->attr&ATOM_ATTR_3BP))
-                               continue;
-
-                       if(ktom==jtom)
-                               continue;
-
-                       if(ktom==&(itom[i]))
-                               continue;
-
-                       moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
-
-                                               } while(list_next(that)!=\
-                                                       L_NO_NEXT_ELEMENT);
-
-                                       }
-                                       
-                                       /* neighbours of i */
+                       /*
+                        * according to mr. nordlund, we dont need to take the 
+                        * sum over all atoms now, as 'this is centered' around
+                        * atom i ...
+                        * i am not quite sure though! there is a not vanishing
+                        * part even if f_c_ik is zero ...
+                        * this analytical potentials suck!
+                        * switching from mc to md to dft soon!
+                        */
+
+                       //              link_cell_neighbour_index(moldyn,
+                       //                 (jtom->r.x+moldyn->dim.x/2)/lc->x,
+                       //                 (jtom->r.y+moldyn->dim.y/2)/lc->y,
+                       //                 (jtom->r.z+moldyn->dim.z/2)/lc->z,
+                       //                 neighbour_j);
+
+//                                     /* neighbours of j */
+//                                     for(k=0;k<lc->countn;k++) {
+//
+//                                             that=&(neighbour_j[k]);
+//                                             list_reset(that);
+//                                     
+//                                             if(that->start==NULL)
+//                                                     continue;
+//
+//                                             bc_ijk=(k<lc->dnlc)?0:1;
+//
+//                                             do {
+//
+//                     ktom=that->current->data;
+//
+//                     if(!(ktom->attr&ATOM_ATTR_3BP))
+//                             continue;
+//
+//                     if(ktom==jtom)
+//                             continue;
+//
+//                     if(ktom==&(itom[i]))
+//                             continue;
+//
+//                     moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+//
+/*                                             } while(list_next(that)!=\ */
+//                                                     L_NO_NEXT_ELEMENT);
+//
+//                                     }
+                       
+                                       /* copy the neighbour lists */
+                                       memcpy(neighbour_i2,neighbour_i,
+                                              27*sizeof(t_list));
+
+                                       /* get neighbours of i */
                                        for(k=0;k<countn;k++) {
 
-                                               that=&(neighbour_i[k]);
+                                               that=&(neighbour_i2[k]);
                                                list_reset(that);
                                        
                                                if(that->start==NULL)
@@ -900,7 +940,9 @@ int potential_force_calc(t_moldyn *moldyn) {
                        if(ktom==&(itom[i]))
                                continue;
 
+printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
                        moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
 
                                                } while(list_next(that)!=\
                                                        L_NO_NEXT_ELEMENT);
@@ -1094,8 +1136,11 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        if(bc) check_per_bound(moldyn,&dist_ij);
 
-       /* save for use in 3bp */ /* REALLY ?!?!?! */
-       exchange->dist_ij=dist_ij;
+       d_ij=v3_norm(&dist_ij);
+
+       /* save for use in 3bp */
+       exchange->dist_ij=dist_ij; /* <- needed ? */
+       exchange->d_ij=d_ij;
 
        /* constants */
        if(num==aj->bnum) {
@@ -1121,11 +1166,6 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                params->exchange.chi=params->chi;
        }
 
-       d_ij=v3_norm(&dist_ij);
-
-       /* save for use in 3bp */
-       exchange->d_ij=d_ij;
-
        if(d_ij>S)
                return 0;
 
@@ -1153,8 +1193,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        /* add forces */
        v3_add(&(ai->f),&(ai->f),&force);
-       /* energy is 0.5 f_r f_c, but we will sum it up twice ... */
-       moldyn->energy+=(0.25*f_r*f_c);
+       /* energy is 0.5 f_r f_c ... */
+       moldyn->energy+=(0.5*f_r*f_c);
 
        /* save for use in 3bp */
        exchange->f_c=f_c;

diff --git a/moldyn.h b/moldyn.h
index 2d4b563..8b40d4c 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -85,7 +85,16 @@ typedef struct s_moldyn {
 
        t_linkcell lc;          /* linked cell list interface */
 
-       double t;               /* temperature */
+       double t_ref;           /* reference temperature */
+       double t;               /* actual temperature */
+
+       double p_ref;           /* reference pressure */
+       double p;               /* actual pressure */
+
+       /* pressure and temperature control (velocity/volume scaling) */
+       unsigned char pt_scale; /* type of p and t scaling */
+       double t_tc;            /* t berendsen control time constant */
+       double p_tc;            /* p berendsen control time constant */
 
        /* simulation schedule */
        t_moldyn_schedule schedule;
@@ -127,8 +136,10 @@ typedef struct s_moldyn {
 #define MOLDYN_2BP                     0x01    /* 2 body */
 #define MOLDYN_3BP                     0x02    /* and 3 body particle pots */
 
-#define VSCALE_INIT_EQUI               0x00    /* initial, eq positions */
-#define VSCALE_BERENDSEN               0x01    /* berendsen control */
+#define T_SCALE_BERENDSEN              0x01    /* berendsen t control */
+#define T_SCALE_DIRECT                 0x02    /* direct t control */
+#define P_SCALE_BERENDSEN              0x04    /* berendsen p control */
+#define P_SCALE_DIRECT                 0x08    /* direct p control */
 
 
 /*
@@ -310,7 +321,8 @@ int moldyn_shutdown(t_moldyn *moldyn);
 
 int set_int_alg(t_moldyn *moldyn,u8 algo);
 int set_cutoff(t_moldyn *moldyn,double cutoff);
-int set_temperature(t_moldyn *moldyn,double t);
+int set_temperature(t_moldyn *moldyn,double t_ref);
+int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc);
 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize);
 int set_nn_dist(t_moldyn *moldyn,double dist);
 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z);
@@ -327,8 +339,8 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
              t_3dvec *r,t_3dvec *v);
 int destroy_atoms(t_moldyn *moldyn);
 
-int thermal_init(t_moldyn *moldyn);
-int scale_velocity(t_moldyn *moldyn,u8 type);
+int thermal_init(t_moldyn *moldyn,u8 equi_init);
+int scale_velocity(t_moldyn *moldyn,u8 equi_init);
 
 double get_e_kin(t_moldyn *moldyn);
 double get_e_pot(t_moldyn *moldyn);

diff --git a/sic.c b/sic.c
index fa1cb66..b10c681 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -23,6 +23,12 @@ int main(int argc,char **argv) {
        t_ho_params ho;
        t_tersoff_mult_params tp;
 
+       /* misc parameters */
+       double tau;
+
+       /* values */
+       tau=1.0e-15;    /* delta t = 1 fs */
+
        /* initialize moldyn */
        printf("[sic] moldyn init\n");
        moldyn_init(&md,argc,argv);
@@ -106,11 +112,15 @@ int main(int argc,char **argv) {
 
        /* set temperature */
        printf("[sic] setting temperature\n");
-       set_temperature(&md,273.0);
+       set_temperature(&md,0.0);
+
+       /* set p/t scaling */
+       printf("[sic] set p/t scaling\n");
+       set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100*tau);
        
-       /* initial thermal fluctuations of particles */
+       /* initial thermal fluctuations of particles (in equilibrium) */
        printf("[sic] thermal init\n");
-       thermal_init(&md);
+       thermal_init(&md,TRUE);
 
        /* create the simulation schedule */
        printf("[sic] adding schedule\n");