From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Thu, 10 Jul 2008 20:11:59 +0000 (+0200)
Subject: secured original albe implementation
X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=7c856cfd07e2e9d57ff0af3172b7f271fc5b22e0;p=physik%2Fposic.git

secured original albe implementation
---

diff --git a/potentials/albe_orig.c b/potentials/albe_orig.c
new file mode 100644
index 0000000..c01ab15
--- /dev/null
+++ b/potentials/albe_orig.c
@@ -0,0 +1,548 @@
+/*
+ * albe.c - albe potential
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <math.h>
+
+#include "../moldyn.h"
+#include "../math/math.h"
+#include "albe.h"
+
+/* create mixed terms from parameters and set them */
+int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
+
+	t_albe_mult_params *p;
+
+	// set cutoff before parameters (actually only necessary for some pots)
+	if(moldyn->cutoff==0.0) {
+		printf("[albe] 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");
+		return -1;
+	}
+
+	/* these are now albe parameters */
+	p=moldyn->pot_params;
+
+	// only 1 combination by now :p
+	switch(element1) {
+		case SI:
+			/* type: silicon */
+			p->S[0]=ALBE_S_SI;
+			p->R[0]=ALBE_R_SI;
+			p->A[0]=ALBE_A_SI;
+			p->B[0]=ALBE_B_SI;
+			p->r0[0]=ALBE_R0_SI;
+			p->lambda[0]=ALBE_LAMBDA_SI;
+			p->mu[0]=ALBE_MU_SI;
+			p->gamma[0]=ALBE_GAMMA_SI;
+			p->c[0]=ALBE_C_SI;
+			p->d[0]=ALBE_D_SI;
+			p->h[0]=ALBE_H_SI;
+			switch(element2) {
+				case C:
+					/* type: carbon */
+					p->S[1]=ALBE_S_C;
+					p->R[1]=ALBE_R_C;
+					p->A[1]=ALBE_A_C;
+					p->B[1]=ALBE_B_C;
+					p->r0[1]=ALBE_R0_C;
+					p->lambda[1]=ALBE_LAMBDA_C;
+					p->mu[1]=ALBE_MU_C;
+					p->gamma[1]=ALBE_GAMMA_C;
+					p->c[1]=ALBE_C_C;
+					p->d[1]=ALBE_D_C;
+					p->h[1]=ALBE_H_C;
+					/* mixed type: silicon carbide */
+					p->Smixed=ALBE_S_SIC;
+					p->Rmixed=ALBE_R_SIC;
+					p->Amixed=ALBE_A_SIC;
+					p->Bmixed=ALBE_B_SIC;
+					p->r0_mixed=ALBE_R0_SIC;
+					p->lambda_m=ALBE_LAMBDA_SIC;
+					p->mu_m=ALBE_MU_SIC;
+					p->gamma_m=ALBE_GAMMA_SIC;
+					p->c_mixed=ALBE_C_SIC;
+					p->d_mixed=ALBE_D_SIC;
+					p->h_mixed=ALBE_H_SIC;
+					break;
+				default:
+					printf("[albe] WARNING: element2\n");
+					return -1;
+			}
+			break;
+		default:
+			printf("[albe] WARNING: element1\n");
+			return -1;
+	}
+
+	printf("[albe] 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("  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);
+	printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
+	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]);
+
+	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) {
+
+	t_albe_mult_params *params;
+	t_albe_exchange *exchange;
+	unsigned char brand;
+	double S2;
+	t_3dvec dist_ij;
+	double d_ij2,d_ij;
+
+	params=moldyn->pot_params;
+	exchange=&(params->exchange);
+
+	/* reset zeta sum */
+	exchange->zeta_ij=0.0;
+
+	/*
+	 * set ij depending values
+	 */
+
+	brand=ai->brand;
+	if(brand==aj->brand) {
+		S2=params->S2[brand];
+	}
+	else {
+		S2=params->S2mixed;
+	}
+
+	/* dist_ij, d_ij2 */
+	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) {
+		moldyn->run3bp=0;
+		return 0;
+	}
+
+	/* d_ij */
+	d_ij=sqrt(d_ij2);
+
+	/* store values */
+	exchange->dist_ij=dist_ij;
+	exchange->d_ij2=d_ij2;
+	exchange->d_ij=d_ij;
+
+	/* reset k counter for first k loop */
+	exchange->kcount=0;
+		
+	return 0;
+}
+
+/* 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) {
+
+	t_albe_mult_params *params;
+	t_albe_exchange *exchange;
+	unsigned char brand;
+	double R,S,S2;
+	t_3dvec dist_ij,dist_ik;
+	double d_ik2,d_ik,d_ij;
+	double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
+	double f_c_ik,df_c_ik;
+	int kcount;
+
+	params=moldyn->pot_params;
+	exchange=&(params->exchange);
+	kcount=exchange->kcount;
+
+	if(kcount>ALBE_MAXN) {
+		printf("FATAL: neighbours = %d\n",kcount);
+		printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
+	}
+
+	/* ik constants */
+	brand=ai->brand;
+	if(brand==ak->brand) {
+		R=params->R[brand];
+		S=params->S[brand];
+		S2=params->S2[brand];
+		/* albe needs i,k depending c,d,h and gamma values */
+		exchange->gamma_i=&(params->gamma[brand]);
+		exchange->c_i=&(params->c[brand]);
+		exchange->d_i=&(params->d[brand]);
+		exchange->h_i=&(params->h[brand]);
+	}
+	else {
+		R=params->Rmixed;
+		S=params->Smixed;
+		S2=params->S2mixed;
+		/* albe needs i,k depending c,d,h and gamma values */
+		exchange->gamma_i=&(params->gamma_m);
+		exchange->c_i=&(params->c_mixed);
+		exchange->d_i=&(params->d_mixed);
+		exchange->h_i=&(params->h_mixed);
+	}
+	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));
+	if(bc) check_per_bound(moldyn,&dist_ik);
+	d_ik2=v3_absolute_square(&dist_ik);
+
+	/* store data for second k loop */
+	exchange->dist_ik[kcount]=dist_ik;
+	exchange->d_ik2[kcount]=d_ik2;
+
+	/* return if not within cutoff */
+	if(d_ik2>S2) {
+		exchange->kcount++;
+		return 0;
+	}
+
+	/* d_ik */
+	d_ik=sqrt(d_ik2);
+
+	/* dist_ij, d_ij */
+	dist_ij=exchange->dist_ij;
+	d_ij=exchange->d_ij;
+
+	/* cos theta */
+	cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
+
+	/* 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);
+	dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
+
+	/* zeta sum += f_c_ik * g_ijk */
+	if(d_ik<=R) {
+		exchange->zeta_ij+=g;
+		f_c_ik=1.0;
+		df_c_ik=0.0;
+	}
+	else {
+		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));
+		exchange->zeta_ij+=f_c_ik*g;
+	}
+
+	/* 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->kcount++;
+
+	return 0;
+}
+
+int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+	t_albe_mult_params *params;
+	t_albe_exchange *exchange;
+	t_3dvec force;
+	double f_a,df_a,b,db,f_c,df_c;
+	double f_r,df_r;
+	double scale;
+	double mu,B;
+	double lambda,A;
+	double d_ij,r0;
+	unsigned char brand;
+	double S,R,s_r,arg;
+	double energy;
+
+	params=moldyn->pot_params;
+	exchange=&(params->exchange);
+
+	brand=aj->brand;
+	if(brand==ai->brand) {
+		S=params->S[brand];
+		R=params->R[brand];
+		B=params->B[brand];
+		A=params->A[brand];
+		r0=params->r0[brand];
+		mu=params->mu[brand];
+		lambda=params->lambda[brand];
+	}
+	else {
+		S=params->Smixed;
+		R=params->Rmixed;
+		B=params->Bmixed;
+		A=params->Amixed;
+		r0=params->r0_mixed;
+		mu=params->mu_m;
+		lambda=params->lambda_m;
+	}
+
+	d_ij=exchange->d_ij;
+
+	/* f_c, df_c */
+	if(d_ij<R) {
+		f_c=1.0;
+		df_c=0.0;
+	}
+	else {
+		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));
+	}
+
+	/* f_a, df_a */
+	f_a=-B*exp(-mu*(d_ij-r0));
+	df_a=mu*f_a/d_ij;
+
+	/* f_r, df_r */
+	f_r=A*exp(-lambda*(d_ij-r0));
+	df_r=lambda*f_r/d_ij;
+
+	/* b, db */
+	if(exchange->zeta_ij==0.0) {
+		b=1.0;
+		db=0.0;
+	}
+	else {
+		b=1.0/sqrt(1.0+exchange->zeta_ij);
+		db=-0.5*b/(1.0+exchange->zeta_ij);
+	}
+
+	/* force contribution for atom i */
+	scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
+	v3_scale(&force,&(exchange->dist_ij),scale);
+	v3_add(&(ai->f),&(ai->f),&force);
+
+	/* force contribution for atom j */
+	v3_scale(&force,&force,-1.0); // dri rij = - drj rij
+	v3_add(&(aj->f),&(aj->f),&force);
+
+	/* virial */
+	virial_calc(aj,&force,&(exchange->dist_ij));
+
+#ifdef DEBUG
+if(moldyn->time>DSTART&&moldyn->time<DEND) {
+	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[0]))
+			printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
+		if(aj==&(moldyn->atom[0]))
+			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
+
+	/* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
+	exchange->pre_dzeta=0.5*f_a*f_c*db;
+
+	/* energy contribution */
+	energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
+	moldyn->energy+=energy;
+	ai->e+=energy;
+
+	/* reset k counter for second k loop */
+	exchange->kcount=0;
+		
+	return 0;
+}
+
+/* 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) {
+
+	t_albe_mult_params *params;
+	t_albe_exchange *exchange;
+	int kcount;
+	t_3dvec dist_ik,dist_ij;
+	double d_ik2,d_ik,d_ij2,d_ij;
+	unsigned char brand;
+	double S2;
+	double g,dg,cos_theta;
+	double pre_dzeta;
+	double f_c_ik,df_c_ik;
+	double dijdik_inv,fcdg,dfcg;
+	t_3dvec dcosdrj,dcosdrk;
+	t_3dvec force,tmp;
+
+	params=moldyn->pot_params;
+	exchange=&(params->exchange);
+	kcount=exchange->kcount;
+
+	if(kcount>ALBE_MAXN)
+		printf("FATAL: neighbours!\n");
+
+	/* d_ik2 */
+	d_ik2=exchange->d_ik2[kcount];
+
+	brand=ak->brand;
+	if(brand==ai->brand)
+		S2=params->S2[brand];
+	else
+		S2=params->S2mixed;
+
+	/* return if d_ik > S */
+	if(d_ik2>S2) {
+		exchange->kcount++;
+		return 0;
+	}
+
+	/* 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(moldyn->time>DSTART&&moldyn->time<DEND) {
+	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
+
+	/* force contribution to atom i */
+	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
+
+	/* virial */
+	virial_calc(ai,&force,&dist_ij);
+
+	/* 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);
+
+	/* force contribution */
+	v3_add(&(ak->f),&(ak->f),&force);
+
+#ifdef DEBUG
+if(moldyn->time>DSTART&&moldyn->time<DEND) {
+	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
+
+	/* force contribution to atom i */
+	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
+
+	/* virial */
+	virial_calc(ai,&force,&dist_ik);
+	
+	/* increase k counter */
+	exchange->kcount++;	
+
+	return 0;
+
+}
+
+int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
+
+	t_albe_mult_params *params;
+	t_3dvec dist;
+	double d;
+	u8 brand;
+
+	v3_sub(&dist,&(jtom->r),&(itom->r));
+	if(bc) check_per_bound(moldyn,&dist);
+	d=v3_absolute_square(&dist);
+
+	params=moldyn->pot_params;
+	brand=itom->brand;
+
+	if(brand==jtom->brand) {
+		if(d<=params->S2[brand])
+			return TRUE;
+	}
+	else {
+		if(d<=params->S2mixed)
+			return TRUE;
+	}
+
+	return FALSE;
+}
diff --git a/potentials/albe_orig.h b/potentials/albe_orig.h
new file mode 100644
index 0000000..15386e7
--- /dev/null
+++ b/potentials/albe_orig.h
@@ -0,0 +1,131 @@
+/*
+ * albe.h - albe potential header file
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#ifndef ALBE_H
+#define ALBE_H
+
+#define ALBE_MAXN	16*27
+
+/* albe exchange type */
+typedef struct s_albe_exchange {
+
+	t_3dvec dist_ij;
+	double d_ij2;
+	double d_ij;
+
+	t_3dvec dist_ik[ALBE_MAXN];
+	double d_ik2[ALBE_MAXN];
+	double d_ik[ALBE_MAXN];
+
+	double f_c_ik[ALBE_MAXN];
+	double df_c_ik[ALBE_MAXN];
+
+	double g[ALBE_MAXN];
+	double dg[ALBE_MAXN];
+	double cos_theta[ALBE_MAXN];
+
+	double *gamma_i;
+	double *c_i;
+	double *d_i;
+	double *h_i;
+
+	double ci2;
+	double di2;
+	double ci2di2;
+
+	double zeta_ij;
+	double pre_dzeta;
+
+	int kcount;
+} t_albe_exchange;
+
+/* albe mult (2!) potential parameters */
+typedef struct s_albe_mult_params {
+	double S[2];		/* albe cutoff radii */
+	double S2[2];		/* albe cutoff radii squared */
+	double R[2];		/* albe cutoff radii */
+	double Smixed;		/* mixed S radius */
+	double S2mixed;		/* mixed S radius squared */
+	double Rmixed;		/* mixed R radius */
+	double A[2];		/* factor of albe attractive part */
+	double B[2];		/* factor of albe repulsive part */
+	double r0[2];		/* r_0 */
+	double Amixed;		/* mixed A factor */
+	double Bmixed;		/* mixed B factor */
+	double r0_mixed;	/* mixed r_0 */
+	double lambda[2];	/* albe lambda */
+	double lambda_m;	/* mixed lambda */
+	double mu[2];		/* albe mu */
+	double mu_m;		/* mixed mu */
+
+	double gamma[2];
+	double gamma_m;
+	double c[2];
+	double c_mixed;
+	double d[2];
+	double d_mixed;
+	double h[2];
+	double h_mixed;
+
+	t_albe_exchange exchange;	/* exchange between 2bp and 3bp calc */
+} t_albe_mult_params;
+
+/* function prototypes */
+int albe_mult_set_params(t_moldyn *moldyn,int element1,int elemnt2);
+int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int albe_mult_3bp_k1(t_moldyn *moldyn,
+                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
+int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int albe_mult_3bp_k2(t_moldyn *moldyn,
+                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
+int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc);
+
+/* albe potential parameter defines */
+
+// silicon
+#define ALBE_R_SI		(2.82-0.14)
+#define ALBE_S_SI		(2.82+0.14)
+#define ALBE_A_SI		(3.24*EV/0.842)
+#define ALBE_B_SI		(-1.842*3.24*EV/0.842)
+#define ALBE_R0_SI		2.232
+#define ALBE_LAMBDA_SI		(1.4761*sqrt(2.0*1.842))
+#define ALBE_MU_SI		(1.4761*sqrt(2.0/1.842))
+#define ALBE_GAMMA_SI		0.114354
+#define ALBE_C_SI		2.00494
+#define ALBE_D_SI		0.81472
+#define ALBE_H_SI		0.259
+#define ALBE_LC_SI		5.429
+
+// carbon
+#define ALBE_R_C		(2.00-0.15)
+#define ALBE_S_C		(2.00+0.15)
+#define ALBE_A_C		(6.00*EV/1.167)
+#define ALBE_B_C		(-2.167*6.00*EV/1.167)
+#define ALBE_R0_C		1.4276
+#define ALBE_LAMBDA_C		(2.0099*sqrt(2.0*2.167))
+#define ALBE_MU_C		(2.0099*sqrt(2.0/2.167))
+#define ALBE_GAMMA_C		0.11233
+#define ALBE_C_C		181.910
+#define ALBE_D_C		6.28433
+#define ALBE_H_C		0.5556
+#define ALBE_LC_C		3.566
+
+// mixed: silicon carbide
+#define ALBE_R_SIC		(2.40-0.20)
+#define ALBE_S_SIC		(2.40+0.20)
+#define ALBE_A_SIC		(4.36*EV/0.847)
+#define ALBE_B_SIC		(-1.847*4.36*EV/0.847)
+#define ALBE_R0_SIC		1.79
+#define ALBE_LAMBDA_SIC		(1.6991*sqrt(2.0*1.847))
+#define ALBE_MU_SIC		(1.6991*sqrt(2.0/1.847))
+#define ALBE_GAMMA_SIC		0.011877
+#define ALBE_C_SIC		273987
+#define ALBE_D_SIC		180.314
+#define ALBE_H_SIC		0.68
+#define ALBE_LC_SIC		4.359
+
+#endif