added MATTONI force field (left parameters by now)

diff --git a/potentials/albe_fast.c b/potentials/albe_fast.c
index 0db5f66..9fc362e 100644 (file)
--- a/potentials/albe_fast.c
+++ b/potentials/albe_fast.c
@@ -452,7 +452,11 @@ int albe_potential_force_calc(t_moldyn *moldyn) {
        }
 
        /* force contribution for atom i */
+#ifdef MATTONI
+       scale=-0.5*(f_c*(df_r-b*df_a)); // - in albe formalism
+#else
        scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
+#endif
        v3_scale(&force,&(dist_ij),scale);
        v3_add(&(ai->f),&(ai->f),&force);
 
@@ -604,10 +608,15 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
        v3_add(&(ai->f),&(ai->f),&force);
 
        /* derivative wrt k */
+#ifdef MATTONI
+       v3_scale(&tmp,&dcosdrk,fcdg);
+       v3_scale(&force,&tmp,pre_dzeta);
+#else
        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);
+#endif
 
        /* force contribution */
        v3_add(&(ktom->f),&(ktom->f),&force);
@@ -1083,7 +1092,11 @@ void *potential_force_thread(void *ptr) {
        }
 
        /* force contribution for atom i */
+#ifdef MATTONI
+       scale=-0.5*(f_c*(df_r-b*df_a)); // - in albe formalism
+#else
        scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
+#endif
        v3_scale(&force,&(dist_ij),scale);
        pthread_mutex_lock(&(amutex[ai->tag])); 
        v3_add(&(ai->f),&(ai->f),&force);
@@ -1249,10 +1262,15 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
        pthread_mutex_unlock(&(amutex[ai->tag]));       
 
        /* derivative wrt k */
+#ifdef MATTONI
+       v3_scale(&tmp,&dcosdrk,fcdg);
+       v3_scale(&force,&tmp,pre_dzeta);
+#else
        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);
+#endif
 
        /* force contribution */
        pthread_mutex_lock(&(amutex[ktom->tag]));