From a70de3dccbf0a01c39c2643818ec86c0b465caba Mon Sep 17 00:00:00 2001 From: hackbard Date: Tue, 22 May 2007 15:39:26 +0000 Subject: [PATCH] added albe potential (still wrong energy!) --- Makefile | 2 + moldyn.h | 61 ++++-- potentials/albe.c | 468 ++++++++++++++++++++++++++++++++++++++++++++++ potentials/albe.h | 84 +++++++++ sic.c | 86 +++++++-- 5 files changed, 674 insertions(+), 27 deletions(-) create mode 100644 potentials/albe.c create mode 100644 potentials/albe.h diff --git a/Makefile b/Makefile index ebad7e3..bb520c4 100644 --- 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 diff --git a/moldyn.h b/moldyn.h index 3bc40c0..86444f4 100644 --- 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 index 0000000..c2acf2f --- /dev/null +++ b/potentials/albe.c @@ -0,0 +1,468 @@ +/* + * albe.c - albe potential + * + * author: Frank Zirkelbach + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include +#include +#include + +#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_ijzeta_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(ajdist_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(ajf),&(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(ajkcount++; + + return 0; + +} diff --git a/potentials/albe.h b/potentials/albe.h new file mode 100644 index 0000000..35a0762 --- /dev/null +++ b/potentials/albe.h @@ -0,0 +1,84 @@ +/* + * albe.h - albe potential header file + * + * author: Frank Zirkelbach + * + */ + +#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 3c815ba..73b2935 100644 --- 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