X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Ftersoff.c;h=5039ff7dd76a6ceb11108172fa61408f6abb5242;hb=1965279af348fbb9b69459d01516bbcd52b6f240;hp=953b839bf893a95c2c6469c98fb8be3fa6c43913;hpb=20409ee0c545235c9246edde2d3cda5de0ddabde;p=physik%2Fposic.git diff --git a/potentials/tersoff.c b/potentials/tersoff.c index 953b839..5039ff7 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -20,9 +20,76 @@ #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; + + // set cutoff before parameters (actually only necessary for some pots) + if(moldyn->cutoff==0.0) { + printf("[tersoff] WARNING: no cutoff!\n"); + return -1; + } + + /* 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 +100,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); @@ -157,7 +224,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* add forces */ v3_add(&(ai->f),&(ai->f),&force); - v3_sub(&(aj->f),&(aj->f),&force); // reason: dri rij = - drj rij + v3_scale(&force,&force,-1.0); // reason: dri rij = - drj rij + v3_add(&(aj->f),&(aj->f),&force); #ifdef DEBUG if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) { @@ -171,7 +239,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { #endif /* virial */ - virial_calc(ai,&force,&dist_ij); + virial_calc(aj,&force,&dist_ij); /* energy 2bp contribution */ energy=f_r*f_c; @@ -417,7 +485,8 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { 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 + v3_scale(&force,&force,-1.0); // dri rij = - drj rij + v3_add(&(aj->f),&(aj->f),&force); #ifdef DEBUG if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) { @@ -434,8 +503,7 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { #endif /* virial */ - if(ajdist_ij)); + virial_calc(aj,&force,&(exchange->dist_ij)); /* dzeta prefactor = - 0.5 f_c f_a db */ exchange->pre_dzeta=-0.5*f_a*f_c*db; @@ -544,9 +612,6 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, } #endif - /* virial */ - //virial_calc(ai,&force,&dist_ij); - /* derivative wrt j */ v3_scale(&force,&dcosdrj,fcdg*pre_dzeta); @@ -562,9 +627,8 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, #endif /* virial */ - //v3_scale(&force,&force,-1.0); - if(ajkcount++; @@ -594,3 +657,30 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, return 0; } + +int tersoff_mult_check_2b_bond(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_tersoff_mult_params *params; + t_3dvec dist; + double d; + u8 brand; + + v3_sub(&dist,&(aj->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + + params=moldyn->pot_params; + brand=ai->brand; + + if(brand==aj->brand) { + if(d<=params->S2[brand]) + return TRUE; + } + else { + if(d<=params->S2mixed) + return TRUE; + } + + return FALSE; +} +