X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=2f274c2da8fb7fcaeec9344f21bec718cb8585a8;hb=8e98c9cbbddf1c9c852874f956b0dbe2f9a4a922;hp=505f4960a18037b6d043ddf0d7596625660b328a;hpb=7aface2e184e527132819f17b0fa18a6e1ba4bc3;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 505f496..2f274c2 100644
--- 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,