X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Falbe.c;h=1ec3938844aa80bcb0043be78133c0ccfd08751a;hb=0d2f9a11030dff3583104dac5d4dcb9f040a1327;hp=677d6e3b6315d48226a4b3da950378b1c0fa865a;hpb=da2d9866e05b1b7a408ecda2b1695e07c30b0533;p=physik%2Fposic.git diff --git a/potentials/albe.c b/potentials/albe.c index 677d6e3..1ec3938 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -20,14 +20,93 @@ #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; + + // set cutoff before parameters (actually only necessary for some pots) + if(moldyn->cutoff==0.0) { + printf("[albe] 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"); + 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"); + p->c2[0]=p->c[0]*p->c[0]; + p->c2[1]=p->c[1]*p->c[1]; + p->c2_mixed=p->c_mixed*p->c_mixed; + p->d2[0]=p->d[0]*p->d[0]; + p->d2[1]=p->d[1]*p->d[1]; + p->d2_mixed=p->d_mixed*p->d_mixed; + p->c2d2[0]=p->c2[0]/p->d2[0]; + p->c2d2[1]=p->c2[1]/p->d2[1]; + p->c2d2_m=p->c2_mixed/p->d2_mixed; + + 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); @@ -191,6 +270,11 @@ int albe_mult_3bp_k1(t_moldyn *moldyn, exchange->zeta_ij+=f_c_ik*g; } +#ifdef DEBUG + if(ai==&(moldyn->atom[DATOM])) + printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik); +#endif + /* store even more data for second k loop */ exchange->g[kcount]=g; exchange->dg[kcount]=dg; @@ -275,31 +359,32 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { db=-0.5*b/(1.0+exchange->zeta_ij); } - /* force contribution */ + /* force contribution for atom i */ scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism 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 + + /* force contribution for atom j */ + v3_scale(&force,&force,-1.0); // dri rij = - drj rij + v3_add(&(aj->f),&(aj->f),&force); + + /* virial */ + virial_calc(ai,&force,&(exchange->dist_ij)); #ifdef DEBUG -if(moldyn->time>DSTART&&moldyn->timeatom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) { 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])) + if(ai==&(moldyn->atom[DATOM])) printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); - if(aj==&(moldyn->atom[0])) + if(aj==&(moldyn->atom[DATOM])) 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); printf(" %f %f %f\n",exchange->zeta_ij,.0,.0); } -} #endif - /* virial */ - virial_calc(ai,&force,&(exchange->dist_ij)); - /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ exchange->pre_dzeta=0.5*f_a*f_c*db; @@ -395,7 +480,6 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, v3_add(&(aj->f),&(aj->f),&force); #ifdef DEBUG -if(moldyn->time>DSTART&&moldyn->timeatom[DATOM])) { 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); @@ -403,16 +487,15 @@ if(moldyn->time>DSTART&&moldyn->timef),&(ai->f),&force); - /* virial */ - 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); @@ -423,7 +506,6 @@ if(moldyn->time>DSTART&&moldyn->timef),&(ak->f),&force); #ifdef DEBUG -if(moldyn->time>DSTART&&moldyn->timeatom[DATOM])) { 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); @@ -431,19 +513,44 @@ if(moldyn->time>DSTART&&moldyn->timef),&(ai->f),&force); - /* virial */ - virial_calc(ai,&force,&dist_ik); - /* increase k counter */ exchange->kcount++; return 0; } + +int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) { + + t_albe_mult_params *params; + t_3dvec dist; + double d; + u8 brand; + + v3_sub(&dist,&(jtom->r),&(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; +}