X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=moldyn.c;h=93c72daa1dccc633c766794a777446d20aacf164;hb=9b3b22db2ac62d0fc45a072749bdb4ad70d10f99;hp=94beb044eacc2211dc290e7e12729de4a7b18ec8;hpb=b3f72740276662c09ff43c294b04409a5bb2e013;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 94beb04..93c72da 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
 
@@ -2223,6 +2229,8 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 
 	int fd;
 	int cnt,size;
+	int fsize;
+	int corr;
 
 	fd=open(file,O_RDONLY);
 	if(fd<0) {
@@ -2230,6 +2238,9 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 		return fd;
 	}
 
+	fsize=lseek(fd,0,SEEK_END);
+	lseek(fd,0,SEEK_SET);
+
 	size=sizeof(t_moldyn);
 
 	while(size) {
@@ -2243,6 +2254,19 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 
 	size=moldyn->count*sizeof(t_atom);
 
+	/* correcting possible atom data offset */
+	corr=0;
+	if(fsize!=sizeof(t_moldyn)+size) {
+		corr=fsize-sizeof(t_moldyn)-size;
+		printf("[moldyn] WARNING: lsf (illegal file size)\n");
+		printf("  moifying offset:\n");
+		printf("  - current pos: %d\n",sizeof(t_moldyn));
+		printf("  - atom size: %d\n",size);
+		printf("  - file size: %d\n",fsize);
+		printf("  => correction: %d\n",corr);
+		lseek(fd,corr,SEEK_CUR);
+	}
+
 	moldyn->atom=(t_atom *)malloc(size);
 	if(moldyn->atom==NULL) {
 		perror("[moldyn] load save file malloc (atoms)");
@@ -2298,6 +2322,42 @@ 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[2]+=d2;
+	}	
+		
+	return 0;
+}
+
 int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 
 	int slots;