X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=6af2dff0b525cd3e96789d2b93fd740e8b043694;hb=9d3e2be6625c7076a7b9243d922195c499e0bf4d;hp=716a087fa90f1b89b073411af22e2d55995aba8a;hpb=15dd7ee084a5db24a3e1a15428cd47e8d61fda5b;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 716a087..6af2dff 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -141,7 +141,7 @@ int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) {
 
 int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) {
 
-	moldyn->func2b=func;
+	moldyn->func2b_post=func;
 	moldyn->pot2b_params=params;
 
 	return 0;
@@ -784,7 +784,6 @@ int velocity_verlet(t_moldyn *moldyn) {
 
 	/* forces depending on chosen potential */
 	potential_force_calc(moldyn);
-	//moldyn->potential_force_function(moldyn);
 
 	for(i=0;i<count;i++) {
 		/* again velocities */
@@ -873,15 +872,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 					   !(jtom->attr&ATOM_ATTR_3BP))
 						continue;
 
-			/*
-			 * according to mr. nordlund, we dont need to take the 
-			 * sum over all atoms now, as 'this is centered' around
-			 * atom i ...
-			 * i am not quite sure though! there is a not vanishing
-			 * part even if f_c_ik is zero ...
-			 * this analytical potentials suck!
-			 * switching from mc to md to dft soon!
-			 */
+			/* neighbourhood of atom j is not needed! */
 
 			//		link_cell_neighbour_index(moldyn,
 			//		   (jtom->r.x+moldyn->dim.x/2)/lc->x,
@@ -948,9 +939,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 			if(ktom==&(itom[i]))
 				continue;
 
-//printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
 			moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
-//printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
 
 						} while(list_next(that)!=\
 						        L_NO_NEXT_ELEMENT);
@@ -960,10 +949,13 @@ int potential_force_calc(t_moldyn *moldyn) {
 				} while(list_next(this)!=L_NO_NEXT_ELEMENT);
 		
 				/* 2bp post function */
-				if(moldyn->func2b_post)
+				if(moldyn->func2b_post) {
+printf("before 2bp post: %.15f %.15f %.15f\n",itom[i].f.x,itom[i].r.x,itom[i].v.x);
 					moldyn->func2b_post(moldyn,
 					                    &(itom[i]),
 					                    jtom,bc_ij);
+printf("after 2bp post: %.15f %.15f %.15f\n",itom[i].f.x,itom[i].r.x,itom[i].v.x);
+				}
 
 			}
 		}
@@ -1156,6 +1148,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	exchange=&(params->exchange);
 
 	exchange->run3bp=0;
+	exchange->run2bp_post=0;
 	
 	/*
 	 * we need: f_c, df_c, f_r, df_r
@@ -1193,14 +1186,13 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		mu=params->mu_m;
 		params->exchange.chi=params->chi;
 	}
-
 	if(d_ij>S)
 		return 0;
 
 	f_r=A*exp(-lambda*d_ij);
 	df_r=-lambda*f_r/d_ij;
 
-	/* f_a, df_a calc + save for 3bp use */
+	/* f_a, df_a calc + save for later use */
 	exchange->f_a=-B*exp(-mu*d_ij);
 	exchange->df_a=-mu*exchange->f_a/d_ij;
 
@@ -1228,8 +1220,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	exchange->f_c=f_c;
 	exchange->df_c=df_c;
 
-	/* enable the run of 3bp function */
+	/* enable the run of 3bp function and 2bp post processing */
 	exchange->run3bp=1;
+	exchange->run2bp_post=1;
 
 	/* reset 3bp sums */
 	exchange->sum1_3bp=0.0;
@@ -1259,6 +1252,10 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	params=moldyn->pot2b_params;
 	exchange=&(params->exchange);
 
+	/* we do not run if f_c_ij was dtected to be 0! */
+	if(!(exchange->run2bp_post))
+		return 0;
+
 	db_ij=&(exchange->db_ij);
 	f_c=exchange->f_c;
 	df_c=exchange->df_c;
@@ -1266,6 +1263,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	df_a=exchange->df_a;
 	betan=exchange->betan;
 	n=*(exchange->n);
+	chi=exchange->chi;
 	dist_ij=&(exchange->dist_ij);
 
 	db_ij_scale1=(1+betan*exchange->sum1_3bp);
@@ -1274,19 +1272,24 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	b_ij=chi*db_ij_scale1*help;
 	db_ij_scale1=-chi/(2*n)*help;
 
+	/* db_ij part */
 	v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
 	v3_scale(db_ij,db_ij,f_a);
 
+	/* df_a part */
 	v3_scale(&temp,dist_ij,b_ij*df_a);
 
+	/* db_ij + df_a part */
 	v3_add(&force,&temp,db_ij);
 	v3_scale(&force,&force,f_c);
 
+	/* df_c part */
 	v3_scale(&temp,dist_ij,f_a*b_ij*df_c);
 
 	/* add energy of 3bp sum */
 	moldyn->energy+=(0.5*f_c*b_ij*f_a);
-	/* add force of 3bp calculation */	
+
+	/* add force of 3bp calculation (all three parts) */
 	v3_add(&(ai->f),&temp,&force);
 
 	return 0;
@@ -1415,15 +1418,12 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 		bracket=0.0;
 		bracket_n_1=0.0;
 		bracket_n=0.0;
-		printf("Foo -> 0: ");
 	}
 	else {
 		bracket=f_c_ik*g;
 		bracket_n_1=pow(bracket,n-1.0);
 		bracket_n=bracket_n_1*bracket;
-		printf("Foo -> 1: ");
 	}
-//printf("%.15f %.15f %.15f\n",bracket_n_1,bracket_n,bracket);
 
 	/* calc of db_ij and the 2 sums */
 	exchange->sum1_3bp+=bracket_n;