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

diff --git a/potentials/albe.c b/potentials/albe.c
index 9d8336e..9f98282 100644
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -105,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);
@@ -301,40 +310,18 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn,
 		dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
 	}
 
-	/* zeta - for albe: ik depending g function */
-//if(ai->tag==0) {
-//	printf("------> %.15f %.15f\n",dj,dk);
-//	printf("------> %.15f %.15f\n",dj,dk);
-//}
-
-	exchange->zeta[j]+=(exchange->f_c[k]*gk);
-	exchange->zeta[k]+=(exchange->f_c[j]*gj);
-
-	/* 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);
+#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
 
-	/* zeta derivatives */
-	dzjj=&(exchange->dzeta[j][j]);
-	dzkk=&(exchange->dzeta[k][k]);
-	dzjk=&(exchange->dzeta[j][k]);
-	dzkj=&(exchange->dzeta[k][j]);
-	v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
-	v3_add(dzjj,dzjj,&tmp);				// j j
-	v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
-	v3_add(dzkk,dzkk,&tmp);				// k k
-	v3_scale(&tmp,&distk,-exchange->df_c[k]*gk);	// dk rik = - di rik
-	v3_add(dzjk,dzjk,&tmp);
-	v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
-	v3_add(dzjk,dzjk,&tmp);				// j k
-	v3_scale(&tmp,&distj,-exchange->df_c[j]*gj);	// dj rij = - di rij
-	v3_add(dzkj,dzkj,&tmp);
-	v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
-	v3_add(dzkj,dzkj,&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;
@@ -440,11 +427,21 @@ printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
 	v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 	v3_add(&(aj->f),&(aj->f),&force);
 
-#ifdef NDEBUG
-if(aj->tag==0) {
-printf("force: %.15f %.15f %.15f | %d %d (ji) %.15f\n",force.x,force.y,force.z,aj->tag,ai->tag,exchange->zeta[j]);
-printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
-}
+	/* virial */
+	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 */
@@ -501,39 +498,89 @@ int albe_mult_i0_j2_k0(t_moldyn *moldyn,
 		return 0;
 	}
 
-	/* k!=j & check whether to run k */
-	k=exchange->kcnt;
-	j=exchange->j2cnt;
-	if((k==j)|(exchange->skip[k])) {
-		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;
 
-#ifdef NDEBUG
-if(ak->tag==0) {
-printf("force: %.15f %.15f %.15f | %d %d %d (k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag);
-printf("    t: %.15f %.15f %.15f\n",ak->f.x,ak->f.y,ak->f.z);
-}
+	/* 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,&(exchange->dist[k]));
+	virial_calc(ai,&force,&dist_ij);
 
+	/* force contribution to atom i */
 	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 %d -- %d(i contr k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag,k);
-printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
-printf("  ## %f\n",exchange->d[k]);
-}
+	/* 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,&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;