From: hackbard Date: Wed, 30 Jul 2008 15:12:21 +0000 (+0200) Subject: introduced albe_orig (much faster!) + small change for c2, d2, c2/d2 ... X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=commitdiff_plain;h=8800ba98083a8e5457b035293999f1c70e780c72 introduced albe_orig (much faster!) + small change for c2, d2, c2/d2 ... --- diff --git a/Makefile b/Makefile index 47670a8..1b19d50 100644 --- a/Makefile +++ b/Makefile @@ -11,8 +11,9 @@ CFLAGS += -DALBE #CFLAGS += -DSTATIC_LISTS +#CFLAGS += -DNDEBUG #CFLAGS += -DDEBUG -#CFLAGS += -DDSTART=19 -DDEND=40 -DDATOM=5832 +#CFLAGS += -DDSTART=-1 -DDEND=3 -DDATOM=0 #CFLAGS += -DVDEBUG LDFLAGS = -lm @@ -20,7 +21,7 @@ LDFLAGS = -lm DEPS = moldyn.o random/random.o list/list.o DEPS += potentials/lennard_jones.o potentials/harmonic_oscillator.o -DEPS += potentials/tersoff.o potentials/albe.o +DEPS += potentials/tersoff.o potentials/albe.o potentials/albe_orig.o ALL = mdrun sic fluctuation_calc postproc pair_correlation_calc diffusion_calc ALL += bond_analyze search_bonds visual_atoms display_atom_data diff --git a/mdrun.c b/mdrun.c index 29d3160..907388a 100644 --- a/mdrun.c +++ b/mdrun.c @@ -11,6 +11,7 @@ #include "potentials/harmonic_oscillator.h" #include "potentials/lennard_jones.h" #include "potentials/albe.h" +#include "potentials/albe_orig.h" #ifdef TERSOFF_ORIG #include "potentials/tersoff_orig.h" #else @@ -196,6 +197,8 @@ int mdrun_parse_config(t_mdrun *mdrun) { if(!strncmp(word[0],"potential",9)) { if(!strncmp(word[1],"albe",4)) mdrun->potential=MOLDYN_POTENTIAL_AM; + if(!strncmp(word[1],"albe_o",6)) + mdrun->potential=MOLDYN_POTENTIAL_AO; if(!strncmp(word[1],"tersoff",7)) mdrun->potential=MOLDYN_POTENTIAL_TM; if(!strncmp(word[1],"ho",2)) @@ -905,6 +908,11 @@ int main(int argc,char **argv) { mdrun.element1, mdrun.element2); break; + case MOLDYN_POTENTIAL_AO: + albe_orig_mult_set_params(&moldyn, + mdrun.element1, + mdrun.element2); + break; case MOLDYN_POTENTIAL_TM: tersoff_mult_set_params(&moldyn, mdrun.element1, diff --git a/moldyn.c b/moldyn.c index 1aff5a8..440a707 100644 --- a/moldyn.c +++ b/moldyn.c @@ -24,6 +24,7 @@ #include "potentials/harmonic_oscillator.h" #include "potentials/lennard_jones.h" #include "potentials/albe.h" +#include "potentials/albe_orig.h" #ifdef TERSOFF_ORIG #include "potentials/tersoff_orig.h" #else @@ -222,15 +223,13 @@ int set_potential(t_moldyn *moldyn,u8 type) { moldyn->func_j1_k1=tersoff_mult_3bp_k2; moldyn->check_2b_bond=tersoff_mult_check_2b_bond; 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; + case MOLDYN_POTENTIAL_AO: + moldyn->func_j1=albe_orig_mult_3bp_j1; + moldyn->func_j1_k0=albe_orig_mult_3bp_k1; + moldyn->func_j1c=albe_orig_mult_3bp_j2; + moldyn->func_j1_k1=albe_orig_mult_3bp_k2; + moldyn->check_2b_bond=albe_orig_mult_check_2b_bond; break; - */ case MOLDYN_POTENTIAL_AM: moldyn->func_i0=albe_mult_i0; moldyn->func_j0=albe_mult_i0_j0; diff --git a/moldyn.h b/moldyn.h index 066ad36..b2d6a34 100644 --- a/moldyn.h +++ b/moldyn.h @@ -307,6 +307,7 @@ typedef struct s_vb { #define MOLDYN_POTENTIAL_LJ 0x01 #define MOLDYN_POTENTIAL_TM 0x02 #define MOLDYN_POTENTIAL_AM 0x03 +#define MOLDYN_POTENTIAL_AO 0x04 #define LOG_TOTAL_ENERGY 0x01 #define LOG_TOTAL_MOMENTUM 0x02 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 * @@ -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->timepot_params; exchange=&(params->exchange); kcount=exchange->kcount; @@ -521,9 +542,10 @@ if(moldyn->time>DSTART&&moldyn->time * */ -#ifndef ALBE_H -#define ALBE_H +#ifndef ALBE_ORIG_H +#define ALBE_ORIG_H -#define ALBE_MAXN 16*27 +/* albe constants */ +#include "albe.h" + +#define ALBE_ORIG_MAXN 16*27 /* albe exchange type */ -typedef struct s_albe_exchange { +typedef struct s_albe_orig_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]; + t_3dvec dist_ik[ALBE_ORIG_MAXN]; + double d_ik2[ALBE_ORIG_MAXN]; + double d_ik[ALBE_ORIG_MAXN]; - double f_c_ik[ALBE_MAXN]; - double df_c_ik[ALBE_MAXN]; + double f_c_ik[ALBE_ORIG_MAXN]; + double df_c_ik[ALBE_ORIG_MAXN]; - double g[ALBE_MAXN]; - double dg[ALBE_MAXN]; - double cos_theta[ALBE_MAXN]; + double g[ALBE_ORIG_MAXN]; + double dg[ALBE_ORIG_MAXN]; + double cos_theta[ALBE_ORIG_MAXN]; double *gamma_i; double *c_i; double *d_i; double *h_i; - double ci2; - double di2; - double ci2di2; + double *ci2; + double *di2; + double *ci2di2; double zeta_ij; double pre_dzeta; int kcount; -} t_albe_exchange; +} t_albe_orig_exchange; /* albe mult (2!) potential parameters */ -typedef struct s_albe_mult_params { +typedef struct s_albe_orig_mult_params { double S[2]; /* albe cutoff radii */ double S2[2]; /* albe cutoff radii squared */ double R[2]; /* albe cutoff radii */ @@ -65,67 +68,33 @@ typedef struct s_albe_mult_params { double gamma[2]; double gamma_m; double c[2]; + double c2[2]; double c_mixed; + double c2_m; double d[2]; + double d2[2]; double d_mixed; + double d2_m; double h[2]; double h_mixed; + double c2d2[2]; + double c2d2_m; - t_albe_exchange exchange; /* exchange between 2bp and 3bp calc */ -} t_albe_mult_params; + t_albe_orig_exchange exchange; /* exchange between 2bp and 3bp calc */ +} t_albe_orig_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); +int albe_orig_mult_set_params(t_moldyn *moldyn,int element1,int elemnt2); +int albe_orig_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int albe_orig_mult_3bp_k1(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); +int albe_orig_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int albe_orig_mult_3bp_k2(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); +int albe_orig_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 +// -> albe.h #endif