X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=12cb49b4f7f1cb4fb8d8da5a80727f085dbdad4f;hb=0d2f9a11030dff3583104dac5d4dcb9f040a1327;hp=7698c24767b90de44a165602077a46714d86d476;hpb=bbc759243283faace6edb73684f963d25c820e6b;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 7698c24..12cb49b 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -258,7 +258,7 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 
 	switch(type) {
 		case MOLDYN_POTENTIAL_TM:
-			moldyn->func1b=tersoff_mult_1bp;
+			//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;
@@ -802,6 +802,8 @@ int del_atom(t_moldyn *moldyn,int tag) {
 			case DEFECT_STYPE_DB_Z:\
 				d_o.z=d_params->od;\
 				d_d.z=d_params->dd;\
+d_d.x=0.9;\
+d_d.y=0.9;\
 				break;\
 			case DEFECT_STYPE_DB_R:\
 				break;\
@@ -2071,11 +2073,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		temperature_calc(moldyn);
 		virial_sum(moldyn);
 		pressure_calc(moldyn);
-		/*
+#ifdef PDEBUG	
 		thermodynamic_pressure_calc(moldyn);
 		printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
 		       moldyn->tp/BAR);
-		*/
+#endif
 
 		/* calculate fluctuations + averages */
 		average_and_fluctuation_calc(moldyn);
@@ -2182,15 +2184,20 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		}
 
 		/* display progress */
+#ifndef PDEBUG
 		if(!(i%10)) {
+#endif
 			/* get current time */
 			gettimeofday(&t2,NULL);
 
 printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
+#ifndef PDEBUG
        moldyn->p/BAR,moldyn->p_avg/BAR,
-       //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#else
+       moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#endif
        moldyn->volume,
        (int)(t2.tv_sec-t1.tv_sec));
 
@@ -2198,7 +2205,9 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
 
 			/* copy over time */
 			t1=t2;
+#ifndef PDEBUG
 		}
+#endif
 
 		/* increase absolute time */
 		moldyn->time+=moldyn->tau;
@@ -2246,6 +2255,7 @@ int velocity_verlet(t_moldyn *moldyn) {
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
 		v3_scale(&delta,&(atom[i].f),h*tau_square);
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
+		//check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
 		check_per_bound(moldyn,&(atom[i].r));
 
 		/* velocities [actually v(t+tau/2)] */
@@ -2261,7 +2271,11 @@ int velocity_verlet(t_moldyn *moldyn) {
 
 	/* forces depending on chosen potential */
 #ifndef ALBE_FAST
-	potential_force_calc(moldyn);
+	// if albe, use albe force calc routine
+	//if(moldyn->func3b_j1==albe_mult_3bp_j1)
+	//	albe_potential_force_calc(moldyn);
+	//else
+		potential_force_calc(moldyn);
 #else
 	albe_potential_force_calc(moldyn);
 #endif
@@ -2713,6 +2727,51 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
 	return 0;
 }
         
+int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a) {
+	
+	double x,y,z;
+	t_3dvec *dim;
+
+	dim=&(moldyn->dim);
+
+	x=dim->x/2;
+	y=dim->y/2;
+	z=dim->z/2;
+
+	if(moldyn->status&MOLDYN_STAT_PBX) {
+		if(a->r.x>=x) {
+			a->pbc[0]+=1;
+			a->r.x-=dim->x;
+		}
+		else if(-a->r.x>x) {
+			a->pbc[0]-=1;
+			a->r.x+=dim->x;
+		}
+	}
+	if(moldyn->status&MOLDYN_STAT_PBY) {
+		if(a->r.y>=y) {
+			a->pbc[1]+=1;
+			a->r.y-=dim->y;
+		}
+		else if(-a->r.y>y) {
+			a->pbc[1]-=1;
+			a->r.y+=dim->y;
+		}
+	}
+	if(moldyn->status&MOLDYN_STAT_PBZ) {
+		if(a->r.z>=z) {
+			a->pbc[2]+=1;
+			a->r.z-=dim->z;
+		}
+		else if(-a->r.z>z) {
+			a->pbc[2]-=1;
+			a->r.z+=dim->z;
+		}
+	}
+
+	return 0;
+}
+        
 /*
  * debugging / critical check functions
  */
@@ -3048,6 +3107,7 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 	int i;
 	t_atom *atom;
 	t_3dvec dist;
+	t_3dvec final_r;
 	double d2;
 	int a_cnt;
 	int b_cnt;
@@ -3061,8 +3121,12 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 
 	for(i=0;i<moldyn->count;i++) {
 
-		v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
-		check_per_bound(moldyn,&dist);
+		/* care for pb crossing */
+		final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x;
+		final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y;
+		final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z;
+		/* calculate distance */
+		v3_sub(&dist,&final_r,&(atom[i].r_0));
 		d2=v3_absolute_square(&dist);
 
 		if(atom[i].brand) {