pot hopefully ok, doing timestep tests for 450C now ...

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index f8bfd3f..5aa6d87 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -767,23 +767,23 @@ int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
        t_moldyn_schedule *schedule;
 
        schedule=&(moldyn->schedule);
-       count=++(schedule->content_count);
+       count=++(schedule->total_sched);
 
-       ptr=realloc(moldyn->schedule.runs,count*sizeof(int));
+       ptr=realloc(schedule->runs,count*sizeof(int));
        if(!ptr) {
                perror("[moldyn] realloc (runs)");
                return -1;
        }
-       moldyn->schedule.runs=ptr;
-       moldyn->schedule.runs[count-1]=runs;
+       schedule->runs=ptr;
+       schedule->runs[count-1]=runs;
 
        ptr=realloc(schedule->tau,count*sizeof(double));
        if(!ptr) {
                perror("[moldyn] realloc (tau)");
                return -1;
        }
-       moldyn->schedule.tau=ptr;
-       moldyn->schedule.tau[count-1]=tau;
+       schedule->tau=ptr;
+       schedule->tau[count-1]=tau;
 
        return 0;
 }
@@ -806,16 +806,16 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) {
 
 int moldyn_integrate(t_moldyn *moldyn) {
 
-       int i,sched;
+       int i;
        unsigned int e,m,s,v;
        t_3dvec p;
-       t_moldyn_schedule *schedule;
+       t_moldyn_schedule *sched;
        t_atom *atom;
        int fd;
        char dir[128];
        double ds;
 
-       schedule=&(moldyn->schedule);
+       sched=&(moldyn->schedule);
        atom=moldyn->atom;
 
        /* initialize linked cell method */
@@ -852,12 +852,12 @@ int moldyn_integrate(t_moldyn *moldyn) {
        moldyn->debug=0;
 
        /* executing the schedule */
-       for(sched=0;sched<moldyn->schedule.content_count;sched++) {
+       for(sched->count=0;sched->count<sched->total_sched;sched->count++) {
 
                /* setting amount of runs and finite time step size */
-               moldyn->tau=schedule->tau[sched];
+               moldyn->tau=sched->tau[sched->count];
                moldyn->tau_square=moldyn->tau*moldyn->tau;
-               moldyn->time_steps=schedule->runs[sched];
+               moldyn->time_steps=sched->runs[sched->count];
 
        /* integration according to schedule */
 
@@ -907,7 +907,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
                                visual_atoms(&(moldyn->vis),moldyn->time,
                                             moldyn->atom,moldyn->count);
                                printf("\rsched: %d, steps: %d, debug: %d",
-                                      sched,i,moldyn->debug);
+                                      sched->count,i,moldyn->debug);
                                fflush(stdout);
                        }
                }
@@ -918,8 +918,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
        }
 
                /* check for hooks */
-               if(schedule->hook)
-                       schedule->hook(moldyn,schedule->hook_params);
+               if(sched->hook)
+                       sched->hook(moldyn,sched->hook_params);
 
                /* get a new info line */
                printf("\n");
@@ -1108,6 +1108,9 @@ int potential_force_calc(t_moldyn *moldyn) {
                }
 
        }
+#ifdef DEBUG
+printf("\n\n");
+#endif
 
        return 0;
 }
@@ -1397,6 +1400,14 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        /* add forces of 2bp (ij, ji) contribution
         * dVij = dVji and we sum up both: no 1/2) */
        v3_add(&(ai->f),&(ai->f),&force);
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+       printf("dVij, dVji (2bp) contrib:\n");
+       printf("%f | %f\n",force.x,ai->f.x);
+       printf("%f | %f\n",force.y,ai->f.y);
+       printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
 
        /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
        moldyn->energy+=(0.5*f_r*f_c);
@@ -1488,6 +1499,14 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        /* add force */
        v3_add(&(ai->f),&(ai->f),&force);
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+       printf("dVij (3bp) contrib:\n");
+       printf("%f | %f\n",force.x,ai->f.x);
+       printf("%f | %f\n",force.y,ai->f.y);
+       printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
 
        /* add energy of 3bp sum */
        moldyn->energy+=(0.5*f_c*b*f_a);
@@ -1516,6 +1535,14 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        /* add force */
        v3_add(&(ai->f),&(ai->f),&force);
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+       printf("dVji (3bp) contrib:\n");
+       printf("%f | %f\n",force.x,ai->f.x);
+       printf("%f | %f\n",force.y,ai->f.y);
+       printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
 
        return 0;
 }
@@ -1777,6 +1804,14 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
                v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
                                                  /* scaled with 0.5 ^ */
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+       printf("dVjk (3bp) contrib:\n");
+       printf("%f | %f\n",temp2.x,ai->f.x);
+       printf("%f | %f\n",temp2.y,ai->f.y);
+       printf("%f | %f\n",temp2.z,ai->f.z);
+}
+#endif
 
        }