X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Ftersoff.c;h=3a0ed1846967e159430cde751d250e2c5347a901;hb=52485a28c5f9812c0f7a2d4e0094cd8000167b24;hp=a4a87e34d38a574e4e1dcb181d27f80eb7fe09e8;hpb=57847266c05bc38218bf162efdb08e8dd40894d7;p=physik%2Fposic.git diff --git a/potentials/tersoff.c b/potentials/tersoff.c index a4a87e3..3a0ed18 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,7 @@ 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"); @@ -156,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]))) { @@ -170,10 +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); /* 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; } @@ -312,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; @@ -340,6 +420,7 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { unsigned char brand; double ni,tmp; double S,R,s_r,arg; + double energy; params=moldyn->pot_params; exchange=&(params->exchange); @@ -414,17 +495,19 @@ 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)); - /* 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+f_r); + 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; @@ -525,9 +608,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); @@ -543,7 +623,7 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, #endif /* virial */ - //v3_scale(&force,&force,-1.0); + v3_scale(&force,&force,-1.0); virial_calc(ai,&force,&dist_ij); /* derivative wrt k */ @@ -564,7 +644,7 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, #endif /* virial */ - //v3_scale(&force,&force,-1.0); + v3_scale(&force,&force,-1.0); virial_calc(ai,&force,&dist_ik); /* increase k counter */ @@ -573,3 +653,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; +} +