added del_atom function

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 92c4e58..41e022a 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -558,6 +558,35 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        return 0;
 }
 
+int del_atom(t_moldyn *moldyn,int tag) {
+
+       t_atom *new,*old;
+       int cnt;
+
+       old=moldyn->atom;
+
+       new=(t_atom *)malloc((moldyn->count-1)*sizeof(t_atom));
+       if(!new) {
+               perror("[moldyn]malloc (del atom)");
+               return -1;
+       }
+
+       for(cnt=0;cnt<tag;cnt++)
+               new[cnt]=old[cnt];
+       
+       for(cnt=tag+1;cnt<moldyn->count;cnt++) {
+               new[cnt-1]=old[cnt];
+               new[cnt-1].tag=cnt-1;
+       }
+
+       moldyn->count-=1;
+       moldyn->atom=new;
+
+       free(old);
+
+       return 0;
+}
+
 /* cubic init */
 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 
@@ -1371,9 +1400,6 @@ int moldyn_integrate(t_moldyn *moldyn) {
        moldyn->tau_square=moldyn->tau*moldyn->tau;
        moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
 
-       /* energy scaling factor */
-       energy_scale=moldyn->count*EV;
-
        /* get current time */
        gettimeofday(&t1,NULL);
 
@@ -1413,6 +1439,9 @@ return 0;
                moldyn->tau_square=moldyn->tau*moldyn->tau;
                moldyn->time_steps=sched->runs[sched->count];
 
+               /* energy scaling factor (might change!) */
+               energy_scale=moldyn->count*EV;
+
        /* integration according to schedule */
 
        for(i=0;i<moldyn->time_steps;i++) {
@@ -1435,7 +1464,7 @@ return 0;
 
                /* check for log & visualization */
                if(e) {
-                       if(!(i%e))
+                       if(!(moldyn->total_steps%e))
                                dprintf(moldyn->efd,
                                        "%f %f %f %f\n",
                                        moldyn->time,moldyn->ekin/energy_scale,
@@ -1443,7 +1472,7 @@ return 0;
                                        get_total_energy(moldyn)/energy_scale);
                }
                if(m) {
-                       if(!(i%m)) {
+                       if(!(moldyn->total_steps%m)) {
                                momentum=get_total_p(moldyn);
                                dprintf(moldyn->mfd,
                                        "%f %f %f %f %f\n",moldyn->time,
@@ -1452,7 +1481,7 @@ return 0;
                        }
                }
                if(p) {
-                       if(!(i%p)) {
+                       if(!(moldyn->total_steps%p)) {
                                dprintf(moldyn->pfd,
                                        "%f %f %f %f %f\n",moldyn->time,
                                         moldyn->p/BAR,moldyn->p_avg/BAR,
@@ -1460,14 +1489,14 @@ return 0;
                        }
                }
                if(t) {
-                       if(!(i%t)) {
+                       if(!(moldyn->total_steps%t)) {
                                dprintf(moldyn->tfd,
                                        "%f %f %f\n",
                                        moldyn->time,moldyn->t,moldyn->t_avg);
                        }
                }
                if(s) {
-                       if(!(i%s)) {
+                       if(!(moldyn->total_steps%s)) {
                                snprintf(dir,128,"%s/s-%07.f.save",
                                         moldyn->vlsdir,moldyn->time);
                                fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
@@ -1482,14 +1511,14 @@ return 0;
                        }       
                }
                if(v) {
-                       if(!(i%v)) {
+                       if(!(moldyn->total_steps%v)) {
                                visual_atoms(&(moldyn->vis),moldyn->time,
                                             moldyn->atom,moldyn->count);
                        }
                }
 
                /* display progress */
-               if(!(i%10)) {
+               if(!(moldyn->total_steps%10)) {
                        /* get current time */
                        gettimeofday(&t2,NULL);