security checkin

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 6736e45..1f6391f 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -681,7 +681,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
        /* zero absolute time */
        moldyn->time=0.0;
 
-       /* debugging, ignre */
+       /* debugging, ignore */
        moldyn->debug=0;
 
        /* executing the schedule */
@@ -741,7 +741,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
                        if(!(i%v)) {
                                visual_atoms(&(moldyn->vis),moldyn->time,
                                             moldyn->atom,moldyn->count);
-                               printf("\rsched: %d, steps: %d, theta: %d",
+                               printf("\rsched: %d, steps: %d, debug: %d",
                                       sched,i,moldyn->debug);
                                fflush(stdout);
                        }
@@ -920,15 +920,14 @@ int potential_force_calc(t_moldyn *moldyn) {
                
                                /* 2bp post function */
                                if(moldyn->func2b_post) {
-printf("DEBUG: vor 2bp post\n");
                                        moldyn->func2b_post(moldyn,
                                                            &(itom[i]),
                                                            jtom,bc_ij);
-printf("DEBUG: nach 2bp post\n");
                                }
 
                        }
                }
+//printf("DEBUG: %.15f \n",itom[i].f.x);
        }
 
        return 0;
@@ -1119,6 +1118,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        /* clear 3bp and 2bp post run */
        exchange->run3bp=0;
        exchange->run2bp_post=0;
+
+       /* reset S > r > R mark */
+       exchange->d_ij_between_rs=0;
        
        /*
         * calc of 2bp contribution of V_ij and dV_ij/ji
@@ -1208,6 +1210,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij));
                /* two body contribution (ij, ji) */
                v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
+               /* tell 3bp that S > r > R */
+               exchange->d_ij_between_rs=1;
        }
 
        /* add forces of 2bp (ij, ji) contribution
@@ -1281,40 +1285,55 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        
        /* Vij and dVij */
        zeta=exchange->zeta_ij;
-       tmp=betaini*pow(zeta,ni-1.0);           /* beta^n * zeta^n-1 */
-       b=(1+zeta*tmp);                         /* 1 + beta^n * zeta^n */
-       db=chi*pow(b,-1.0/(2*ni)-1);            /* chi * (...)^(-1/2n - 1) */
-       b=db*b;                                 /* b_ij */
-       db*=-0.5*tmp;                           /* db_ij */
-       v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
-       v3_scale(&temp,dist_ij,df_a*b);
-       v3_add(&force,&force,&temp);
-       v3_scale(&force,&force,f_c);
-       v3_scale(&temp,dist_ij,df_c*b*f_a);
-       v3_add(&force,&force,&temp);
+       if(zeta==0.0) {
+               moldyn->debug++;                /* just for debugging ... */
+               db=0.0;
+               b=chi;
+       }
+       else {
+               tmp=betaini*pow(zeta,ni-1.0);           /* beta^n * zeta^n-1 */
+               b=(1+zeta*tmp);                         /* 1 + beta^n zeta^n */
+               db=chi*pow(b,-1.0/(2*ni)-1);            /* x(...)^(-1/2n - 1) */
+               b=db*b;                                 /* b_ij */
+               db*=-0.5*tmp;                           /* db_ij */
+               v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
+               v3_scale(&temp,dist_ij,df_a*b);
+               v3_add(&force,&force,&temp);
+               v3_scale(&force,&force,f_c);
+               v3_scale(&temp,dist_ij,df_c*b*f_a);
+               v3_add(&force,&force,&temp);
+               v3_scale(&force,&force,-0.5);
+
+               /* add force */
+               v3_add(&(ai->f),&(ai->f),&force);
+       }
 
        /* add energy of 3bp sum */
        moldyn->energy+=(0.5*f_c*b*f_a);
 
-       /* add force (sub, as F = - dVij) */
-       v3_sub(&(ai->f),&(ai->f),&force);
-
        /* dVji */
        zeta=exchange->zeta_ji;
-       tmp=betajnj*pow(zeta,nj-1.0);           /* beta^n * zeta^n-1 */
-       b=(1+zeta*tmp);                         /* 1 + beta^n * zeta^n */
-       db=chi*pow(b,-1.0/(2*nj)-1);            /* chi * (...)^(-1/2n - 1) */
-       b=db*b;                                 /* b_ij */
-       db*=-0.5*tmp;                           /* db_ij */
-       v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
-       v3_scale(&temp,dist_ij,df_a*b);
-       v3_add(&force,&force,&temp);
-       v3_scale(&force,&force,f_c);
-       v3_scale(&temp,dist_ij,df_c*b*f_a);
-       v3_add(&force,&force,&temp);
-
-       /* add force (sub, as F = - dVji) */
-       v3_sub(&(ai->f),&(ai->f),&force);
+       if(zeta==0.0) {
+               // do nothing ... (db=0, b=chi)
+               moldyn->debug++;
+       }
+       else {
+               tmp=betajnj*pow(zeta,nj-1.0);           /* beta^n * zeta^n-1 */
+               b=(1+zeta*tmp);                         /* 1 + beta^n zeta^n */
+               db=chi*pow(b,-1.0/(2*nj)-1);            /* x(...)^(-1/2n - 1) */
+               b=db*b;                                 /* b_ij */
+               db*=-0.5*tmp;                           /* db_ij */
+               v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
+               v3_scale(&temp,dist_ij,df_a*b);
+               v3_add(&force,&force,&temp);
+               v3_scale(&force,&force,f_c);
+               v3_scale(&temp,dist_ij,0.5*df_c*b*f_a);
+               v3_add(&force,&force,&temp);
+               v3_scale(&force,&force,-0.5);
+
+               /* add force */
+               v3_sub(&(ai->f),&(ai->f),&force);
+       }
 
        return 0;
 }
@@ -1569,8 +1588,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                /* betajnj * zeta_jk ^ nj-1 */
                tmp=exchange->betajnj*pow(zeta,(n-1.0));
                tmp=-chi/2.0*pow(1+tmp*zeta,-1.0/(2.0*n)-1)*tmp;
-               v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk);
+               v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
                v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */
+                                                 /* scaled with 0.5 ^ */
        }
 
        return 0;