X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Falbe_orig.c;h=46e160630cf160dfec4a912f27796d444e84d2b6;hb=8800ba98083a8e5457b035293999f1c70e780c72;hp=c01ab15106d51834980fed6f9c0836c9336efcb1;hpb=2b6a06cdf69b0d8f1dafc0e21b89c6bb4bfc192c;p=physik%2Fposic.git

diff --git a/potentials/albe_orig.c b/potentials/albe_orig.c
index c01ab15..46e1606 100644
--- a/potentials/albe_orig.c
+++ b/potentials/albe_orig.c
@@ -1,5 +1,5 @@
 /*
- * albe.c - albe potential
+ * albe_orig.c - albe potential
  *
  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
  *
@@ -17,23 +17,23 @@
 
 #include "../moldyn.h"
 #include "../math/math.h"
-#include "albe.h"
+#include "albe_orig.h"
 
 /* create mixed terms from parameters and set them */
-int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
+int albe_orig_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 
-	t_albe_mult_params *p;
+	t_albe_orig_mult_params *p;
 
 	// set cutoff before parameters (actually only necessary for some pots)
 	if(moldyn->cutoff==0.0) {
-		printf("[albe] WARNING: no cutoff!\n");
+		printf("[albe orig] WARNING: no cutoff!\n");
 		return -1;
 	}
 
 	/* alloc mem for potential parameters */
 	moldyn->pot_params=malloc(sizeof(t_albe_mult_params));
 	if(moldyn->pot_params==NULL) {
-		perror("[albe] pot params alloc");
+		perror("[albe orig] pot params alloc");
 		return -1;
 	}
 
@@ -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,25 +85,29 @@ 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_m=p->c_mixed*p->c_mixed;
 					p->d_mixed=ALBE_D_SIC;
+					p->d2_m=p->d_mixed*p->d_mixed;
+					p->c2d2_m=p->c2_m/p->d2_m;
 					p->h_mixed=ALBE_H_SIC;
 					break;
 				default:
-					printf("[albe] WARNING: element2\n");
+					printf("[albe orig] WARNING: element2");
+					printf("\n");
 					return -1;
 			}
 			break;
 		default:
-			printf("[albe] WARNING: element1\n");
+			printf("[albe orig] WARNING: element1\n");
 			return -1;
 	}
 
-	printf("[albe] parameter completion\n");
+	printf("[albe orig] parameter completion\n");
 	p->S2[0]=p->S[0]*p->S[0];
 	p->S2[1]=p->S[1]*p->S[1];
 	p->S2mixed=p->Smixed*p->Smixed;
 
-	printf("[albe] mult parameter info:\n");
+	printf("[albe orig] mult parameter info:\n");
 	printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
 	printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
 	printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
@@ -105,19 +115,19 @@ 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("  h      | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
 
 	return 0;
 }
 
 /* albe 3 body potential function (first ij loop) */
-int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+int albe_orig_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
-	t_albe_mult_params *params;
-	t_albe_exchange *exchange;
+	t_albe_orig_mult_params *params;
+	t_albe_orig_exchange *exchange;
 	unsigned char brand;
 	double S2;
 	t_3dvec dist_ij;
@@ -167,11 +177,11 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 }
 
 /* albe 3 body potential function (first k loop) */
-int albe_mult_3bp_k1(t_moldyn *moldyn,
-                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+int albe_orig_mult_3bp_k1(t_moldyn *moldyn,
+                          t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
-	t_albe_mult_params *params;
-	t_albe_exchange *exchange;
+	t_albe_orig_mult_params *params;
+	t_albe_orig_exchange *exchange;
 	unsigned char brand;
 	double R,S,S2;
 	t_3dvec dist_ij,dist_ik;
@@ -180,11 +190,15 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
 	double f_c_ik,df_c_ik;
 	int kcount;
 
+	/* continue if aj equals ak */
+	if(aj==ak)
+		return 0;
+
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 	kcount=exchange->kcount;
 
-	if(kcount>ALBE_MAXN) {
+	if(kcount>ALBE_ORIG_MAXN) {
 		printf("FATAL: neighbours = %d\n",kcount);
 		printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
 	}
@@ -200,6 +214,9 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
 		exchange->c_i=&(params->c[brand]);
 		exchange->d_i=&(params->d[brand]);
 		exchange->h_i=&(params->h[brand]);
+		exchange->ci2=&(params->c2[brand]);
+		exchange->di2=&(params->d2[brand]);
+		exchange->ci2di2=&(params->c2d2[brand]);
 	}
 	else {
 		R=params->Rmixed;
@@ -210,10 +227,10 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
 		exchange->c_i=&(params->c_mixed);
 		exchange->d_i=&(params->d_mixed);
 		exchange->h_i=&(params->h_mixed);
+		exchange->ci2=&(params->c2_m);
+		exchange->di2=&(params->d2_m);
+		exchange->ci2di2=&(params->c2d2_m);
 	}
-	exchange->ci2=*(exchange->c_i)**(exchange->c_i);
-	exchange->di2=*(exchange->d_i)**(exchange->d_i);
-	exchange->ci2di2=exchange->ci2/exchange->di2;
 
 	/* dist_ik, d_ik2 */
 	v3_sub(&dist_ik,&(ak->r),&(ai->r));
@@ -242,9 +259,9 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
 
 	/* g_ijk */
 	h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
-	d2_h_cos2=exchange->di2+(h_cos*h_cos);
-	frac=exchange->ci2/d2_h_cos2;
-	g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
+	d2_h_cos2=*(exchange->di2)+(h_cos*h_cos);
+	frac=*(exchange->ci2)/d2_h_cos2;
+	g=*(exchange->gamma_i)*(1.0+*(exchange->ci2di2)-frac);
 	dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
 
 	/* zeta sum += f_c_ik * g_ijk */
@@ -275,10 +292,10 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
 	return 0;
 }
 
-int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+int albe_orig_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
-	t_albe_mult_params *params;
-	t_albe_exchange *exchange;
+	t_albe_orig_mult_params *params;
+	t_albe_orig_exchange *exchange;
 	t_3dvec force;
 	double f_a,df_a,b,db,f_c,df_c;
 	double f_r,df_r;
@@ -388,11 +405,11 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 }
 
 /* albe 3 body potential function (second k loop) */
-int albe_mult_3bp_k2(t_moldyn *moldyn,
-                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+int albe_orig_mult_3bp_k2(t_moldyn *moldyn,
+                          t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
-	t_albe_mult_params *params;
-	t_albe_exchange *exchange;
+	t_albe_orig_mult_params *params;
+	t_albe_orig_exchange *exchange;
 	int kcount;
 	t_3dvec dist_ik,dist_ij;
 	double d_ik2,d_ik,d_ij2,d_ij;
@@ -405,6 +422,10 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 	t_3dvec dcosdrj,dcosdrk;
 	t_3dvec force,tmp;
 
+	/* continue if aj equals ak */
+	if(aj==ak)
+		return 0;
+
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 	kcount=exchange->kcount;
@@ -521,9 +542,10 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 
 }
 
-int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
+int albe_orig_mult_check_2b_bond(t_moldyn *moldyn,
+                                 t_atom *itom,t_atom *jtom,u8 bc) {
 
-	t_albe_mult_params *params;
+	t_albe_orig_mult_params *params;
 	t_3dvec dist;
 	double d;
 	u8 brand;