Merge branch 'leadoff'

Conflicts:
	Makefile
	moldyn.c
	potentials/albe.c
	potentials/albe.h
	runmd

 ... hopefully fixed! Just a test to merge leadoff back to master! :)

diff --cc mdrun.c
Simple merge

diff --cc moldyn.c
index 010b8ed,33224a6..1d6e0b7
--- 1/moldyn.c
--- 2/moldyn.c
+++ b/moldyn.c
@@@ -216,28 -263,18 +264,28 @@@ int set_potential(t_moldyn *moldyn,u8 t
  
        switch(type) {
                case MOLDYN_POTENTIAL_TM:
-                       moldyn->func_i0=tersoff_mult_1bp;
-                       moldyn->func_j1=tersoff_mult_3bp_j1;
-                       moldyn->func_j1_k0=tersoff_mult_3bp_k1;
-                       moldyn->func_j1c=tersoff_mult_3bp_j2;
-                       moldyn->func_j1_k1=tersoff_mult_3bp_k2;
+                       //moldyn->func1b=tersoff_mult_1bp;
+                       moldyn->func3b_j1=tersoff_mult_3bp_j1;
+                       moldyn->func3b_k1=tersoff_mult_3bp_k1;
+                       moldyn->func3b_j2=tersoff_mult_3bp_j2;
+                       moldyn->func3b_k2=tersoff_mult_3bp_k2;
                        moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
                        break;
 +              case MOLDYN_POTENTIAL_AO:
 +                      moldyn->func_j1=albe_orig_mult_3bp_j1;
 +                      moldyn->func_j1_k0=albe_orig_mult_3bp_k1;
 +                      moldyn->func_j1c=albe_orig_mult_3bp_j2;
 +                      moldyn->func_j1_k1=albe_orig_mult_3bp_k2;
 +                      moldyn->check_2b_bond=albe_orig_mult_check_2b_bond;
 +                      break;
                case MOLDYN_POTENTIAL_AM:
 -                      moldyn->func3b_j1=albe_mult_3bp_j1;
 -                      moldyn->func3b_k1=albe_mult_3bp_k1;
 -                      moldyn->func3b_j2=albe_mult_3bp_j2;
 -                      moldyn->func3b_k2=albe_mult_3bp_k2;
 +                      moldyn->func_i0=albe_mult_i0;
 +                      moldyn->func_j0=albe_mult_i0_j0;
 +                      moldyn->func_j0_k0=albe_mult_i0_j0_k0;
 +                      moldyn->func_j0e=albe_mult_i0_j1;
 +                      moldyn->func_j1=albe_mult_i0_j2;
 +                      moldyn->func_j1_k0=albe_mult_i0_j2_k0;
 +                      moldyn->func_j1c=albe_mult_i0_j3;
                        moldyn->check_2b_bond=albe_mult_check_2b_bond;
                        break;
                case MOLDYN_POTENTIAL_HO:
@@@ -1908,86 -2536,17 +2553,74 @@@ int potential_force_calc(t_moldyn *mold
                                        if(jtom==&(itom[i]))
                                                continue;
  
 +                                      /* reset 3bp run */
 +                                      moldyn->run3bp=1;
 +
                                        if((jtom->attr&ATOM_ATTR_2BP)&
                                           (itom[i].attr&ATOM_ATTR_2BP)) {
 -                                              moldyn->func2b(moldyn,
 -                                                             &(itom[i]),
 -                                                             jtom,
 -                                                             bc_ij);
 +                                              moldyn->func_j0(moldyn,
 +                                                              &(itom[i]),
 +                                                              jtom,
 +                                                              bc_ij);
 +                                      }
 +
 +                                      /* 3 body potential/force */
 +
 +                                      /* in j loop, 3bp run can be skipped */
 +                                      if(!(moldyn->run3bp))
 +                                              continue;
 +
 +                                      if(!(itom[i].attr&ATOM_ATTR_3BP))
 +                                              continue;
 +
 +                                      if(!(jtom->attr&ATOM_ATTR_3BP))
 +                                              continue;
 +
 +                                      if(moldyn->func_j0_k0==NULL)
 +                                              continue;
 +
 +                              /* first loop over atoms k in first j loop */
 +                              for(k=0;k<27;k++) {
 +
 +                                      bc_ik=(k<dnlc)?0:1;
 +#ifdef STATIC_LISTS
 +                                      q=0;
 +
 +                                      while(neighbour_i[j][q]!=0) {
 +
 +                                              ktom=&(atom[neighbour_i[k][q]]);
 +                                              q++;
 +#else
 +                                      that=&(neighbour_i2[k]);
 +                                      list_reset_f(that);
 +                                      
 +                                      if(that->start==NULL)
 +                                              continue;
 +
 +                                      do {
 +                                              ktom=that->current->data;
 +#endif
 +
 +                                              if(!(ktom->attr&ATOM_ATTR_3BP))
 +                                                      continue;
 +
 +                                              //if(ktom==jtom)
 +                                              //      continue;
 +
 +                                              if(ktom==&(itom[i]))
 +                                                      continue;
 +
 +                                              moldyn->func_j0_k0(moldyn,
 +                                                                 &(itom[i]),
 +                                                                 jtom,
 +                                                                 ktom,
 +                                                                 bc_ik|bc_ij);
 +#ifdef STATIC_LISTS
                                        }
- #else
-                                       } while(list_next_f(that)!=\
-                                               L_NO_NEXT_ELEMENT);
- #endif
- 
-                               }
- 
-                               /* finish of first j loop after first k loop */
-                               if(moldyn->func_j0e)
-                                       moldyn->func_j0e(moldyn,
-                                                        &(itom[i]),
-                                                        jtom,
-                                                        bc_ij);
- 
  #ifdef STATIC_LISTS
                                }
+ #elif LOWMEM_LISTS
+                               }
  #else
                                } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
  #endif

diff --cc moldyn.h
Simple merge

diff --cc potentials/albe.c
index 9d8336e,1ec3938..9f98282
--- 1/potentials/albe.c
--- 2/potentials/albe.c
+++ b/potentials/albe.c
@@@ -271,88 -239,52 +280,66 @@@ int albe_mult_i0_j0_k0(t_moldyn *moldyn
                return 0;
        }
  
 -      /* d_ik */
 -      d_ik=sqrt(d_ik2);
 -
 -      /* dist_ij, d_ij */
 -      dist_ij=exchange->dist_ij;
 -      d_ij=exchange->d_ij;
 +      /* distances */
 +      distj=exchange->dist[j];
 +      distk=exchange->dist[k];
 +      dj=exchange->d[j];
 +      dk=exchange->d[k];
 +      djdk_inv=1.0/(dj*dk);
  
        /* cos theta */
 -      cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 -
 -      /* g_ijk */
 -      h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
 -      d2_h_cos2=exchange->di2+(h_cos*h_cos);
 -      frac=exchange->ci2/d2_h_cos2;
 -      g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
 -      dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
 -
 -      /* zeta sum += f_c_ik * g_ijk */
 -      if(d_ik<=R) {
 -              exchange->zeta_ij+=g;
 -              f_c_ik=1.0;
 -              df_c_ik=0.0;
 +      cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
 +
 +      /* g(cos(theta)) ij and ik values */
 +      h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism
 +      d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j);
 +      frac_j=*exchange->c2_[j]/d2_h_cos2_j;
 +      gj=1.0+*exchange->c2d2_[j]-frac_j;
 +      gj*=*(exchange->gamma_[j]);
 +      dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe 
 +      if(ak->brand==aj->brand) {
 +              gk=gj;
 +              dgk=dgj;
        }
        else {
 -              s_r=S-R;
 -              arg=M_PI*(d_ik-R)/s_r;
 -              f_c_ik=0.5+0.5*cos(arg);
 -              df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 -              exchange->zeta_ij+=f_c_ik*g;
 +              h_cos_k=*(exchange->h_[k])+cos_theta;
 +              d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k);
 +              frac_k=*exchange->c2_[k]/d2_h_cos2_k;
 +              gk=1.0+*exchange->c2d2_[k]-frac_k;
 +              gk*=*(exchange->gamma_[k]);
 +              dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
        }
  
-       /* zeta - for albe: ik depending g function */
- //if(ai->tag==0) {
- //    printf("------> %.15f %.15f\n",dj,dk);
- //    printf("------> %.15f %.15f\n",dj,dk);
- //}
- 
-       exchange->zeta[j]+=(exchange->f_c[k]*gk);
-       exchange->zeta[k]+=(exchange->f_c[j]*gj);
- 
-       /* cos theta derivatives */
-       v3_scale(&dcosdrj,&distk,djdk_inv);             // j
-       v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]);
-       v3_add(&dcosdrj,&dcosdrj,&tmp);
-       v3_scale(&dcosdrk,&distj,djdk_inv);             // k
-       v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]);
-       v3_add(&dcosdrk,&dcosdrk,&tmp);
+ #ifdef DEBUG
+       if(ai==&(moldyn->atom[DATOM])) 
+               printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik);
+ #endif
  
-       /* zeta derivatives */
-       dzjj=&(exchange->dzeta[j][j]);
-       dzkk=&(exchange->dzeta[k][k]);
-       dzjk=&(exchange->dzeta[j][k]);
-       dzkj=&(exchange->dzeta[k][j]);
-       v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
-       v3_add(dzjj,dzjj,&tmp);                         // j j
-       v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
-       v3_add(dzkk,dzkk,&tmp);                         // k k
-       v3_scale(&tmp,&distk,-exchange->df_c[k]*gk);    // dk rik = - di rik
-       v3_add(dzjk,dzjk,&tmp);
-       v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
-       v3_add(dzjk,dzjk,&tmp);                         // j k
-       v3_scale(&tmp,&distj,-exchange->df_c[j]*gj);    // dj rij = - di rij
-       v3_add(dzkj,dzkj,&tmp);
-       v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
-       v3_add(dzkj,dzkj,&tmp);                         // k j
+       /* store even more data for second k loop */
+       exchange->g[kcount]=g;
+       exchange->dg[kcount]=dg;
+       exchange->d_ik[kcount]=d_ik;
+       exchange->cos_theta[kcount]=cos_theta;
+       exchange->f_c_ik[kcount]=f_c_ik;
+       exchange->df_c_ik[kcount]=df_c_ik;
  
        /* increase k counter */
 -      exchange->kcount++;
 +      exchange->kcnt+=1;
 +              
 +      return 0;
 +}
 +
 +/* first j loop within first i loop */
 +int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 +
 +      t_albe_mult_params *params;
 +      t_albe_exchange *exchange;
 +
 +      params=moldyn->pot_params;
 +      exchange=&(params->exchange);
 +
 +      /* increase j counter */
 +      exchange->jcnt+=1;
  
        return 0;
  }
@@@ -440,16 -368,23 +427,26 @@@ printf("    t: %.15f %.15f %.15f\n",ai-
        v3_scale(&force,&force,-1.0); // dri rij = - drj rij
        v3_add(&(aj->f),&(aj->f),&force);
  
- #ifdef NDEBUG
- if(aj->tag==0) {
- printf("force: %.15f %.15f %.15f | %d %d (ji) %.15f\n",force.x,force.y,force.z,aj->tag,ai->tag,exchange->zeta[j]);
- printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
- }
+       /* virial */
+       virial_calc(ai,&force,&(exchange->dist_ij));
+ 
+ #ifdef DEBUG
+       if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
+               printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
+               printf("  adding %f %f %f\n",force.x,force.y,force.z);
+               if(ai==&(moldyn->atom[DATOM]))
+                       printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
+               if(aj==&(moldyn->atom[DATOM]))
+                       printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+               printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
+                                                   f_c,b,f_a,f_r);
+               printf("          %f %f %f\n",exchange->zeta_ij,.0,.0);
+       }
  #endif
  
 +      /* virial */
 +      virial_calc(ai,&force,dist);
 +
        /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
        exchange->pre_dzeta=0.5*f_a*f_c*db;
  
@@@ -526,32 -496,37 +555,50 @@@ int albe_mult_i0_j2_k0(t_moldyn *moldyn
        v3_scale(&force,&force,-1.0);
        v3_add(&(ai->f),&(ai->f),&force);
  
- #ifdef NDEBUG
- if(ai->tag==0) {
- printf("force: %.15f %.15f %.15f | %d %d %d -- %d(i contr k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag,k);
- printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
- printf("  ## %f\n",exchange->d[k]);
- }
+       /* derivative wrt k */
+       v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
+       v3_scale(&tmp,&dcosdrk,fcdg);
+       v3_add(&force,&force,&tmp);
+       v3_scale(&force,&force,pre_dzeta);
+ 
 -      /* force contribution */
 -      v3_add(&(ak->f),&(ak->f),&force);
++      v3_scale(&force,&force,-1.0);
++      v3_add(&(ai->f),&(ai->f),&force);
+ 
+ #ifdef DEBUG
+       if(ak==&(moldyn->atom[DATOM])) {
+               printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
+               printf("  adding %f %f %f\n",force.x,force.y,force.z);
+               printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
+               printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
+               printf("    d ij ik = %f %f\n",d_ij,d_ik);
+       }
  #endif
  
+       /* virial */
+       virial_calc(ai,&force,&dist_ik);
+       
+       /* force contribution to atom i */
+       v3_scale(&force,&force,-1.0);
+       v3_add(&(ai->f),&(ai->f),&force);
+ 
        /* increase k counter */
 -      exchange->kcount++;     
 +      exchange->kcnt+=1;
  
        return 0;
 +}
 +
 +int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
  
 +      t_albe_mult_params *params;
 +      t_albe_exchange *exchange;
 +
 +      params=moldyn->pot_params;
 +      exchange=&(params->exchange);
 +
 +      /* increase j counter */
 +      exchange->j2cnt+=1;
 +
 +      return 0;
  }
  
  int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {

diff --cc potentials/albe.h
index 39c4df8,89cdae2..9687655
--- 1/potentials/albe.h
--- 2/potentials/albe.h
+++ b/potentials/albe.h
@@@ -62,17 -65,17 +62,19 @@@ typedef struct s_albe_mult_params 
        double gamma[2];
        double gamma_m;
        double c[2];
 +      double c2[2];
        double c_mixed;
+       double c2[2];
        double c2_mixed;
        double d[2];
 +      double d2[2];
        double d_mixed;
+       double d2[2];
        double d2_mixed;
-       double c2d2[2];
-       double c2d2_m;
        double h[2];
        double h_mixed;
+       double c2d2[2];
+       double c2d2_m;
  
        t_albe_exchange exchange;       /* exchange between 2bp and 3bp calc */
  } t_albe_mult_params;