test

diff --git a/moldyn.c b/moldyn.c
index 505f496..2f274c2 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -41,6 +41,7 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 
 int moldyn_shutdown(t_moldyn *moldyn) {
 
+       moldyn_log_shutdown(moldyn);
        link_cell_shutdown(moldyn);
        moldyn_log_shutdown(moldyn);
        rand_close(&(moldyn->random));
@@ -382,15 +383,19 @@ t_3dvec get_total_p(t_moldyn *moldyn) {
        return p_total;
 }
 
-double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
+double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
 
        double tau;
 
-       tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
-       tau*=1.0E-9;
-       if(tau<moldyn->tau)
-               printf("[moldyn] warning: time step  (%f > %.15f)\n",
-                      moldyn->tau,tau);
+       /* nn_dist is the nearest neighbour distance */
+
+       if(moldyn->t==5.0) {
+               printf("[moldyn] i do not estimate timesteps below %f K!\n",
+                      MOLDYN_CRITICAL_EST_TEMP);
+               return 23.42;
+       }
+
+       tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
 
        return tau;     
 }
@@ -744,6 +749,7 @@ int potential_force_calc(t_moldyn *moldyn) {
        /* reset energy */
        moldyn->energy=0.0;
 
+printf("DEBUG: count = %d\n",count);
        for(i=0;i<count;i++) {
        
                /* reset force */
@@ -755,7 +761,8 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                /* 2 body pair potential/force */
                if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
-               
+       
+printf("DEBUG: processing atom %d\n",i);
                        link_cell_neighbour_index(moldyn,
                                (atom[i].r.x+moldyn->dim.x/2)/lc->x,
                                (atom[i].r.y+moldyn->dim.y/2)/lc->y,
@@ -765,6 +772,7 @@ int potential_force_calc(t_moldyn *moldyn) {
                        countn=lc->countn;
                        dnlc=lc->dnlc;
 
+printf("DEBUG: countn = %d - dnslc = %d\n",countn,dnlc);
                        for(j=0;j<countn;j++) {
 
                                this=&(neighbour[j]);
@@ -783,6 +791,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                        if((btom->attr&ATOM_ATTR_2BP)&
                                           (atom[i].attr&ATOM_ATTR_2BP))
+printf("DEBUG: calling func2b\n");
                                                moldyn->func2b(moldyn,
                                                               &(atom[i]),
                                                               btom,

diff --git a/moldyn.h b/moldyn.h
index 58302fb..947467c 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -118,12 +118,6 @@ typedef struct s_moldyn {
        t_random random;        /* random interface */
 } t_moldyn;
 
-#define MOLDYN_LVSTAT_TOTAL_E          0x01
-#define MOLDYN_LVSTAT_TOTAL_M          0x02
-#define MOLDYN_LVSTAT_SAVE             0x04
-#define MOLDYN_LVSTAT_VISUAL           0x08
-#define MOLDYN_LVSTAT_INITIALIZED      0x10
-
 #define MOLDYN_STAT_PBX                        0x08    /* periodic boudaries in x */
 #define MOLDYN_STAT_PBY                        0x10    /* y */
 #define MOLDYN_STAT_PBZ                        0x20    /* and z direction */
@@ -231,6 +225,8 @@ typedef struct s_tersoff_mult_params {
 #define MOLDYN_CUTOFF                  1.0e-9
 #define MOLDYN_RUNS                    1000000
 
+#define MOLDYN_CRITICAL_EST_TEMP       5.0
+
 #define MOLDYN_INTEGRATE_VERLET                0x00
 #define MOLDYN_INTEGRATE_DEFAULT       MOLDYN_INTEGRATE_VERLET
 
@@ -307,7 +303,7 @@ double get_e_pot(t_moldyn *moldyn);
 double get_total_energy(t_moldyn *moldyn);
 t_3dvec get_total_p(t_moldyn *moldyn);
 
-double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t);
+double estimate_time_step(t_moldyn *moldyn,double nn_dist);
 
 int link_cell_init(t_moldyn *moldyn);
 int link_cell_update(t_moldyn *moldyn);

diff --git a/sic.c b/sic.c
index 9c98857..3924375 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -30,6 +30,9 @@ int main(int argc,char **argv) {
        /* misc variables, mainly to initialize stuff */
        t_3dvec r,v;
 
+       /* temperature */
+       double t;
+
        /* initialize moldyn */
        printf("[sic] moldyn init\n");
        moldyn_init(&md,argc,argv);
@@ -91,12 +94,12 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        printf("[sic] adding schedule\n");
-       moldyn_add_schedule(&md,1000,1.0e-15);
+       moldyn_add_schedule(&md,10000,1.0e-12);
 
        /* activate logging */
        printf("[sic] activate logging\n");
        moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",100);
-       moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",50);
+       moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",100);
 
        /*
         * let's do the actual md algorithm now