X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Falbe.c;h=9f9828280acc5e0c908e44d0aaa0ce4dc5ba6222;hb=6e6d7126ea9a845f11637d8e1b8eb2b570ac4dc9;hp=9dfba13ed9993857b30dfef146e422147c6a44d7;hpb=959801961cfdd51e080e6500f51647b3397d336c;p=physik%2Fposic.git

diff --git a/potentials/albe.c b/potentials/albe.c
index 9dfba13..9f98282 100644
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -53,7 +53,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 			p->mu[0]=ALBE_MU_SI;
 			p->gamma[0]=ALBE_GAMMA_SI;
 			p->c[0]=ALBE_C_SI;
+			p->c2[0]=p->c[0]*p->c[0];
 			p->d[0]=ALBE_D_SI;
+			p->d2[0]=p->d[0]*p->d[0];
+			p->c2d2[0]=p->c2[0]/p->d2[0];
 			p->h[0]=ALBE_H_SI;
 			switch(element2) {
 				case C:
@@ -67,7 +70,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 					p->mu[1]=ALBE_MU_C;
 					p->gamma[1]=ALBE_GAMMA_C;
 					p->c[1]=ALBE_C_C;
+					p->c2[1]=p->c[1]*p->c[1];
 					p->d[1]=ALBE_D_C;
+					p->d2[1]=p->d[1]*p->d[1];
+					p->c2d2[1]=p->c2[1]/p->d2[1];
 					p->h[1]=ALBE_H_C;
 					/* mixed type: silicon carbide */
 					p->Smixed=ALBE_S_SIC;
@@ -79,7 +85,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 					p->mu_m=ALBE_MU_SIC;
 					p->gamma_m=ALBE_GAMMA_SIC;
 					p->c_mixed=ALBE_C_SIC;
+					p->c2_mixed=p->c_mixed*p->c_mixed;
 					p->d_mixed=ALBE_D_SIC;
+					p->d2_mixed=p->d_mixed*p->d_mixed;
+					p->c2d2_m=p->c2_mixed/p->d2_mixed;
 					p->h_mixed=ALBE_H_SIC;
 					break;
 				default:
@@ -96,6 +105,15 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 	p->S2[0]=p->S[0]*p->S[0];
 	p->S2[1]=p->S[1]*p->S[1];
 	p->S2mixed=p->Smixed*p->Smixed;
+	p->c2[0]=p->c[0]*p->c[0];
+	p->c2[1]=p->c[1]*p->c[1];
+	p->c2_mixed=p->c_mixed*p->c_mixed;
+	p->d2[0]=p->d[0]*p->d[0];
+	p->d2[1]=p->d[1]*p->d[1];
+	p->d2_mixed=p->d_mixed*p->d_mixed;
+	p->c2d2[0]=p->c2[0]/p->d2[0];
+	p->c2d2[1]=p->c2[1]/p->d2[1];
+	p->c2d2_m=p->c2_mixed/p->d2_mixed;
 
 	printf("[albe] mult parameter info:\n");
 	printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
@@ -105,10 +123,13 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 	printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
 	                                  p->lambda_m);
 	printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
-	printf("  gamma  | %f | %f\n",p->gamma[0],p->gamma[1]);
-	printf("  c      | %f | %f\n",p->c[0],p->c[1]);
-	printf("  d      | %f | %f\n",p->d[0],p->d[1]);
-	printf("  h      | %f | %f\n",p->h[0],p->h[1]);
+	printf("  gamma  | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m);
+	printf("  c      | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed);
+	printf("  d      | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed);
+	printf("  c2     | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed);
+	printf("  d2     | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed);
+	printf("  c2d2   | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m);
+	printf("  h      | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
 
 	return 0;
 }
@@ -116,21 +137,18 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 /* first i loop */
 int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
 
-	int i;
 	t_albe_mult_params *params;
 	t_albe_exchange *exchange;
 	
+	int i;
+
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 
 	/* zero exchange values */
-	memset(exchange->dist,0,ALBE_MAXN*sizeof(t_3dvec));
-	memset(exchange->d2,0,ALBE_MAXN*sizeof(double));
-	memset(exchange->d,0,ALBE_MAXN*sizeof(double));
 	memset(exchange->zeta,0,ALBE_MAXN*sizeof(double));
 	for(i=0;i<ALBE_MAXN;i++)
-		memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(double));
-	memset(exchange->skip,0,ALBE_MAXN*sizeof(u8));
+		memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(t_3dvec));
 	exchange->jcnt=0;
 	exchange->j2cnt=0;
 
@@ -142,6 +160,7 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
 	t_albe_mult_params *params;
 	t_albe_exchange *exchange;
+
 	double S2,S,R,d2,d,s_r,arg;
 	t_3dvec dist;
 	int j;
@@ -164,8 +183,8 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
 	/* dist_ij, d_ij2 */
 	v3_sub(&dist,&(aj->r),&(ai->r));
-	exchange->dist[j]=dist;
 	if(bc) check_per_bound(moldyn,&dist);
+	exchange->dist[j]=dist;
 	d2=v3_absolute_square(&dist);
 	exchange->d2[j]=d2;
 
@@ -176,16 +195,20 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		exchange->jcnt+=1;
 		return 0;
 	}
+	exchange->skip[j]=0;
 
 	/* more ij depending values */
 	if(brand==aj->brand) {
 		R=params->R[brand];
 		S=params->S[brand];
-		/* set albe needs i,(j/k) depending c,d,h and gamma values */
+		/* albe needs i,(j/k) depending c,d,h and gamma values */
 		exchange->gamma_[j]=&(params->gamma[brand]);
 		exchange->c_[j]=&(params->c[brand]);
 		exchange->d_[j]=&(params->d[brand]);
 		exchange->h_[j]=&(params->h[brand]);
+		exchange->c2_[j]=&(params->c2[brand]);
+		exchange->d2_[j]=&(params->d2[brand]);
+		exchange->c2d2_[j]=&(params->c2d2[brand]);
 	}
 	else {
 		R=params->Rmixed;
@@ -195,13 +218,13 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		exchange->c_[j]=&(params->c_mixed);
 		exchange->d_[j]=&(params->d_mixed);
 		exchange->h_[j]=&(params->h_mixed);
+		exchange->c2_[j]=&(params->c2_mixed);
+		exchange->d2_[j]=&(params->d2_mixed);
+		exchange->c2d2_[j]=&(params->c2d2_m);
 	}
-	exchange->c2_[j]=*exchange->c_[j]**exchange->c_[j];
-	exchange->d2_[j]=*exchange->d_[j]**exchange->d_[j];
-	exchange->c2d2_[j]=exchange->c2_[j]/exchange->d2_[j];
 
 	/* d_ij */
-	d=sqrt(exchange->d2[j]);
+	d=sqrt(d2);
 	exchange->d[j]=d;
 	
 	/* f_c, df_c */
@@ -228,16 +251,23 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
 
 	t_albe_mult_params *params;
 	t_albe_exchange *exchange;
+
 	int j,k;
 	t_3dvec distj,distk;
 	double dj,dk,djdk_inv,cos_theta;
 	double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j;
 	double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k;
-	t_3dvec dcosdrj,dcosdrk,tmp,**dzeta;
+	t_3dvec dcosdrj,dcosdrk,tmp;
+	t_3dvec *dzjj,*dzkk,*dzjk,*dzkj;
 
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 
+	if(aj==ak) {
+		exchange->kcnt+=1;
+		return 0;
+	}
+
 	/* k<j & check whether to run k */
 	j=exchange->jcnt;
 	k=exchange->kcnt;
@@ -249,7 +279,7 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
 		exchange->kcnt+=1;
 		return 0;
 	}
-	
+
 	/* distances */
 	distj=exchange->dist[j];
 	distk=exchange->dist[k];
@@ -260,11 +290,11 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
 	/* cos theta */
 	cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
 
-	/* g(cos(theta)) - in albe: ik-depending values! */
+	/* g(cos(theta)) ij and ik values */
 	h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism
-	d2_h_cos2_j=exchange->d2_[j]+(h_cos_j*h_cos_j);
-	frac_j=exchange->c2_[j]/d2_h_cos2_j;
-	gj=1.0+exchange->c2d2_[j]-frac_j;
+	d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j);
+	frac_j=*exchange->c2_[j]/d2_h_cos2_j;
+	gj=1.0+*exchange->c2d2_[j]-frac_j;
 	gj*=*(exchange->gamma_[j]);
 	dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe 
 	if(ak->brand==aj->brand) {
@@ -273,39 +303,25 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
 	}
 	else {
 		h_cos_k=*(exchange->h_[k])+cos_theta;
-		d2_h_cos2_k=exchange->d2_[k]+(h_cos_k*h_cos_k);
-		frac_k=exchange->c2_[k]/d2_h_cos2_k;
-		gk=1.0+exchange->c2d2_[k]-frac_k;
+		d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k);
+		frac_k=*exchange->c2_[k]/d2_h_cos2_k;
+		gk=1.0+*exchange->c2d2_[k]-frac_k;
 		gk*=*(exchange->gamma_[k]);
 		dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
 	}
 
-	/* zeta */
-	exchange->zeta[j]+=(exchange->f_c[k]*gk);
-	exchange->zeta[k]+=(exchange->f_c[j]*gj);
+#ifdef DEBUG
+	if(ai==&(moldyn->atom[DATOM])) 
+		printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik);
+#endif
 
-	/* cos theta derivatives */
-	v3_scale(&dcosdrj,&distk,djdk_inv);		// j
-	v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]);
-	v3_add(&dcosdrj,&dcosdrj,&tmp);
-	v3_scale(&dcosdrk,&distj,djdk_inv);		// k
-	v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]);
-	v3_add(&dcosdrk,&dcosdrk,&tmp);
-
-	/* zeta derivatives */
-	dzeta=exchange->dzeta;
-	v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
-	v3_add(&dzeta[j][j],&dzeta[j][j],&tmp);		// j j
-	v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
-	v3_add(&dzeta[k][k],&dzeta[k][k],&tmp);		// k k
-	v3_scale(&tmp,&distk,exchange->df_c[k]*gk/dk);
-	v3_add(&dzeta[j][k],&dzeta[j][k],&tmp);
-	v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
-	v3_add(&dzeta[j][k],&dzeta[j][k],&tmp);		// j k
-	v3_scale(&tmp,&distj,exchange->df_c[j]*gj/dj);
-	v3_add(&dzeta[k][j],&dzeta[k][j],&tmp);
-	v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
-	v3_add(&dzeta[k][j],&dzeta[k][j],&tmp);		// k j
+	/* store even more data for second k loop */
+	exchange->g[kcount]=g;
+	exchange->dg[kcount]=dg;
+	exchange->d_ik[kcount]=d_ik;
+	exchange->cos_theta[kcount]=cos_theta;
+	exchange->f_c_ik[kcount]=f_c_ik;
+	exchange->df_c_ik[kcount]=df_c_ik;
 
 	/* increase k counter */
 	exchange->kcnt+=1;
@@ -400,25 +416,64 @@ int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	v3_scale(&force,dist,scale);
 	v3_add(&(ai->f),&(ai->f),&force);
 
+#ifdef NDEBUG
+if(ai->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]);
+printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
+}
+#endif
+
 	/* force contribution for atom j due to ij bond */
 	v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 	v3_add(&(aj->f),&(aj->f),&force);
 
 	/* virial */
-	virial_calc(aj,&force,dist);
+	virial_calc(ai,&force,&(exchange->dist_ij));
+
+#ifdef DEBUG
+	if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
+		printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
+		printf("  adding %f %f %f\n",force.x,force.y,force.z);
+		if(ai==&(moldyn->atom[DATOM]))
+			printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
+		if(aj==&(moldyn->atom[DATOM]))
+			printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+		printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
+		                                    f_c,b,f_a,f_r);
+		printf("          %f %f %f\n",exchange->zeta_ij,.0,.0);
+	}
+#endif
+
+	/* virial */
+	virial_calc(ai,&force,dist);
 
 	/* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
-	exchange->pre_dzeta=0.5*f_a*exchange->f_c[j]*db;
+	exchange->pre_dzeta=0.5*f_a*f_c*db;
 
 	/* force contribution (drj derivative) */
 	v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta);
 	v3_add(&(aj->f),&(aj->f),&force);
-	v3_scale(&force,&force,-1.0);
-	v3_add(&(ai->f),&(ai->f),&force);
+
+#ifdef NDEBUG
+if(aj->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag);
+printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
+}
+#endif
 
 	/* virial */
 	virial_calc(ai,&force,dist);
 
+	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
+
+#ifdef NDEBUG
+if(ai->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (i contr j der)\n",force.x,force.y,force.z,ai->tag,aj->tag);
+printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
+}
+#endif
+
 	/* reset k counter for second k loop */
 	exchange->kcnt=0;
 		
@@ -438,22 +493,93 @@ int albe_mult_i0_j2_k0(t_moldyn *moldyn,
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 
-	/* k!=j & check whether to run k */
-	k=exchange->kcnt;
-	j=exchange->j2cnt;
-	if((k==j)|(exchange->skip[k])) {
+	if(aj==ak) {
 		exchange->kcnt+=1;
 		return 0;
 	}
-	
-	/* force contribution (drk derivative) */
-	v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta);
-	v3_add(&(ak->f),&(ak->f),&force);
+
+	/* prefactor dzeta */
+	pre_dzeta=exchange->pre_dzeta;
+
+	/* dist_ik, d_ik */
+	dist_ik=exchange->dist_ik[kcount];
+	d_ik=exchange->d_ik[kcount];
+
+	/* f_c_ik, df_c_ik */
+	f_c_ik=exchange->f_c_ik[kcount];
+	df_c_ik=exchange->df_c_ik[kcount];
+
+	/* dist_ij, d_ij2, d_ij */
+	dist_ij=exchange->dist_ij;
+	d_ij2=exchange->d_ij2;
+	d_ij=exchange->d_ij;
+
+	/* g, dg, cos_theta */
+	g=exchange->g[kcount];
+	dg=exchange->dg[kcount];
+	cos_theta=exchange->cos_theta[kcount];
+
+	/* cos_theta derivatives wrt j,k */
+	dijdik_inv=1.0/(d_ij*d_ik);
+	v3_scale(&dcosdrj,&dist_ik,dijdik_inv);		// j
+	v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
+	v3_add(&dcosdrj,&dcosdrj,&tmp);
+	v3_scale(&dcosdrk,&dist_ij,dijdik_inv);		// k
+	v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
+	v3_add(&dcosdrk,&dcosdrk,&tmp);
+
+	/* f_c_ik * dg, df_c_ik * g */
+	fcdg=f_c_ik*dg;
+	dfcg=df_c_ik*g;
+
+	/* derivative wrt j */
+	v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
+
+	/* force contribution */
+	v3_add(&(aj->f),&(aj->f),&force);
+
+#ifdef DEBUG
+	if(aj==&(moldyn->atom[DATOM])) {
+		printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
+		printf("  adding %f %f %f\n",force.x,force.y,force.z);
+		printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+		printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
+		printf("    d ij ik = %f %f\n",d_ij,d_ik);
+	}
+#endif
+
+	/* virial */
+	virial_calc(ai,&force,&dist_ij);
+
+	/* force contribution to atom i */
 	v3_scale(&force,&force,-1.0);
 	v3_add(&(ai->f),&(ai->f),&force);
 
+	/* derivative wrt k */
+	v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
+	v3_scale(&tmp,&dcosdrk,fcdg);
+	v3_add(&force,&force,&tmp);
+	v3_scale(&force,&force,pre_dzeta);
+
+	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
+
+#ifdef DEBUG
+	if(ak==&(moldyn->atom[DATOM])) {
+		printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
+		printf("  adding %f %f %f\n",force.x,force.y,force.z);
+		printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
+		printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
+		printf("    d ij ik = %f %f\n",d_ij,d_ik);
+	}
+#endif
+
 	/* virial */
-	virial_calc(ai,&force,&(exchange->dist[k]));
+	virial_calc(ai,&force,&dist_ik);
+	
+	/* force contribution to atom i */
+	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
 
 	/* increase k counter */
 	exchange->kcnt+=1;