]> hackdaworld.org Git - physik/posic.git/commitdiff
added albe potential (still wrong energy!)
authorhackbard <hackbard>
Tue, 22 May 2007 15:39:26 +0000 (15:39 +0000)
committerhackbard <hackbard>
Tue, 22 May 2007 15:39:26 +0000 (15:39 +0000)
Makefile
moldyn.h
potentials/albe.c [new file with mode: 0644]
potentials/albe.h [new file with mode: 0644]
sic.c

index ebad7e3edb0ea3083fc8e0ac6755a26fda7193e9..bb520c4a81c214d234dbf00d18a3267ca24a0aa5 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -7,12 +7,14 @@ CFLAGS+=-g
 #CFLAGS+=-DDEBUG
 #CFLAGS+=-DVDEBUG
 #CFLAGS+=-DTERSOFF_ORIG
+CFLAGS+=-DALBE
 LDFLAGS=-lm
 
 SOURCE=moldyn.c visual/visual.c random/random.c list/list.c
 
 #POT_SRC=potentials/tersoff_orig.c
 POT_SRC=potentials/tersoff.c
+POT_SRC+=potentials/albe.c
 POT_SRC+= potentials/lennard_jones.c potentials/harmonic_oscillator.c
 
 all: sic
index 3bc40c0cd86235e282dd0daf79b46f9488ad3519..86444f4a6ea3ef5e68ce0e4de4bfff171d5c4b2d 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -242,38 +242,39 @@ typedef struct s_moldyn {
 #define ONE_THIRD              (1.0/3.0)
 
 #define C                      0x06
-#define LC_C                   (0.3567e-9*METER)               /* A */
+//#define LC_C                 3.567                           /* A */
+#define LC_C                   3.560                           /* A */
 #define M_C                    12.011                          /* amu */
 
 #define SI                     0x0e
-#define LC_SI                  (0.543105e-9*METER)             /* A */
+#define LC_SI                  5.43105                         /* A */
 #define M_SI                   28.08553                        /* amu */
 
-#define LC_3C_SIC              (0.43596e-9*METER)              /* A */
+#define LC_3C_SIC              4.3596                          /* A */
 
 #define LJ_SIGMA_SI            ((0.25*sqrt(3.0)*LC_SI)/1.122462)       /* A */
 //#define LJ_SIGMA_SI          (LC_SI/1.122462)                        /* A */
 //#define LJ_SIGMA_SI          (0.5*sqrt(2.0)*LC_SI/1.122462)                  /* A */
 #define LJ_EPSILON_SI          (2.1678*EV)                             /* NA */
 
-#define TM_R_SI                        (2.7e-10*METER)                 /* A */
-#define TM_S_SI                        (3.0e-10*METER)                 /* A */
+#define TM_R_SI                        2.7                             /* A */
+#define TM_S_SI                        3.0                             /* A */
 #define TM_A_SI                        (1830.8*EV)                     /* NA */
 #define TM_B_SI                        (471.18*EV)                     /* NA */
-#define TM_LAMBDA_SI           (2.4799e10/METER)               /* 1/A */
-#define TM_MU_SI               (1.7322e10/METER)               /* 1/A */
+#define TM_LAMBDA_SI           2.4799                          /* 1/A */
+#define TM_MU_SI               1.7322                          /* 1/A */
 #define TM_BETA_SI             1.1000e-6
 #define TM_N_SI                        0.78734
 #define TM_C_SI                        1.0039e5
 #define TM_D_SI                        16.217
 #define TM_H_SI                        -0.59825
 
-#define TM_R_C                 (1.8e-10*METER)                 /* A */
-#define TM_S_C                 (2.1e-10*METER)                 /* A */
+#define TM_R_C                 1.8                             /* A */
+#define TM_S_C                 2.1                             /* A */
 #define TM_A_C                 (1393.6*EV)                     /* NA */
 #define TM_B_C                 (346.7*EV)                      /* NA */
-#define TM_LAMBDA_C            (3.4879e10/METER)               /* 1/A */
-#define TM_MU_C                        (2.2119e10/METER)               /* 1/A */
+#define TM_LAMBDA_C            3.4879                          /* 1/A */
+#define TM_MU_C                        2.2119                          /* 1/A */
 #define TM_BETA_C              1.5724e-7
 #define TM_N_C                 0.72751
 #define TM_C_C                 3.8049e4
@@ -284,6 +285,44 @@ typedef struct s_moldyn {
 
 #define TM_LC_3C_SIC           (0.432e-9*METER)                /* A */
 
+#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_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_R_SIC             (2.40-0.20)
+#define ALBE_S_SIC             (2.40+0.10)
+#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_CHI_SIC           1.0
+
 /*
  * lattice constants
  */
diff --git a/potentials/albe.c b/potentials/albe.c
new file mode 100644 (file)
index 0000000..c2acf2f
--- /dev/null
@@ -0,0 +1,468 @@
+/*
+ * 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_complete_params(t_albe_mult_params *p) {
+
+       printf("[moldyn] 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("[moldyn] 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 1 body part */
+int albe_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
+
+       int brand;
+       t_albe_mult_params *params;
+       t_albe_exchange *exchange;
+       
+       brand=ai->brand;
+       params=moldyn->pot_params;
+       exchange=&(params->exchange);
+
+       /*
+        * simple: point constant parameters only depending on atom i to
+        *         their right 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]);
+
+       exchange->ci2=params->c[brand]*params->c[brand];
+       exchange->di2=params->d[brand]*params->d[brand];
+       exchange->ci2di2=exchange->ci2/exchange->di2;
+
+       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];
+       }
+       else {
+               R=params->Rmixed;
+               S=params->Smixed;
+               S2=params->S2mixed;
+       }
+
+       /* 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;
+       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;
+
+       /* 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;
+
+       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 */
+       scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_a));
+       v3_scale(&force,&(exchange->dist_ij),scale);
+       v3_add(&(ai->f),&(ai->f),&force);
+       v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
+
+#ifdef DEBUG
+       if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
+               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);
+       }
+#endif
+
+       /* virial */
+       if(aj<ai)
+               virial_calc(ai,&force,&(exchange->dist_ij));
+
+       /* dzeta prefactor = - 0.5 f_c f_a db */
+       exchange->pre_dzeta=-0.5*f_a*f_c*db;
+
+       /* energy contribution */
+       moldyn->energy+=0.5*f_c*(f_r+b*f_a);
+
+       /* 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 dcosdri,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 i,j,k */
+       dijdik_inv=1.0/(d_ij*d_ik);
+       v3_scale(&dcosdrj,&dist_ik,dijdik_inv);
+       v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
+       v3_add(&dcosdrj,&dcosdrj,&tmp);
+       v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
+       v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
+       v3_add(&dcosdrk,&dcosdrk,&tmp);
+       v3_add(&dcosdri,&dcosdrj,&dcosdrk);
+       v3_scale(&dcosdri,&dcosdri,-1.0);
+
+       /* f_c_ik * dg, df_c_ik * g */
+       fcdg=f_c_ik*dg;
+       dfcg=df_c_ik*g;
+
+       /* derivative wrt i */
+       v3_scale(&force,&dist_ik,dfcg);
+       v3_scale(&tmp,&dcosdri,fcdg);
+       v3_add(&force,&force,&tmp);
+       v3_scale(&force,&force,pre_dzeta);
+
+       /* force contribution */
+       v3_add(&(ai->f),&(ai->f),&force);
+       
+#ifdef DEBUG
+       if(ai==&(moldyn->atom[0])) {
+               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 i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
+       }
+#endif
+
+       /* virial */
+       //virial_calc(ai,&force,&dist_ij);
+
+       /* 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[0])) {
+               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);
+       }
+#endif
+
+       /* virial */
+       //v3_scale(&force,&force,-1.0);
+       if(aj<ai)
+               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(ak==&(moldyn->atom[0])) {
+               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);
+       }
+#endif
+
+       /* virial */
+       //v3_scale(&force,&force,-1.0);
+       if(aj<ai)
+               virial_calc(ai,&force,&dist_ik);
+       
+       /* increase k counter */
+       exchange->kcount++;     
+
+       return 0;
+
+}
diff --git a/potentials/albe.h b/potentials/albe.h
new file mode 100644 (file)
index 0000000..35a0762
--- /dev/null
@@ -0,0 +1,84 @@
+/*
+ * 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 c[2];
+       double d[2];
+       double h[2];
+
+       t_albe_exchange exchange;       /* exchange between 2bp and 3bp calc */
+} t_albe_mult_params;
+
+/* function prototypes */
+int albe_mult_complete_params(t_albe_mult_params *p);
+int albe_mult_1bp(t_moldyn *moldyn,t_atom *ai);
+int albe_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+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);
+
+#endif
diff --git a/sic.c b/sic.c
index 3c815bac57a9ce9d0325bb4e43523177e25ad79e..73b2935d1633dd1e2de111cde4f4f4cc91892a35 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -13,6 +13,7 @@
 /* potential */
 #include "potentials/harmonic_oscillator.h"
 #include "potentials/lennard_jones.h"
+#include "potentials/albe.h"
 
 #ifdef TERSOFF_ORIG
 #include "potentials/tersoff_orig.h"
@@ -84,6 +85,7 @@ int main(int argc,char **argv) {
        t_lj_params lj;
        t_ho_params ho;
        t_tersoff_mult_params tp;
+       t_albe_mult_params ap;
 
        /* atom injection counter */
        int inject;
@@ -105,6 +107,12 @@ int main(int argc,char **argv) {
        set_potential3b_j1(&md,tersoff_mult_2bp);
        set_potential3b_k1(&md,tersoff_mult_3bp);
        set_potential3b_j2(&md,tersoff_mult_post_2bp);
+#elif ALBE
+       set_potential1b(&md,albe_mult_1bp);
+       set_potential3b_j1(&md,albe_mult_3bp_j1);
+       set_potential3b_k1(&md,albe_mult_3bp_k1);
+       set_potential3b_j2(&md,albe_mult_3bp_j2);
+       set_potential3b_k2(&md,albe_mult_3bp_k2);
 #else
        set_potential3b_j1(&md,tersoff_mult_3bp_j1);
        set_potential3b_k1(&md,tersoff_mult_3bp_k1);
@@ -113,12 +121,18 @@ int main(int argc,char **argv) {
 #endif
        //set_potential2b(&md,lennard_jones);
        //set_potential2b(&md,harmonic_oscillator);
+
+#ifdef ALBE
+       set_potential_params(&md,&ap);
+#else
        set_potential_params(&md,&tp);
+#endif
        //set_potential_params(&md,&lj);
        //set_potential_params(&md,&ho);
 
        /* cutoff radius */
-       set_cutoff(&md,TM_S_SI);
+       set_cutoff(&md,ALBE_S_SI);
+       //set_cutoff(&md,TM_S_SI);
        //set_cutoff(&md,LC_SI*sqrt(3.0));
        //set_cutoff(&md,2.0*LC_SI);
 
@@ -169,9 +183,47 @@ int main(int argc,char **argv) {
 
        tersoff_mult_complete_params(&tp);
 
+       /*
+         * albe mult potential parameters for SiC
+        */
+       ap.S[0]=ALBE_S_SI;
+       ap.R[0]=ALBE_R_SI;
+       ap.A[0]=ALBE_A_SI;
+       ap.B[0]=ALBE_B_SI;
+       ap.r0[0]=ALBE_R0_SI;
+       ap.lambda[0]=ALBE_LAMBDA_SI;
+       ap.mu[0]=ALBE_MU_SI;
+       ap.gamma[0]=ALBE_GAMMA_SI;
+       ap.c[0]=ALBE_C_SI;
+       ap.d[0]=ALBE_D_SI;
+       ap.h[0]=ALBE_H_SI;
+
+       ap.S[1]=ALBE_S_C;
+       ap.R[1]=ALBE_R_C;
+       ap.A[1]=ALBE_A_C;
+       ap.B[1]=ALBE_B_C;
+       ap.r0[1]=ALBE_R0_C;
+       ap.lambda[1]=ALBE_LAMBDA_C;
+       ap.mu[1]=ALBE_MU_C;
+       ap.gamma[1]=ALBE_GAMMA_C;
+       ap.c[1]=ALBE_C_C;
+       ap.d[1]=ALBE_D_C;
+       ap.h[1]=ALBE_H_C;
+
+       ap.Smixed=ALBE_S_SIC;
+       ap.Rmixed=ALBE_R_SIC;
+       ap.Amixed=ALBE_A_SIC;
+       ap.Bmixed=ALBE_B_SIC;
+       ap.r0_mixed=ALBE_R0_SIC;
+       ap.lambda_m=ALBE_LAMBDA_SIC;
+       ap.mu_m=ALBE_MU_SIC;
+
+       albe_mult_complete_params(&ap);
+
        /* set (initial) dimensions of simulation volume */
-       //set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
-       set_dim(&md,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,TRUE);
+       set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
+       //set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE);
+       //set_dim(&md,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,6*TM_LC_3C_SIC,TRUE);
 
        /* set periodic boundary conditions in all directions */
        set_pbc(&md,TRUE,TRUE,TRUE);
@@ -179,21 +231,23 @@ int main(int argc,char **argv) {
        /* create the lattice / place atoms */
        //create_lattice(&md,CUBIC,LC_SI,SI,M_SI,
        //create_lattice(&md,FCC,LC_SI,SI,M_SI,
-       //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+       create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
        //create_lattice(&md,DIAMOND,LC_C,C,M_C,
-       //               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
        //               ATOM_ATTR_2BP|ATOM_ATTR_HB,
+                      0,6,6,6,NULL);
        //               1,6,6,6,NULL);
 
        /* create centered zinc blende lattice */
-       r.x=0.5*0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x;
-       create_lattice(&md,FCC,TM_LC_3C_SIC,SI,M_SI,
-                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
-                      0,6,6,6,&r);
-       r.x+=0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x;
-       create_lattice(&md,FCC,TM_LC_3C_SIC,C,M_C,
-                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
-                      1,6,6,6,&r);
+       //r.x=0.5*0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x;
+       //create_lattice(&md,FCC,TM_LC_3C_SIC,SI,M_SI,
+       //               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+       //               0,6,6,6,&r);
+       //r.x+=0.25*TM_LC_3C_SIC; r.y=r.x; r.z=r.x;
+       //create_lattice(&md,FCC,TM_LC_3C_SIC,C,M_C,
+       //               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+       //               1,6,6,6,&r);
+
        moldyn_bc_check(&md);
 
        /* testing configuration */
@@ -243,7 +297,7 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        /* initial configuration */
-       moldyn_add_schedule(&md,1000,1.0);
+       moldyn_add_schedule(&md,100,1.0);
        /* adding atoms */
        //for(inject=0;inject<INJECT;inject++) {
        //      /* injecting atom and run with enabled t scaling */
@@ -261,8 +315,8 @@ int main(int argc,char **argv) {
        moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
        moldyn_set_log(&md,LOG_TEMPERATURE,1);
        moldyn_set_log(&md,LOG_PRESSURE,1);
-       moldyn_set_log(&md,VISUAL_STEP,10);
-       moldyn_set_log(&md,SAVE_STEP,10);
+       moldyn_set_log(&md,VISUAL_STEP,1);
+       moldyn_set_log(&md,SAVE_STEP,1);
        moldyn_set_log(&md,CREATE_REPORT,0);
 
        /*