X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Ftersoff.c;h=5039ff7dd76a6ceb11108172fa61408f6abb5242;hb=1965279af348fbb9b69459d01516bbcd52b6f240;hp=b6533e19bf90c3709b5915198c6a2964be2efc41;hpb=7fcb896467bdf236f60c4c3784a9d51214737fcb;p=physik%2Fposic.git diff --git a/potentials/tersoff.c b/potentials/tersoff.c index b6533e1..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); @@ -93,6 +160,10 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { int brand; double s_r; double arg; + double energy; + + printf("WARNING! - tersoff_mult_2bp is obsolete.\n"); + printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n"); /* use newtons third law */ if(aif),&(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]))) { @@ -167,17 +239,13 @@ 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); - //ai->virial.xx-=force.x*dist_ij.x; - //ai->virial.yy-=force.y*dist_ij.y; - //ai->virial.zz-=force.z*dist_ij.z; - //ai->virial.xy-=force.x*dist_ij.y; - //ai->virial.xz-=force.x*dist_ij.z; - //ai->virial.yz-=force.y*dist_ij.z; + virial_calc(aj,&force,&dist_ij); /* energy 2bp contribution */ - moldyn->energy+=f_r*f_c; + energy=f_r*f_c; + moldyn->energy+=energy; + ai->e+=0.5*energy; + aj->e+=0.5*energy; return 0; } @@ -316,6 +384,14 @@ int tersoff_mult_3bp_k1(t_moldyn *moldyn, exchange->zeta_ij+=f_c_ik*g; } +#ifdef DEBUG + if((ai==&(moldyn->atom[0]))| + (aj==&(moldyn->atom[864]))| + (ak==&(moldyn->atom[1003]))) { + printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos); + } +#endif + /* store even more data for second k loop */ exchange->g[kcount]=g; exchange->dg[kcount]=dg; @@ -336,11 +412,15 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_tersoff_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,chi; + double lambda,A; double d_ij; unsigned char brand; double ni,tmp; double S,R,s_r,arg; + double energy; params=moldyn->pot_params; exchange=&(params->exchange); @@ -350,14 +430,18 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { S=params->S[brand]; R=params->R[brand]; B=params->B[brand]; + A=params->A[brand]; mu=params->mu[brand]; + lambda=params->lambda[brand]; chi=1.0; } else { S=params->Smixed; R=params->Rmixed; B=params->Bmixed; + A=params->Amixed; mu=params->mu_m; + lambda=params->lambda_m; chi=params->chi; } @@ -379,6 +463,10 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { f_a=-B*exp(-mu*d_ij); df_a=mu*f_a/d_ij; + /* f_r, df_r */ + f_r=A*exp(-lambda*d_ij); + df_r=lambda*f_r/d_ij; + /* b, db */ if(exchange->zeta_ij==0.0) { b=chi; @@ -394,10 +482,11 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { } /* force contribution */ - v3_scale(&force,&(exchange->dist_ij),df_a*f_c+f_a*df_c); - v3_scale(&force,&force,-0.5*b); + 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]))) { @@ -407,18 +496,22 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { 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); + printf(" %f %f %f\n",exchange->zeta_ij,.0,.0); } #endif /* virial */ - //virial_calc(ai,&force,&(exchange->dist_ij)); - //virial_calc(aj,&force,&(exchange->dist_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; /* energy contribution */ - moldyn->energy+=0.5*f_c*b*f_a; + energy=0.5*f_c*(b*f_a+f_r); + moldyn->energy+=energy; + ai->e+=energy; /* reset k counter for second k loop */ exchange->kcount=0; @@ -519,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); @@ -537,10 +627,11 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, #endif /* virial */ - //virial_calc(aj,&force,&dist_ij); + v3_scale(&force,&force,-1.0); + virial_calc(ai,&force,&dist_ij); /* derivative wrt k */ - v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rij = - drj rij + v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik v3_scale(&tmp,&dcosdrk,fcdg); v3_add(&force,&force,&tmp); v3_scale(&force,&force,pre_dzeta); @@ -557,7 +648,8 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, #endif /* virial */ - virial_calc(ak,&force,&dist_ik); + v3_scale(&force,&force,-1.0); + virial_calc(ai,&force,&dist_ik); /* increase k counter */ exchange->kcount++; @@ -565,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; +} +