X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=faa6b7c3ea6354a7f1ab08a51cca462d567fd405;hb=2d562a1946322ae5b70bdce180321f32adbeab62;hp=f8bfd3f5face19085ac447cbcbfd9f5fcd4fbede;hpb=9c2172b5fce9688fa6be8d8a30745f9ef0b31419;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index f8bfd3f..faa6b7c 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -519,7 +519,8 @@ int scale_volume(t_moldyn *moldyn) {
 
 	t_atom *atom;
 	t_3dvec *dim,*vdim;
-	double virial,scale;
+	double scale,v;
+	t_virial virial;
 	t_linkcell *lc;
 	int i;
 
@@ -528,14 +529,25 @@ int scale_volume(t_moldyn *moldyn) {
 	vdim=&(moldyn->vis.dim);
 	lc=&(moldyn->lc);
 
-	for(i=0;i<moldyn->count;i++)
-		virial+=v3_norm(&(atom[i].virial));
+	memset(&virial,0,sizeof(t_virial));
+
+	for(i=0;i<moldyn->count;i++) {
+		virial.xx+=atom[i].virial.xx;
+		virial.yy+=atom[i].virial.yy;
+		virial.zz+=atom[i].virial.zz;
+		virial.xy+=atom[i].virial.xy;
+		virial.xz+=atom[i].virial.xz;
+		virial.yz+=atom[i].virial.yz;
+	}
+
+	/* just a guess so far ... */
+	v=virial.xx+virial.yy+virial.zz;
 
-printf("%f\n",virial);
+printf("%f\n",v);
 	/* get pressure from virial */
-	moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial;
+	moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+ONE_THIRD*v;
 	moldyn->p/=moldyn->volume;
-printf("%f\n",moldyn->p/(ATM));
+printf("%f | %f\n",moldyn->p/(ATM),moldyn->p_ref/ATM);
 
 	/* scale factor */
 	if(moldyn->pt_scale&P_SCALE_BERENDSEN)
@@ -767,23 +779,23 @@ int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
 	t_moldyn_schedule *schedule;
 
 	schedule=&(moldyn->schedule);
-	count=++(schedule->content_count);
+	count=++(schedule->total_sched);
 
-	ptr=realloc(moldyn->schedule.runs,count*sizeof(int));
+	ptr=realloc(schedule->runs,count*sizeof(int));
 	if(!ptr) {
 		perror("[moldyn] realloc (runs)");
 		return -1;
 	}
-	moldyn->schedule.runs=ptr;
-	moldyn->schedule.runs[count-1]=runs;
+	schedule->runs=ptr;
+	schedule->runs[count-1]=runs;
 
 	ptr=realloc(schedule->tau,count*sizeof(double));
 	if(!ptr) {
 		perror("[moldyn] realloc (tau)");
 		return -1;
 	}
-	moldyn->schedule.tau=ptr;
-	moldyn->schedule.tau[count-1]=tau;
+	schedule->tau=ptr;
+	schedule->tau[count-1]=tau;
 
 	return 0;
 }
@@ -806,16 +818,16 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) {
 
 int moldyn_integrate(t_moldyn *moldyn) {
 
-	int i,sched;
+	int i;
 	unsigned int e,m,s,v;
 	t_3dvec p;
-	t_moldyn_schedule *schedule;
+	t_moldyn_schedule *sched;
 	t_atom *atom;
 	int fd;
 	char dir[128];
 	double ds;
 
-	schedule=&(moldyn->schedule);
+	sched=&(moldyn->schedule);
 	atom=moldyn->atom;
 
 	/* initialize linked cell method */
@@ -852,12 +864,12 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	moldyn->debug=0;
 
 	/* executing the schedule */
-	for(sched=0;sched<moldyn->schedule.content_count;sched++) {
+	for(sched->count=0;sched->count<sched->total_sched;sched->count++) {
 
 		/* setting amount of runs and finite time step size */
-		moldyn->tau=schedule->tau[sched];
+		moldyn->tau=sched->tau[sched->count];
 		moldyn->tau_square=moldyn->tau*moldyn->tau;
-		moldyn->time_steps=schedule->runs[sched];
+		moldyn->time_steps=sched->runs[sched->count];
 
 	/* integration according to schedule */
 
@@ -907,7 +919,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 				visual_atoms(&(moldyn->vis),moldyn->time,
 				             moldyn->atom,moldyn->count);
 				printf("\rsched: %d, steps: %d, debug: %d",
-				       sched,i,moldyn->debug);
+				       sched->count,i,moldyn->debug);
 				fflush(stdout);
 			}
 		}
@@ -918,8 +930,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	}
 
 		/* check for hooks */
-		if(schedule->hook)
-			schedule->hook(moldyn,schedule->hook_params);
+		if(sched->hook)
+			sched->hook(moldyn,sched->hook_params);
 
 		/* get a new info line */
 		printf("\n");
@@ -984,6 +996,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 	int i,j,k,count;
 	t_atom *itom,*jtom,*ktom;
+	t_virial *virial;
 	t_linkcell *lc;
 	t_list neighbour_i[27];
 	t_list neighbour_i2[27];
@@ -1005,7 +1018,13 @@ int potential_force_calc(t_moldyn *moldyn) {
 		v3_zero(&(itom[i].f));
 
 		/* reset viral of atom i */
-		v3_zero(&(itom[i].virial));
+		virial=&(itom[i].virial);
+		virial->xx=0.0;
+		virial->yy=0.0;
+		virial->zz=0.0;
+		virial->xy=0.0;
+		virial->xz=0.0;
+		virial->yz=0.0;
 
 		/* reset site energy */
 		itom[i].e=0.0;
@@ -1108,6 +1127,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 		}
 
 	}
+#ifdef DEBUG
+printf("\n\n");
+#endif
 
 	return 0;
 }
@@ -1221,7 +1243,10 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
 
 	printf("[moldyn] tersoff parameter completion\n");
+	p->S2[0]=p->S[0]*p->S[0];
+	p->S2[1]=p->S[1]*p->S[1];
 	p->Smixed=sqrt(p->S[0]*p->S[1]);
+	p->S2mixed=p->Smixed*p->Smixed;
 	p->Rmixed=sqrt(p->R[0]*p->R[1]);
 	p->Amixed=sqrt(p->A[0]*p->A[1]);
 	p->Bmixed=sqrt(p->B[0]*p->B[1]);
@@ -1282,8 +1307,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	t_tersoff_mult_params *params;
 	t_tersoff_exchange *exchange;
 	t_3dvec dist_ij,force;
-	double d_ij;
-	double A,B,R,S,lambda,mu;
+	double d_ij,d_ij2;
+	double A,B,R,S,S2,lambda,mu;
 	double f_r,df_r;
 	double f_c,df_c;
 	int brand;
@@ -1314,18 +1339,10 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	 *
 	 */
 
-	/* dist_ij, d_ij */
-	v3_sub(&dist_ij,&(aj->r),&(ai->r));
-	if(bc) check_per_bound(moldyn,&dist_ij);
-	d_ij=v3_norm(&dist_ij);
-
-	/* save for use in 3bp */
-	exchange->d_ij=d_ij;
-	exchange->dist_ij=dist_ij;
-
 	/* constants */
 	if(brand==ai->brand) {
 		S=params->S[brand];
+		S2=params->S2[brand];
 		R=params->R[brand];
 		A=params->A[brand];
 		B=params->B[brand];
@@ -1335,6 +1352,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	}
 	else {
 		S=params->Smixed;
+		S2=params->S2mixed;
 		R=params->Rmixed;
 		A=params->Amixed;
 		B=params->Bmixed;
@@ -1343,10 +1361,23 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		params->exchange.chi=params->chi;
 	}
 
-	/* if d_ij > S => no force & potential energy contribution */
-	if(d_ij>S)
+	/* dist_ij, d_ij */
+	v3_sub(&dist_ij,&(aj->r),&(ai->r));
+	if(bc) check_per_bound(moldyn,&dist_ij);
+	d_ij2=v3_absolute_square(&dist_ij);
+
+	/* if d_ij2 > S2 => no force & potential energy contribution */
+	if(d_ij2>S2)
 		return 0;
 
+	/* now we will need the distance */
+	//d_ij=v3_norm(&dist_ij);
+	d_ij=sqrt(d_ij2);
+
+	/* save for use in 3bp */
+	exchange->d_ij=d_ij;
+	exchange->dist_ij=dist_ij;
+
 	/* more constants */
 	exchange->beta_j=&(params->beta[brand]);
 	exchange->n_j=&(params->n[brand]);
@@ -1386,7 +1417,6 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		s_r=S-R;
 		arg=M_PI*(d_ij-R)/s_r;
 		f_c=0.5+0.5*cos(arg);
-		//df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij));	/* MARK! */
 		df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 		/* two body contribution (ij, ji) */
 		v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
@@ -1398,6 +1428,23 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	 * dVij = dVji and we sum up both: no 1/2) */
 	v3_add(&(ai->f),&(ai->f),&force);
 
+	/* virial */
+	ai->virial.xx-=force.x*dist_ij.x;
+	ai->virial.yy-=force.y*dist_ij.y;
+	ai->virial.zz-=force.z*dist_ij.z;
+	ai->virial.xy-=force.x*dist_ij.y;
+	ai->virial.xz-=force.x*dist_ij.z;
+	ai->virial.yz-=force.y*dist_ij.z;
+
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+	printf("dVij, dVji (2bp) contrib:\n");
+	printf("%f | %f\n",force.x,ai->f.x);
+	printf("%f | %f\n",force.y,ai->f.y);
+	printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
+
 	/* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
 	moldyn->energy+=(0.5*f_r*f_c);
 
@@ -1489,6 +1536,23 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	/* add force */
 	v3_add(&(ai->f),&(ai->f),&force);
 
+	/* virial */
+	ai->virial.xx-=force.x*dist_ij->x;
+	ai->virial.yy-=force.y*dist_ij->y;
+	ai->virial.zz-=force.z*dist_ij->z;
+	ai->virial.xy-=force.x*dist_ij->y;
+	ai->virial.xz-=force.x*dist_ij->z;
+	ai->virial.yz-=force.y*dist_ij->z;
+
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+	printf("dVij (3bp) contrib:\n");
+	printf("%f | %f\n",force.x,ai->f.x);
+	printf("%f | %f\n",force.y,ai->f.y);
+	printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
+
 	/* add energy of 3bp sum */
 	moldyn->energy+=(0.5*f_c*b*f_a);
 
@@ -1517,6 +1581,23 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	/* add force */
 	v3_add(&(ai->f),&(ai->f),&force);
 
+	/* virial - plus sign, as dist_ij = - dist_ji - (really??) */
+	ai->virial.xx+=force.x*dist_ij->x;
+	ai->virial.yy+=force.y*dist_ij->y;
+	ai->virial.zz+=force.z*dist_ij->z;
+	ai->virial.xy+=force.x*dist_ij->y;
+	ai->virial.xz+=force.x*dist_ij->z;
+	ai->virial.yz+=force.y*dist_ij->z;
+
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+	printf("dVji (3bp) contrib:\n");
+	printf("%f | %f\n",force.x,ai->f.x);
+	printf("%f | %f\n",force.y,ai->f.y);
+	printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
+
 	return 0;
 }
 
@@ -1658,7 +1739,6 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 			s_r=S-R;
 			arg=M_PI*(d_ik-R)/s_r;
 			f_c_ik=0.5+0.5*cos(arg);
-			//df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); /* MARK */
 			df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 
 			/* zeta_ij */
@@ -1778,6 +1858,23 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 		v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
 		                                  /* scaled with 0.5 ^ */
 
+		/* virial */
+		ai->virial.xx-=temp2.x*dist_jk.x;
+		ai->virial.yy-=temp2.y*dist_jk.y;
+		ai->virial.zz-=temp2.z*dist_jk.z;
+		ai->virial.xy-=temp2.x*dist_jk.y;
+		ai->virial.xz-=temp2.x*dist_jk.z;
+		ai->virial.yz-=temp2.y*dist_jk.z;
+
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+	printf("dVjk (3bp) contrib:\n");
+	printf("%f | %f\n",temp2.x,ai->f.x);
+	printf("%f | %f\n",temp2.y,ai->f.y);
+	printf("%f | %f\n",temp2.z,ai->f.z);
+}
+#endif
+
 	}
 
 	return 0;