From 4c2140b0f76fb191bdd9b9c2a329877eb0aae531 Mon Sep 17 00:00:00 2001 From: hackbard Date: Tue, 15 Apr 2008 13:06:57 +0200 Subject: [PATCH] fixed bond analyze, introduced more comfortable set potential function, adapted other functions --- bond_analyze.c | 12 +++- moldyn.c | 124 ++++++++++++++++++++-------------------- moldyn.h | 89 ++-------------------------- pair_correlation_calc.c | 6 +- potentials/albe.c | 96 ++++++++++++++++++++++++++++++- potentials/albe.h | 47 ++++++++++++++- potentials/tersoff.c | 67 +++++++++++++++++++++- potentials/tersoff.h | 34 ++++++++++- sic.c | 87 ++-------------------------- 9 files changed, 319 insertions(+), 243 deletions(-) diff --git a/bond_analyze.c b/bond_analyze.c index 0feffef..cc1c318 100644 --- a/bond_analyze.c +++ b/bond_analyze.c @@ -15,6 +15,7 @@ //#include #include "moldyn.h" +#include "potentials/albe.h" int usage(char *prog) { @@ -28,7 +29,7 @@ int main(int argc,char **argv) { t_moldyn moldyn; int ret; - double quality; + double quality[2]; if(argc!=2) { usage(argv[0]); @@ -44,9 +45,14 @@ int main(int argc,char **argv) { return ret; } - bond_analyze(&moldyn,&quality); + /* potential setting */ + set_potential(&moldyn,MOLDYN_POTENTIAL_AM); + albe_mult_set_params(&moldyn,SI,C); - printf("[bond analyze] quality = %f\n",quality); + /* analyzing ... */ + bond_analyze(&moldyn,quality); + + printf("[bond analyze] quality = %f | %f\n",quality[0],quality[1]); moldyn_free_save_file(&moldyn); diff --git a/moldyn.c b/moldyn.c index 84cf700..def5079 100644 --- a/moldyn.c +++ b/moldyn.c @@ -20,6 +20,17 @@ #include "moldyn.h" #include "report/report.h" +/* potential includes */ +#include "potentials/harmonic_oscillator.h" +#include "potentials/lennard_jones.h" +#include "potentials/albe.h" +#ifdef TERSOFF_ORIG +#include "potentials/tersoff_orig.h" +#else +#include "potentials/tersoff.h" +#endif + + /* * global variables, pse and atom colors (only needed here) */ @@ -228,58 +239,29 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) { return 0; } -int set_potential1b(t_moldyn *moldyn,pf_func1b func) { - - moldyn->func1b=func; - - return 0; -} - -int set_potential2b(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func2b=func; - - return 0; -} - -int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j1=func; - - return 0; -} - -int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j2=func; - - return 0; -} - -int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) { - - moldyn->func3b_j3=func; - - return 0; -} - -int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) { - - moldyn->func3b_k1=func; - - return 0; -} - -int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) { - - moldyn->func3b_k2=func; - - return 0; -} - -int set_potential_params(t_moldyn *moldyn,void *params) { +int set_potential(t_moldyn *moldyn,u8 type) { - moldyn->pot_params=params; + switch(type) { + case MOLDYN_POTENTIAL_TM: + moldyn->func1b=tersoff_mult_1bp; + moldyn->func3b_j1=tersoff_mult_3bp_j1; + moldyn->func3b_k1=tersoff_mult_3bp_k1; + moldyn->func3b_j2=tersoff_mult_3bp_j2; + moldyn->func3b_k2=tersoff_mult_3bp_k2; + // missing: check 2b bond func + break; + case MOLDYN_POTENTIAL_AM: + moldyn->func3b_j1=albe_mult_3bp_j1; + moldyn->func3b_k1=albe_mult_3bp_k1; + moldyn->func3b_j2=albe_mult_3bp_j2; + moldyn->func3b_k2=albe_mult_3bp_k2; + moldyn->check_2b_bond=albe_mult_check_2b_bond; + break; + default: + printf("[moldyn] set potential: unknown type %02x\n", + type); + return -1; + } return 0; } @@ -2508,8 +2490,8 @@ int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom, if(bc) check_per_bound(moldyn,&dist); d=v3_absolute_square(&dist); - /* ignore if greater or equal 2 times cutoff */ - if(d>=4.0*moldyn->cutoff_square) + /* ignore if greater cutoff */ + if(d>moldyn->cutoff_square) return 0; /* fill the slots */ @@ -2548,7 +2530,7 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { int i; pcc.dr=dr; - pcc.o1=2.0*moldyn->cutoff/dr; + pcc.o1=moldyn->cutoff/dr; pcc.o2=2*pcc.o1; if(pcc.o1*dr<=moldyn->cutoff) @@ -2612,9 +2594,15 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, d=v3_absolute_square(&dist); /* ignore if greater or equal cutoff */ - if(d>=moldyn->cutoff_square) + if(d>moldyn->cutoff_square) return 0; + /* check for potential bond */ + if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) + return 0; + + d=sqrt(d); + /* now count this bonding ... */ ba=data; @@ -2641,6 +2629,7 @@ int bond_analyze(t_moldyn *moldyn,double *quality) { // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total int qcnt; + int ccnt,cset; t_ba ba; int i; t_atom *atom; @@ -2660,8 +2649,9 @@ int bond_analyze(t_moldyn *moldyn,double *quality) { memset(ba.bcnt,0,moldyn->count*sizeof(int)); ba.tcnt=0; - qcnt=0; + ccnt=0; + cset=0; atom=moldyn->atom; @@ -2673,17 +2663,25 @@ int bond_analyze(t_moldyn *moldyn,double *quality) { qcnt+=4; } else { - if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) + if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) { qcnt+=4; + ccnt+=1; + } + cset+=1; } } -printf("%d %d\n",qcnt,ba.tcnt); - if(quality) - *quality=1.0*qcnt/ba.tcnt; - else - printf("[moldyn] bond analyze: quality = %f\n", - 1.0*qcnt/ba.tcnt); + printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset); + printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt); + + if(quality) { + quality[0]=1.0*ccnt/cset; + quality[1]=1.0*qcnt/ba.tcnt; + } + else { + printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt); + printf("[moldyn] bond analyze: tot_q=%f\n",1.0*qcnt/ba.tcnt); + } return 0; } diff --git a/moldyn.h b/moldyn.h index 1647c1c..de6157c 100644 --- a/moldyn.h +++ b/moldyn.h @@ -217,6 +217,10 @@ typedef struct s_moldyn { t_random random; /* random interface */ double debug; /* debugging stuff, ignore */ + + /* potential 2 body check function */ + int (*check_2b_bond)(struct s_moldyn *moldyn, + t_atom *itom,t_atom *jtom,u8 bc); } t_moldyn; typedef struct s_pcc { @@ -288,6 +292,7 @@ typedef struct s_ba { #define MOLDYN_POTENTIAL_HO 0x00 #define MOLDYN_POTENTIAL_LJ 0x01 #define MOLDYN_POTENTIAL_TM 0x02 +#define MOLDYN_POTENTIAL_AM 0x03 #define LOG_TOTAL_ENERGY 0x01 #define LOG_TOTAL_MOMENTUM 0x02 @@ -330,77 +335,6 @@ typedef struct s_ba { //#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.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.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.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.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 -#define TM_D_C 4.384 -#define TM_H_C -0.57058 - -#define TM_CHI_SIC 0.9776 - -#define TM_LC_SIC 4.32 /* 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_LC_SI 5.429 - -#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 - -#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 - - /* * lattice types */ @@ -415,10 +349,6 @@ typedef struct s_ba { * */ -typedef int (*pf_func1b)(t_moldyn *,t_atom *); -typedef int (*pf_func2b)(t_moldyn *,t_atom *,t_atom *,u8); -typedef int (*pf_func3b)(t_moldyn *,t_atom *,t_atom *,t_atom *,u8); - int moldyn_init(t_moldyn *moldyn,int argc,char **argv); int moldyn_shutdown(t_moldyn *moldyn); @@ -431,14 +361,7 @@ int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc); int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize); int set_nn_dist(t_moldyn *moldyn,double dist); int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z); -int set_potential1b(t_moldyn *moldyn,pf_func1b func); -int set_potential2b(t_moldyn *moldyn,pf_func2b func); -int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func); -int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func); -int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func); -int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func); -int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func); -int set_potential_params(t_moldyn *moldyn,void *params); +int set_potential(t_moldyn *moldyn,u8 type); int set_avg_skip(t_moldyn *moldyn,int skip); diff --git a/pair_correlation_calc.c b/pair_correlation_calc.c index c7716a3..56369f6 100644 --- a/pair_correlation_calc.c +++ b/pair_correlation_calc.c @@ -48,11 +48,11 @@ int main(int argc,char **argv) { return ret; } - //moldyn.cutoff*=2; - //moldyn.cutoff_square*=4; + moldyn.cutoff*=2; + moldyn.cutoff_square*=4; dr=atof(argv[2]); - slots=2.0*moldyn.cutoff/dr; + slots=moldyn.cutoff/dr; printf("[pair corr calc]\n"); printf(" slots: %d\n",slots); printf(" cutoff: %f\n",moldyn.cutoff); diff --git a/potentials/albe.c b/potentials/albe.c index d157bae..ddfb3c6 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -20,14 +20,78 @@ #include "albe.h" /* create mixed terms from parameters and set them */ -int albe_mult_complete_params(t_albe_mult_params *p) { +int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) { - printf("[moldyn] albe parameter completion\n"); + t_albe_mult_params *p; + + /* 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("[moldyn] albe mult parameter info:\n"); + 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); @@ -450,3 +514,29 @@ if(moldyn->time>DSTART&&moldyn->timer),&(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.h b/potentials/albe.h index 09b96b1..15386e7 100644 --- a/potentials/albe.h +++ b/potentials/albe.h @@ -75,12 +75,57 @@ typedef struct s_albe_mult_params { } t_albe_mult_params; /* function prototypes */ -int albe_mult_complete_params(t_albe_mult_params *p); +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 diff --git a/potentials/tersoff.c b/potentials/tersoff.c index 953b839..b83bfd6 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -20,9 +20,70 @@ #include "tersoff.h" /* create mixed terms from parameters and set them */ -int tersoff_mult_complete_params(t_tersoff_mult_params *p) { +int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2) { - printf("[moldyn] tersoff parameter completion\n"); + t_tersoff_mult_params *p; + + /* alloc mem for potential parameters */ + moldyn->pot_params=malloc(sizeof(t_tersoff_mult_params)); + if(moldyn->pot_params==NULL) { + perror("[tersoff] pot params alloc"); + return -1; + } + + /* these are now tersoff parameters */ + p=moldyn->pot_params; + + // only 1 combination by now :p + switch(element1) { + case SI: + /* type: silicon */ + p->S[0]=TM_S_SI; + p->R[0]=TM_R_SI; + p->A[0]=TM_A_SI; + p->B[0]=TM_B_SI; + p->lambda[0]=TM_LAMBDA_SI; + p->mu[0]=TM_MU_SI; + p->beta[0]=TM_BETA_SI; + p->n[0]=TM_N_SI; + p->c[0]=TM_C_SI; + p->d[0]=TM_D_SI; + p->h[0]=TM_H_SI; + switch(element2) { + case C: + p->chi=TM_CHI_SIC; + break; + default: + printf("[tersoff] WARNING: element2\n"); + return -1; + } + break; + default: + printf("[tersoff] WARNING: element1\n"); + return -1; + } + + switch(element2) { + case C: + /* type carbon */ + p->S[1]=TM_S_C; + p->R[1]=TM_R_C; + p->A[1]=TM_A_C; + p->B[1]=TM_B_C; + p->lambda[1]=TM_LAMBDA_C; + p->mu[1]=TM_MU_C; + p->beta[1]=TM_BETA_C; + p->n[1]=TM_N_C; + p->c[1]=TM_C_C; + p->d[1]=TM_D_C; + p->h[1]=TM_H_C; + break; + default: + printf("[tersoff] WARNING: element1\n"); + return -1; + } + + printf("[tersoff] parameter completion\n"); p->S2[0]=p->S[0]*p->S[0]; p->S2[1]=p->S[1]*p->S[1]; p->Smixed=sqrt(p->S[0]*p->S[1]); @@ -33,7 +94,7 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) { p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]); p->mu_m=0.5*(p->mu[0]+p->mu[1]); - printf("[moldyn] tersoff mult parameter info:\n"); + printf("[tersoff] 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); diff --git a/potentials/tersoff.h b/potentials/tersoff.h index 52e848f..4e12957 100644 --- a/potentials/tersoff.h +++ b/potentials/tersoff.h @@ -74,7 +74,7 @@ typedef struct s_tersoff_mult_params { } t_tersoff_mult_params; /* function prototypes */ -int tersoff_mult_complete_params(t_tersoff_mult_params *p); +int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2); int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai); int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int tersoff_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); @@ -84,4 +84,36 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int tersoff_mult_3bp_k2(t_moldyn *moldyn, t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); +/* tersoff potential paramter defines */ + +// silicon +#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.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 + +// carbon +#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.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 +#define TM_D_C 4.384 +#define TM_H_C -0.57058 + +// mixed: silicon carbide +#define TM_CHI_SIC 0.9776 +#define TM_LC_SIC 4.32 /* A */ + #endif diff --git a/sic.c b/sic.c index dea4d4b..0e14bb1 100644 --- a/sic.c +++ b/sic.c @@ -249,10 +249,6 @@ int main(int argc,char **argv) { /* hook parameter structure */ t_hp hookparam; - /* potential parameters */ - t_tersoff_mult_params tp; - t_albe_mult_params ap; - /* testing location & velocity vector */ t_3dvec r,v; memset(&r,0,sizeof(t_3dvec)); @@ -266,22 +262,9 @@ int main(int argc,char **argv) { /* choose potential */ #ifdef ALBE - 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_potential1b(&md,tersoff_mult_1bp); - set_potential3b_j1(&md,tersoff_mult_3bp_j1); - set_potential3b_k1(&md,tersoff_mult_3bp_k1); - set_potential3b_j2(&md,tersoff_mult_3bp_j2); - set_potential3b_k2(&md,tersoff_mult_3bp_k2); -#endif - -#ifdef ALBE - set_potential_params(&md,&ap); + set_potential(&md,MOLDYN_POTENTIAL_AM); #else - set_potential_params(&md,&tp); + set_potential(&md,MOLDYN_POTENTIAL_TM); #endif /* cutoff radius & bondlen */ @@ -302,74 +285,12 @@ int main(int argc,char **argv) { /* * tersoff mult potential parameters for SiC */ - tp.S[0]=TM_S_SI; - tp.R[0]=TM_R_SI; - tp.A[0]=TM_A_SI; - tp.B[0]=TM_B_SI; - tp.lambda[0]=TM_LAMBDA_SI; - tp.mu[0]=TM_MU_SI; - tp.beta[0]=TM_BETA_SI; - tp.n[0]=TM_N_SI; - tp.c[0]=TM_C_SI; - tp.d[0]=TM_D_SI; - tp.h[0]=TM_H_SI; - - tp.S[1]=TM_S_C; - tp.R[1]=TM_R_C; - tp.A[1]=TM_A_C; - tp.B[1]=TM_B_C; - tp.lambda[1]=TM_LAMBDA_C; - tp.mu[1]=TM_MU_C; - tp.beta[1]=TM_BETA_C; - tp.n[1]=TM_N_C; - tp.c[1]=TM_C_C; - tp.d[1]=TM_D_C; - tp.h[1]=TM_H_C; - - tp.chi=TM_CHI_SIC; - - tersoff_mult_complete_params(&tp); + tersoff_mult_set_params(&md,SI,C); /* * 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; - ap.gamma_m=ALBE_GAMMA_SIC; - ap.c_mixed=ALBE_C_SIC; - ap.d_mixed=ALBE_D_SIC; - ap.h_mixed=ALBE_H_SIC; - - albe_mult_complete_params(&ap); + albe_mult_set_params(&md,SI,C); /* set (initial) dimensions of simulation volume */ #ifdef ALBE -- 2.20.1