X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=5aa6d87bb4490ad4585d9a196ec6b4d3867fc511;hb=c0a8b254109929fba10795e644187c51742108a8;hp=c25fed3ae9d4c21587336bc2cf321d0bd96b5bbb;hpb=785f05e0c0ddb57428dc71b0a3f7798673799871;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index c25fed3..5aa6d87 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -230,33 +230,38 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
  */
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
-                   u8 attr,u8 bnum,int a,int b,int c) {
+                   u8 attr,u8 brand,int a,int b,int c) {
 
-	int count;
+	int new,count;
 	int ret;
 	t_3dvec origin;
+	void *ptr;
+	t_atom *atom;
 
-	count=a*b*c;
+	new=a*b*c;
+	count=moldyn->count;
 
 	/* how many atoms do we expect */
-	if(type==FCC) count*=4;
-	if(type==DIAMOND) count*=8;
+	if(type==FCC) new*=4;
+	if(type==DIAMOND) new*=8;
 
 	/* allocate space for atoms */
-	moldyn->atom=malloc(count*sizeof(t_atom));
-	if(moldyn->atom==NULL) {
-		perror("malloc (atoms)");
+	ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
+	if(!ptr) {
+		perror("[moldyn] realloc (create lattice)");
 		return -1;
 	}
-
+	moldyn->atom=ptr;
+	atom=&(moldyn->atom[count]);
+		
 	v3_zero(&origin);
 
 	switch(type) {
 		case FCC:
-			ret=fcc_init(a,b,c,lc,moldyn->atom,&origin);
+			ret=fcc_init(a,b,c,lc,atom,&origin);
 			break;
 		case DIAMOND:
-			ret=diamond_init(a,b,c,lc,moldyn->atom,&origin);
+			ret=diamond_init(a,b,c,lc,atom,&origin);
 			break;
 		default:
 			printf("unknown lattice type (%02x)\n",type);
@@ -264,24 +269,24 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	}
 
 	/* debug */
-	if(ret!=count) {
+	if(ret!=new) {
 		printf("[moldyn] creating lattice failed\n");
 		printf("  amount of atoms\n");
-		printf("  - expected: %d\n",count);
+		printf("  - expected: %d\n",new);
 		printf("  - created: %d\n",ret);
 		return -1;
 	}
 
-	moldyn->count=count;
-	printf("[moldyn] created lattice with %d atoms\n",count);
+	moldyn->count+=new;
+	printf("[moldyn] created lattice with %d atoms\n",new);
 
-	while(count) {
-		count-=1;
-		moldyn->atom[count].element=element;
-		moldyn->atom[count].mass=mass;
-		moldyn->atom[count].attr=attr;
-		moldyn->atom[count].bnum=bnum;
-		check_per_bound(moldyn,&(moldyn->atom[count].r));
+	for(ret=0;ret<new;ret++) {
+		atom[ret].element=element;
+		atom[ret].mass=mass;
+		atom[ret].attr=attr;
+		atom[ret].brand=brand;
+		atom[ret].tag=count+ret;
+		check_per_bound(moldyn,&(atom[ret].r));
 	}
 
 	return ret;
@@ -371,7 +376,7 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 	return count;
 }
 
-int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
+int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
              t_3dvec *r,t_3dvec *v) {
 
 	t_atom *atom;
@@ -379,9 +384,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
 	int count;
 	
 	atom=moldyn->atom;
-	count=++(moldyn->count);
+	count=(moldyn->count)++;
 
-	ptr=realloc(atom,count*sizeof(t_atom));
+	ptr=realloc(atom,(count+1)*sizeof(t_atom));
 	if(!ptr) {
 		perror("[moldyn] realloc (add atom)");
 		return -1;
@@ -389,12 +394,13 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
 	moldyn->atom=ptr;
 
 	atom=moldyn->atom;
-	atom[count-1].r=*r;
-	atom[count-1].v=*v;
-	atom[count-1].element=element;
-	atom[count-1].mass=mass;
-	atom[count-1].bnum=bnum;
-	atom[count-1].attr=attr;
+	atom[count].r=*r;
+	atom[count].v=*v;
+	atom[count].element=element;
+	atom[count].mass=mass;
+	atom[count].brand=brand;
+	atom[count].tag=count;
+	atom[count].attr=attr;
 
 	return 0;
 }
@@ -761,23 +767,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;
 }
@@ -800,16 +806,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 */
@@ -846,12 +852,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 */
 
@@ -901,7 +907,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);
 			}
 		}
@@ -912,8 +918,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");
@@ -1001,6 +1007,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 		/* reset viral of atom i */
 		v3_zero(&(itom[i].virial));
 
+		/* reset site energy */
+		itom[i].e=0.0;
+
 		/* single particle potential/force */
 		if(itom[i].attr&ATOM_ATTR_1BP)
 			moldyn->func1b(moldyn,&(itom[i]));
@@ -1099,6 +1108,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 		}
 
 	}
+#ifdef DEBUG
+printf("\n\n");
+#endif
 
 	return 0;
 }
@@ -1240,11 +1252,11 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
 /* tersoff 1 body part */
 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
 
-	int num;
+	int brand;
 	t_tersoff_mult_params *params;
 	t_tersoff_exchange *exchange;
 	
-	num=ai->bnum;
+	brand=ai->brand;
 	params=moldyn->pot1b_params;
 	exchange=&(params->exchange);
 
@@ -1253,15 +1265,15 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
 	 *         their right values
 	 */
 
-	exchange->beta_i=&(params->beta[num]);
-	exchange->n_i=&(params->n[num]);
-	exchange->c_i=&(params->c[num]);
-	exchange->d_i=&(params->d[num]);
-	exchange->h_i=&(params->h[num]);
+	exchange->beta_i=&(params->beta[brand]);
+	exchange->n_i=&(params->n[brand]);
+	exchange->c_i=&(params->c[brand]);
+	exchange->d_i=&(params->d[brand]);
+	exchange->h_i=&(params->h[brand]);
 
 	exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
-	exchange->ci2=params->c[num]*params->c[num];
-	exchange->di2=params->d[num]*params->d[num];
+	exchange->ci2=params->c[brand]*params->c[brand];
+	exchange->di2=params->d[brand]*params->d[brand];
 	exchange->ci2di2=exchange->ci2/exchange->di2;
 
 	return 0;
@@ -1277,12 +1289,12 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	double A,B,R,S,lambda,mu;
 	double f_r,df_r;
 	double f_c,df_c;
-	int num;
+	int brand;
 	double s_r;
 	double arg;
 
 	params=moldyn->pot2b_params;
-	num=aj->bnum;
+	brand=aj->brand;
 	exchange=&(params->exchange);
 
 	/* clear 3bp and 2bp post run */
@@ -1315,13 +1327,13 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	exchange->dist_ij=dist_ij;
 
 	/* constants */
-	if(num==ai->bnum) {
-		S=params->S[num];
-		R=params->R[num];
-		A=params->A[num];
-		B=params->B[num];
-		lambda=params->lambda[num];
-		mu=params->mu[num];
+	if(brand==ai->brand) {
+		S=params->S[brand];
+		R=params->R[brand];
+		A=params->A[brand];
+		B=params->B[brand];
+		lambda=params->lambda[brand];
+		mu=params->mu[brand];
 		exchange->chi=1.0;
 	}
 	else {
@@ -1339,12 +1351,12 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		return 0;
 
 	/* more constants */
-	exchange->beta_j=&(params->beta[num]);
-	exchange->n_j=&(params->n[num]);
-	exchange->c_j=&(params->c[num]);
-	exchange->d_j=&(params->d[num]);
-	exchange->h_j=&(params->h[num]);
-	if(num==ai->bnum) {
+	exchange->beta_j=&(params->beta[brand]);
+	exchange->n_j=&(params->n[brand]);
+	exchange->c_j=&(params->c[brand]);
+	exchange->d_j=&(params->d[brand]);
+	exchange->h_j=&(params->h[brand]);
+	if(brand==ai->brand) {
 		exchange->betajnj=exchange->betaini;
 		exchange->cj2=exchange->ci2;
 		exchange->dj2=exchange->di2;
@@ -1352,8 +1364,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	}
 	else {
 		exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
-		exchange->cj2=params->c[num]*params->c[num];
-		exchange->dj2=params->d[num]*params->d[num];
+		exchange->cj2=params->c[brand]*params->c[brand];
+		exchange->dj2=params->d[brand]*params->d[brand];
 		exchange->cj2dj2=exchange->cj2/exchange->dj2;
 	}
 
@@ -1388,6 +1400,14 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	/* add forces of 2bp (ij, ji) contribution
 	 * dVij = dVji and we sum up both: no 1/2) */
 	v3_add(&(ai->f),&(ai->f),&force);
+#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);
@@ -1479,6 +1499,14 @@ 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);
+#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);
@@ -1507,6 +1535,14 @@ 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);
+#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;
 }
@@ -1533,7 +1569,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	double h_cos,d2_h_cos2;
 	double frac,g,zeta,chi;
 	double tmp;
-	int num;
+	int brand;
 
 	params=moldyn->pot3b_params;
 	exchange=&(params->exchange);
@@ -1585,10 +1621,10 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	d_ik=v3_norm(&dist_ik);
 
 	/* ik constants */
-	num=ai->bnum;
-	if(num==ak->bnum) {
-		R=params->R[num];
-		S=params->S[num];
+	brand=ai->brand;
+	if(brand==ak->brand) {
+		R=params->R[brand];
+		S=params->S[brand];
 	}
 	else {
 		R=params->Rmixed;
@@ -1671,12 +1707,12 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	d_jk=v3_norm(&dist_jk);
 
 	/* jk constants */
-	num=aj->bnum;
-	if(num==ak->bnum) {
-		R=params->R[num];
-		S=params->S[num];
-		B=params->B[num];
-		mu=params->mu[num];
+	brand=aj->brand;
+	if(brand==ak->brand) {
+		R=params->R[brand];
+		S=params->S[brand];
+		B=params->B[brand];
+		mu=params->mu[brand];
 		chi=1.0;
 	}
 	else {
@@ -1768,6 +1804,14 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 		v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
 		v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
 		                                  /* scaled with 0.5 ^ */
+#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
 
 	}
 
@@ -1784,13 +1828,13 @@ int moldyn_bc_check(t_moldyn *moldyn) {
 	t_atom *atom;
 	t_3dvec *dim;
 	int i;
-double x;
-u8 byte;
-int j,k;
+	double x;
+	u8 byte;
+	int j,k;
 
 	atom=moldyn->atom;
 	dim=&(moldyn->dim);
-x=dim->x/2;
+	x=dim->x/2;
 
 	for(i=0;i<moldyn->count;i++) {
 		if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {