scaling to fs, angstrom, amu done, probs with velocity scaling!

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 1f6391f..9d9901c 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -385,12 +385,17 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
        else
                if(moldyn->pt_scale&T_SCALE_BERENDSEN)
                        scale=1.0+moldyn->tau*(scale-1.0)/moldyn->t_tc;
+printf("scale=%f\n",scale);
        scale=sqrt(scale);
+printf("debug: %f %f %f %f \n",scale,moldyn->t_ref,moldyn->t,moldyn->tau);
 
        /* velocity scaling */
-       for(i=0;i<moldyn->count;i++)
+       for(i=0;i<moldyn->count;i++) {
+printf("vorher: %f\n",atom[i].v.x);
                if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
                        v3_scale(&(atom[i].v),&(atom[i].v),scale);
+printf("nachher: %f\n",atom[i].v.x);
+       }
 
        return 0;
 }
@@ -447,12 +452,6 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
 
        /* nn_dist is the nearest neighbour distance */
 
-       if(moldyn->t==5.0) {
-               printf("[moldyn] i do not estimate timesteps below %f K!\n",
-                      MOLDYN_CRITICAL_EST_TEMP);
-               return 23.42;
-       }
-
        tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
 
        return tau;     
@@ -813,6 +812,7 @@ int velocity_verlet(t_moldyn *moldyn) {
 /* generic potential and force calculation */
 
 int potential_force_calc(t_moldyn *moldyn) {
+printf("start pot force calc\n");
 
        int i,j,k,count;
        t_atom *itom,*jtom,*ktom;
@@ -833,6 +833,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 
        /* get energy and force of every atom */
        for(i=0;i<count;i++) {
+printf("atom %d: %f\n",i,itom[i].r.x);
+printf("atom %d: %f\n",i,itom[i].v.x);
+printf("atom %d: %f\n",i,itom[i].f.x);
 
                /* reset force */
                v3_zero(&(itom[i].f));
@@ -927,9 +930,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                        }
                }
-//printf("DEBUG: %.15f \n",itom[i].f.x);
        }
 
+printf("end pot force calc\n");
        return 0;
 }
 
@@ -1050,8 +1053,8 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
        p->mu_m=0.5*(p->mu[0]+p->mu[1]);
 
        printf("[moldyn] tersoff mult parameter info:\n");
-       printf("  S (m)  | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed);
-       printf("  R (m)  | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed);
+       printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
+       printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
        printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
        printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
        printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
@@ -1189,7 +1192,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */
        f_r=A*exp(-lambda*d_ij);
-       df_r=-lambda*f_r/d_ij;
+       df_r=lambda*f_r/d_ij;
 
        /* f_a, df_a calc (again, same for ij and ji) | save for later use! */
        exchange->f_a=-B*exp(-mu*d_ij);
@@ -1289,6 +1292,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                moldyn->debug++;                /* just for debugging ... */
                db=0.0;
                b=chi;
+               v3_scale(&force,dist_ij,df_a*b*f_c);
        }
        else {
                tmp=betaini*pow(zeta,ni-1.0);           /* beta^n * zeta^n-1 */
@@ -1300,13 +1304,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                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);
        }
+       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);
@@ -1314,8 +1318,9 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        /* dVji */
        zeta=exchange->zeta_ji;
        if(zeta==0.0) {
-               // do nothing ... (db=0, b=chi)
                moldyn->debug++;
+               b=chi;
+               v3_scale(&force,dist_ij,df_a*b*f_c);
        }
        else {
                tmp=betajnj*pow(zeta,nj-1.0);           /* beta^n * zeta^n-1 */
@@ -1327,13 +1332,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                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);
        }
+       v3_scale(&temp,dist_ij,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;
 }
@@ -1596,3 +1601,32 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        return 0;
 }
 
+
+/*
+ * debugging / critical check functions
+ */
+
+int moldyn_bc_check(t_moldyn *moldyn) {
+
+       t_atom *atom;
+       t_3dvec *dim;
+       int i;
+
+       atom=moldyn->atom;
+       dim=&(moldyn->dim);
+
+       for(i=0;i<moldyn->count;i++) {
+               if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
+                       printf("FATAL: atom %d: x: %.20f (%.20f)\n",
+                              i,atom[i].r.x*1e10,dim->x/2*1e10);
+               if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
+                       printf("FATAL: atom %d: y: %.20f (%.20f)\n",
+                              i,atom[i].r.y*1e10,dim->y/2*1e10);
+               if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
+                       printf("FATAL: atom %d: z: %.20f (%.20f)\n",
+                              i,atom[i].r.z*1e10,dim->z/2*1e10);
+       }
+
+       return 0;
+}
+