added del_atoms, playing with Ef calc, some more TODO

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index bc4fdff..7698c24 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -736,6 +736,9 @@ int del_atom(t_moldyn *moldyn,int tag) {
 
        t_atom *new,*old;
        int cnt;
+#if defined LOWMEM_LISTS || defined PTHREADS
+       void *ptr;
+#endif
 
        old=moldyn->atom;
 
@@ -758,6 +761,28 @@ int del_atom(t_moldyn *moldyn,int tag) {
 
        free(old);
 
+#ifdef LOWMEM_LISTS
+       ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int));
+       if(!ptr) {
+               perror("[moldyn] list realloc (del atom)");
+               return -1;
+       }
+       moldyn->lc.subcell->list=ptr;
+       // update lists
+       link_cell_update(moldyn);
+#endif
+
+#ifdef PTHREADS
+       ptr=realloc(amutex,moldyn->count*sizeof(pthread_mutex_t));
+       if(!ptr) {
+               perror("[moldyn] mutex realloc (add atom)");
+               return -1;
+       }
+       amutex=ptr;
+       pthread_mutex_destroy(&(amutex[moldyn->count+1]));
+#endif
+
+
        return 0;
 }
 
@@ -2065,10 +2090,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
                if(e) {
                        if(!(moldyn->total_steps%e))
                                dprintf(moldyn->efd,
-                                       "%f %f %f %f\n",
+                                       "%f %f %f %f %f %f\n",
                                        moldyn->time,moldyn->ekin/energy_scale,
                                        moldyn->energy/energy_scale,
-                                       get_total_energy(moldyn)/energy_scale);
+                                       get_total_energy(moldyn)/energy_scale,
+                                       moldyn->ekin/EV,moldyn->energy/EV);
                }
                if(m) {
                        if(!(moldyn->total_steps%m)) {
@@ -3203,9 +3229,8 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
        ba=data;
 
        /* increase total bond counter
-        * ... double counting!
         */
-       ba->tcnt+=2;
+       ba->tcnt+=1;
 
        if(itom->brand==0)
                ba->acnt[jtom->tag]+=1;
@@ -3222,10 +3247,11 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 
 int bond_analyze(t_moldyn *moldyn,double *quality) {
 
-       // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
-
-       int qcnt;
-       int ccnt,cset;
+       int qcnt4;
+       int qcnt3;
+       int ccnt4;
+       int ccnt3;
+       int bcnt;
        t_ba ba;
        int i;
        t_atom *atom;
@@ -3245,9 +3271,9 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
        memset(ba.bcnt,0,moldyn->count*sizeof(int));
 
        ba.tcnt=0;
-       qcnt=0;
-       ccnt=0;
-       cset=0;
+       qcnt4=0; qcnt3=0;
+       ccnt4=0; ccnt3=0;
+       bcnt=0;
 
        atom=moldyn->atom;
 
@@ -3256,27 +3282,30 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
        for(i=0;i<moldyn->count;i++) {
                if(atom[i].brand==0) {
                        if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
-                               qcnt+=4;
+                               qcnt4+=4;
+                       if((ba.acnt[i]==0)&(ba.bcnt[i]==3))
+                               qcnt3+=3;
                }
                else {
                        if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
-                               qcnt+=4;
-                               ccnt+=1;
+                               qcnt4+=4;
+                               ccnt4+=1;
                        }
-                       cset+=1;
+                       if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) {
+                               qcnt3+=4;
+                               ccnt3+=1;
+                       }
+                       bcnt+=1;
                }
        }
 
-       printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
-       printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
-
        if(quality) {
-               quality[0]=1.0*ccnt/cset;
-               quality[1]=1.0*qcnt/ba.tcnt;
+               quality[0]=1.0*ccnt4/bcnt;
+               quality[1]=1.0*ccnt3/bcnt;
        }
        else {
-               printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
-               printf("[moldyn] bond analyze:   tot_q=%f\n",1.0*qcnt/ba.tcnt);
+               printf("[moldyn] bond analyze: %f %f\n",
+                      1.0*ccnt4/bcnt,1.0*ccnt3/bcnt);
        }
 
        return 0;