X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=moldyn.c;h=7496223be1207d14f86ccc8d7010befdb9d8b725;hb=da2d9866e05b1b7a408ecda2b1695e07c30b0533;hp=e2c84a524928a9ab12f7b18cf53cdbca18843208;hpb=19bf7f2ba36f79406db4efbf1201481a846762ac;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index e2c84a5..7496223 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -1738,6 +1738,9 @@ int velocity_verlet(t_moldyn *moldyn) {
 	tau_square=moldyn->tau_square;
 
 	for(i=0;i<count;i++) {
+		/* check whether fixed atom */
+		if(atom[i].attr&ATOM_ATTR_FP)
+			continue;
 		/* new positions */
 		h=0.5/atom[i].mass;
 		v3_scale(&delta,&(atom[i].v),tau);
@@ -1761,6 +1764,9 @@ int velocity_verlet(t_moldyn *moldyn) {
 	potential_force_calc(moldyn);
 
 	for(i=0;i<count;i++) {
+		/* check whether fixed atom */
+		if(atom[i].attr&ATOM_ATTR_FP)
+			continue;
 		/* again velocities [actually v(t+tau)] */
 		v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 		v3_add(&(atom[i].v),&(atom[i].v),&delta);
@@ -2088,9 +2094,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 	//printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
 	if(moldyn->time>DSTART&&moldyn->time<DEND) {
 		printf("force:\n");
-		printf("  x: %0.40f\n",moldyn->atom[5832].f.x);
-		printf("  y: %0.40f\n",moldyn->atom[5832].f.y);
-		printf("  z: %0.40f\n",moldyn->atom[5832].f.z);
+		printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
+		printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
+		printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
 	}
 #endif
 
@@ -2316,6 +2322,47 @@ int pair_correlation_init(t_moldyn *moldyn,double dr) {
 	return 0;
 }
 
+int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
+
+	int i;
+	t_atom *atom;
+	t_3dvec dist;
+	double d2;
+	int a_cnt;
+	int b_cnt;
+
+	atom=moldyn->atom;
+	dc[0]=0;
+	dc[1]=0;
+	dc[2]=0;
+	a_cnt=0;
+	b_cnt=0;
+
+	for(i=0;i<moldyn->count;i++) {
+
+		v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
+		check_per_bound(moldyn,&dist);
+		d2=v3_absolute_square(&dist);
+
+		if(atom[i].brand) {
+			b_cnt+=1;
+			dc[1]+=d2;
+		}
+		else {
+			a_cnt+=1;
+			dc[0]+=d2;
+		}
+
+		dc[2]+=d2;
+	}
+
+	dc[0]*=(1.0/(6.0*moldyn->time*a_cnt));
+	dc[1]*=(1.0/(6.0*moldyn->time*b_cnt));
+	dc[2]*=(1.0/(6.0*moldyn->time*moldyn->count));
+		
+	return 0;
+}
+
 int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 
 	int slots;
@@ -2333,13 +2380,13 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 	unsigned char bc;
 	t_3dvec dist;
 	double d;
-	//double norm;
+	double norm;
 	int o,s;
 	unsigned char ibrand;
 
 	lc=&(moldyn->lc);
 
-	slots=moldyn->cutoff/dr;
+	slots=2.0*moldyn->cutoff/dr;
 	o=2*slots;
 
 	if(slots*dr<=moldyn->cutoff)
@@ -2416,7 +2463,7 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 				d=v3_absolute_square(&dist);
 
 				/* ignore if greater or equal cutoff */
-				if(d>=moldyn->cutoff_square)
+				if(d>=4.0*moldyn->cutoff_square)
 					continue;
 
 				/* fill the slots */
@@ -2452,16 +2499,17 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 		}
 	}
 
-	/* normalization 
+	/* normalization */
 	for(i=1;i<slots;i++) {
-		 // normalization: 4 pi r r dr
+		 // normalization: 4 pi r^2 dr
 		 // here: not double counting pairs -> 2 pi r r dr
-		norm=2*M_PI*moldyn->count*(i*dr*i*dr)*dr;
+		 // ... and actually it's a constant times r^2
+		norm=i*i*dr*dr;
 		stat[i]/=norm;
 		stat[slots+i]/=norm;
 		stat[o+i]/=norm;
 	}
-	*/
+	/* */
 
 	if(ptr==NULL) {
 		/* todo: store/print pair correlation function */