From: hackbard Date: Mon, 20 Feb 2012 15:55:56 +0000 (+0100) Subject: Merge branch 'leadoff' X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=commitdiff_plain;h=6e6d7126ea9a845f11637d8e1b8eb2b570ac4dc9;hp=98df942a2412e6ace72ae221f2e967f24bba7108 Merge branch 'leadoff' Conflicts: Makefile moldyn.c potentials/albe.c potentials/albe.h runmd ... hopefully fixed! Just a test to merge leadoff back to master! :) --- diff --git a/mdrun.c b/mdrun.c index fb1cbe0..26fbc87 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 @@ -234,6 +235,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)) @@ -1653,6 +1656,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 33224a6..1d6e0b7 100644 --- a/moldyn.c +++ b/moldyn.c @@ -34,6 +34,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 @@ -270,19 +271,29 @@ int set_potential(t_moldyn *moldyn,u8 type) { moldyn->func3b_k2=tersoff_mult_3bp_k2; moldyn->check_2b_bond=tersoff_mult_check_2b_bond; break; + 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->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->func_i0=albe_mult_i0; + moldyn->func_j0=albe_mult_i0_j0; + moldyn->func_j0_k0=albe_mult_i0_j0_k0; + moldyn->func_j0e=albe_mult_i0_j1; + moldyn->func_j1=albe_mult_i0_j2; + moldyn->func_j1_k0=albe_mult_i0_j2_k0; + moldyn->func_j1c=albe_mult_i0_j3; moldyn->check_2b_bond=albe_mult_check_2b_bond; break; case MOLDYN_POTENTIAL_HO: - moldyn->func2b=harmonic_oscillator; + moldyn->func_j0=harmonic_oscillator; moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond; break; case MOLDYN_POTENTIAL_LJ: - moldyn->func2b=lennard_jones; + moldyn->func_j0=lennard_jones; moldyn->check_2b_bond=lennard_jones_check_2b_bond; break; default: @@ -2487,8 +2498,8 @@ int potential_force_calc(t_moldyn *moldyn) { /* single particle potential/force */ if(itom[i].attr&ATOM_ATTR_1BP) - if(moldyn->func1b) - moldyn->func1b(moldyn,&(itom[i])); + if(moldyn->func_i0) + moldyn->func_i0(moldyn,&(itom[i])); if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP))) continue; @@ -2503,8 +2514,14 @@ int potential_force_calc(t_moldyn *moldyn) { dnlc=lc->dnlc; +#ifndef STATIC_LISTS + /* check for later 3 body interaction */ + if(itom[i].attr&ATOM_ATTR_3BP) + memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); +#endif + /* first loop over atoms j */ - if(moldyn->func2b) { + if(moldyn->func_j0) { for(j=0;j<27;j++) { bc_ij=(jrun3bp=1; + if((jtom->attr&ATOM_ATTR_2BP)& (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); + moldyn->func_j0(moldyn, + &(itom[i]), + jtom, + bc_ij); + } + + /* 3 body potential/force */ + + /* in j loop, 3bp run can be skipped */ + if(!(moldyn->run3bp)) + continue; + + if(!(itom[i].attr&ATOM_ATTR_3BP)) + continue; + + if(!(jtom->attr&ATOM_ATTR_3BP)) + continue; + + if(moldyn->func_j0_k0==NULL) + continue; + + /* first loop over atoms k in first j loop */ + for(k=0;k<27;k++) { + + bc_ik=(kstart==NULL) + continue; + + do { + ktom=that->current->data; +#endif + + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; + + //if(ktom==jtom) + // continue; + + if(ktom==&(itom[i])) + continue; + + moldyn->func_j0_k0(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); +#ifdef STATIC_LISTS } #ifdef STATIC_LISTS } @@ -2554,7 +2628,7 @@ int potential_force_calc(t_moldyn *moldyn) { } } - /* 3 body potential/force */ + /* continued 3 body potential/force */ if(!(itom[i].attr&ATOM_ATTR_3BP)) continue; @@ -2607,18 +2681,18 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset 3bp run */ moldyn->run3bp=1; - if(moldyn->func3b_j1) - moldyn->func3b_j1(moldyn, - &(itom[i]), - jtom, - bc_ij); + if(moldyn->func_j1) + moldyn->func_j1(moldyn, + &(itom[i]), + jtom, + bc_ij); - /* in first j loop, 3bp run can be skipped */ + /* in j loop, 3bp run can be skipped */ if(!(moldyn->run3bp)) continue; - /* first loop over atoms k */ - if(moldyn->func3b_k1) { + /* first loop over atoms k in second j loop */ + if(moldyn->func_j1_k0) { for(k=0;k<27;k++) { @@ -2651,8 +2725,8 @@ int potential_force_calc(t_moldyn *moldyn) { if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; @@ -2676,14 +2750,15 @@ int potential_force_calc(t_moldyn *moldyn) { } - if(moldyn->func3b_j2) - moldyn->func3b_j2(moldyn, - &(itom[i]), - jtom, - bc_ij); + /* continued j loop after first k loop */ + if(moldyn->func_j1c) + moldyn->func_j1c(moldyn, + &(itom[i]), + jtom, + bc_ij); /* second loop over atoms k */ - if(moldyn->func3b_k2) { + if(moldyn->func_j1_k1) { for(k=0;k<27;k++) { @@ -2716,17 +2791,17 @@ int potential_force_calc(t_moldyn *moldyn) { if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; - moldyn->func3b_k2(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func_j1_k1(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); #ifdef STATIC_LISTS } @@ -2741,11 +2816,11 @@ int potential_force_calc(t_moldyn *moldyn) { } - /* 2bp post function */ - if(moldyn->func3b_j3) { - moldyn->func3b_j3(moldyn, - &(itom[i]), - jtom,bc_ij); + /* finish of j loop after second k loop */ + if(moldyn->func_j1e) { + moldyn->func_j1e(moldyn, + &(itom[i]), + jtom,bc_ij); } #ifdef STATIC_LISTS } @@ -2756,7 +2831,7 @@ int potential_force_calc(t_moldyn *moldyn) { #endif } - + #ifdef DEBUG //printf("\n\n"); #endif diff --git a/moldyn.h b/moldyn.h index a1866a9..0670f3d 100644 --- a/moldyn.h +++ b/moldyn.h @@ -113,15 +113,18 @@ typedef struct s_moldyn { double volume; /* volume of sim cell (dim.x*dim.y*dim.z) */ /* potential force function and parameter pointers */ - int (*func1b)(struct s_moldyn *moldyn,t_atom *ai); - int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j2)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j3)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_k1)(struct s_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); - int (*func3b_k2)(struct s_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_i0)(struct s_moldyn *moldyn,t_atom *ai); + int (*func_j0)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j0_k0)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j0e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1_k0)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j1c)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1_k1)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j1e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); void *pot_params; unsigned char run3bp; @@ -358,6 +361,7 @@ typedef struct s_offset_params { #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.c b/potentials/albe.c index 1ec3938..9f98282 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -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,7 +85,10 @@ 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_mixed=p->c_mixed*p->c_mixed; p->d_mixed=ALBE_D_SIC; + p->d2_mixed=p->d_mixed*p->d_mixed; + p->c2d2_m=p->c2_mixed/p->d2_mixed; p->h_mixed=ALBE_H_SIC; break; default: @@ -114,34 +123,56 @@ 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(" c2 | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed); + printf(" d2 | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed); + printf(" c2d2 | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m); + 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) { +/* first i loop */ +int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) { t_albe_mult_params *params; t_albe_exchange *exchange; - unsigned char brand; - double S2; - t_3dvec dist_ij; - double d_ij2,d_ij; + + int i; params=moldyn->pot_params; exchange=&(params->exchange); - /* reset zeta sum */ - exchange->zeta_ij=0.0; + /* zero exchange values */ + memset(exchange->zeta,0,ALBE_MAXN*sizeof(double)); + for(i=0;idzeta[i],0,ALBE_MAXN*sizeof(t_3dvec)); + exchange->jcnt=0; + exchange->j2cnt=0; - /* - * set ij depending values - */ + return 0; +} +/* first j loop within first i loop */ +int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_albe_mult_params *params; + t_albe_exchange *exchange; + + double S2,S,R,d2,d,s_r,arg; + t_3dvec dist; + int j; + u8 brand; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* get j counter */ + j=exchange->jcnt; + + /* set ij depending values */ brand=ai->brand; if(brand==aj->brand) { S2=params->S2[brand]; @@ -151,123 +182,132 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { } /* dist_ij, d_ij2 */ - v3_sub(&dist_ij,&(aj->r),&(ai->r)); - if(bc) check_per_bound(moldyn,&dist_ij); - d_ij2=v3_absolute_square(&dist_ij); + v3_sub(&dist,&(aj->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&dist); + exchange->dist[j]=dist; + d2=v3_absolute_square(&dist); + exchange->d2[j]=d2; /* if d_ij2 > S2 => no force & potential energy contribution */ - if(d_ij2>S2) { + if(d2>S2) { moldyn->run3bp=0; + exchange->skip[j]=1; + exchange->jcnt+=1; return 0; } + exchange->skip[j]=0; + + /* more ij depending values */ + if(brand==aj->brand) { + R=params->R[brand]; + S=params->S[brand]; + /* albe needs i,(j/k) depending c,d,h and gamma values */ + exchange->gamma_[j]=&(params->gamma[brand]); + exchange->c_[j]=&(params->c[brand]); + exchange->d_[j]=&(params->d[brand]); + exchange->h_[j]=&(params->h[brand]); + exchange->c2_[j]=&(params->c2[brand]); + exchange->d2_[j]=&(params->d2[brand]); + exchange->c2d2_[j]=&(params->c2d2[brand]); + } + else { + R=params->Rmixed; + S=params->Smixed; + /* albe needs i,(j/k) depending c,d,h and gamma values */ + exchange->gamma_[j]=&(params->gamma_m); + exchange->c_[j]=&(params->c_mixed); + exchange->d_[j]=&(params->d_mixed); + exchange->h_[j]=&(params->h_mixed); + exchange->c2_[j]=&(params->c2_mixed); + exchange->d2_[j]=&(params->d2_mixed); + exchange->c2d2_[j]=&(params->c2d2_m); + } /* d_ij */ - d_ij=sqrt(d_ij2); + d=sqrt(d2); + exchange->d[j]=d; + + /* f_c, df_c */ + if(df_c[j]=1.0; + exchange->df_c[j]=0.0; + } + else { + s_r=S-R; + arg=M_PI*(d-R)/s_r; + exchange->f_c[j]=0.5+0.5*cos(arg); + exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d)); + } - /* store values */ - exchange->dist_ij=dist_ij; - exchange->d_ij2=d_ij2; - exchange->d_ij=d_ij; + /* reset k counter */ + exchange->kcnt=0; - /* reset k counter for first k loop */ - exchange->kcount=0; - return 0; } -/* 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) { +/* first k loop within first j loop within first i loop */ +int albe_mult_i0_j0_k0(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; - unsigned char brand; - double R,S,S2; - t_3dvec dist_ij,dist_ik; - double d_ik2,d_ik,d_ij; - double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg; - double f_c_ik,df_c_ik; - int kcount; + + int j,k; + t_3dvec distj,distk; + double dj,dk,djdk_inv,cos_theta; + double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j; + double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k; + t_3dvec dcosdrj,dcosdrk,tmp; + t_3dvec *dzjj,*dzkk,*dzjk,*dzkj; params=moldyn->pot_params; exchange=&(params->exchange); - kcount=exchange->kcount; - if(kcount>ALBE_MAXN) { - printf("FATAL: neighbours = %d\n",kcount); - printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag); + if(aj==ak) { + exchange->kcnt+=1; + return 0; } - /* ik constants */ - brand=ai->brand; - if(brand==ak->brand) { - R=params->R[brand]; - S=params->S[brand]; - S2=params->S2[brand]; - /* albe needs i,k depending c,d,h and gamma values */ - exchange->gamma_i=&(params->gamma[brand]); - exchange->c_i=&(params->c[brand]); - exchange->d_i=&(params->d[brand]); - exchange->h_i=&(params->h[brand]); - } - else { - R=params->Rmixed; - S=params->Smixed; - S2=params->S2mixed; - /* albe needs i,k depending c,d,h and gamma values */ - exchange->gamma_i=&(params->gamma_m); - exchange->c_i=&(params->c_mixed); - exchange->d_i=&(params->d_mixed); - exchange->h_i=&(params->h_mixed); + /* kjcnt; + k=exchange->kcnt; + if(k>=ALBE_MAXN) { + printf("FATAL: too many neighbours! (%d)\n",k); + printf(" atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag); } - 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)); - if(bc) check_per_bound(moldyn,&dist_ik); - d_ik2=v3_absolute_square(&dist_ik); - - /* store data for second k loop */ - exchange->dist_ik[kcount]=dist_ik; - exchange->d_ik2[kcount]=d_ik2; - - /* return if not within cutoff */ - if(d_ik2>S2) { - exchange->kcount++; + if((k>=j)|(exchange->skip[k])) { + exchange->kcnt+=1; return 0; } - /* d_ik */ - d_ik=sqrt(d_ik2); - - /* dist_ij, d_ij */ - dist_ij=exchange->dist_ij; - d_ij=exchange->d_ij; + /* distances */ + distj=exchange->dist[j]; + distk=exchange->dist[k]; + dj=exchange->d[j]; + dk=exchange->d[k]; + djdk_inv=1.0/(dj*dk); /* cos theta */ - cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); - - /* 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); - dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f.. - - /* zeta sum += f_c_ik * g_ijk */ - if(d_ik<=R) { - exchange->zeta_ij+=g; - f_c_ik=1.0; - df_c_ik=0.0; + cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv; + + /* g(cos(theta)) ij and ik values */ + h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism + d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j); + frac_j=*exchange->c2_[j]/d2_h_cos2_j; + gj=1.0+*exchange->c2d2_[j]-frac_j; + gj*=*(exchange->gamma_[j]); + dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe + if(ak->brand==aj->brand) { + gk=gj; + dgk=dgj; } else { - s_r=S-R; - arg=M_PI*(d_ik-R)/s_r; - f_c_ik=0.5+0.5*cos(arg); - df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik)); - exchange->zeta_ij+=f_c_ik*g; + h_cos_k=*(exchange->h_[k])+cos_theta; + d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k); + frac_k=*exchange->c2_[k]/d2_h_cos2_k; + gk=1.0+*exchange->c2d2_[k]-frac_k; + gk*=*(exchange->gamma_[k]); + dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k; } #ifdef DEBUG @@ -284,33 +324,62 @@ int albe_mult_3bp_k1(t_moldyn *moldyn, exchange->df_c_ik[kcount]=df_c_ik; /* increase k counter */ - exchange->kcount++; + exchange->kcnt+=1; + + return 0; +} + +/* first j loop within first i loop */ +int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_albe_mult_params *params; + t_albe_exchange *exchange; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* increase j counter */ + exchange->jcnt+=1; return 0; } -int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { +/* second j loop within first i loop */ +int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_albe_mult_params *params; t_albe_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; - double lambda,A; - double d_ij,r0; - unsigned char brand; - double S,R,s_r,arg; + + int j; + double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db; + double A,B,mu,lambda,r0; double energy; + t_3dvec *dist,force; + double scale; + u8 brand; params=moldyn->pot_params; exchange=&(params->exchange); + /* get j counter */ + j=exchange->j2cnt; + + /* skip if j not within cutoff */ + if(exchange->skip[j]) { + moldyn->run3bp=0; + exchange->j2cnt+=1; + return 0; + } + + /* distance */ + d=exchange->d[j]; + dist=&(exchange->dist[j]); + f_c=exchange->f_c[j]; + df_c=exchange->df_c[j]; + + /* determine parameters to calculate fa, dfa, fr, dfr */ brand=aj->brand; if(brand==ai->brand) { - S=params->S[brand]; - R=params->R[brand]; B=params->B[brand]; A=params->A[brand]; r0=params->r0[brand]; @@ -318,8 +387,6 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { lambda=params->lambda[brand]; } else { - S=params->Smixed; - R=params->Rmixed; B=params->Bmixed; A=params->Amixed; r0=params->r0_mixed; @@ -327,44 +394,36 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { lambda=params->lambda_m; } - d_ij=exchange->d_ij; - - /* f_c, df_c */ - if(d_ijzeta_ij==0.0) { - b=1.0; - db=0.0; - } - else { - b=1.0/sqrt(1.0+exchange->zeta_ij); - db=-0.5*b/(1.0+exchange->zeta_ij); - } + b=1.0/sqrt(1.0+exchange->zeta[j]); + db=-0.5*b/(1.0+exchange->zeta[j]); + + /* energy contribution */ + energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism + moldyn->energy+=energy; + ai->e+=energy; - /* force contribution for atom i */ + /* force contribution for atom i due to ij bond */ 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_scale(&force,dist,scale); v3_add(&(ai->f),&(ai->f),&force); - /* force contribution for atom j */ +#ifdef NDEBUG +if(ai->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]); +printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z); +} +#endif + + /* force contribution for atom j due to ij bond */ v3_scale(&force,&force,-1.0); // dri rij = - drj rij v3_add(&(aj->f),&(aj->f),&force); @@ -385,57 +444,57 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { } #endif + /* virial */ + virial_calc(ai,&force,dist); + /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ exchange->pre_dzeta=0.5*f_a*f_c*db; - /* energy contribution */ - energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism - moldyn->energy+=energy; - ai->e+=energy; + /* force contribution (drj derivative) */ + v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta); + v3_add(&(aj->f),&(aj->f),&force); + +#ifdef NDEBUG +if(aj->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag); +printf(" t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z); +} +#endif + + /* virial */ + virial_calc(ai,&force,dist); + + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); + +#ifdef NDEBUG +if(ai->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (i contr j der)\n",force.x,force.y,force.z,ai->tag,aj->tag); +printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z); +} +#endif /* reset k counter for second k loop */ - exchange->kcount=0; + exchange->kcnt=0; return 0; } -/* albe 3 body potential function (second k loop) */ -int albe_mult_3bp_k2(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { +/* second k loop within second j loop within first i loop */ +int albe_mult_i0_j2_k0(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; - int kcount; - t_3dvec dist_ik,dist_ij; - double d_ik2,d_ik,d_ij2,d_ij; - unsigned char brand; - double S2; - double g,dg,cos_theta; - double pre_dzeta; - double f_c_ik,df_c_ik; - double dijdik_inv,fcdg,dfcg; - t_3dvec dcosdrj,dcosdrk; - t_3dvec force,tmp; + + int j,k; + t_3dvec force; params=moldyn->pot_params; exchange=&(params->exchange); - kcount=exchange->kcount; - - if(kcount>ALBE_MAXN) - printf("FATAL: neighbours!\n"); - - /* d_ik2 */ - d_ik2=exchange->d_ik2[kcount]; - brand=ak->brand; - if(brand==ai->brand) - S2=params->S2[brand]; - else - S2=params->S2mixed; - - /* return if d_ik > S */ - if(d_ik2>S2) { - exchange->kcount++; + if(aj==ak) { + exchange->kcnt+=1; return 0; } @@ -502,8 +561,8 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, v3_add(&force,&force,&tmp); v3_scale(&force,&force,pre_dzeta); - /* force contribution */ - v3_add(&(ak->f),&(ak->f),&force); + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); #ifdef DEBUG if(ak==&(moldyn->atom[DATOM])) { @@ -523,10 +582,23 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, v3_add(&(ai->f),&(ai->f),&force); /* increase k counter */ - exchange->kcount++; + exchange->kcnt+=1; return 0; +} + +int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + t_albe_mult_params *params; + t_albe_exchange *exchange; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* increase j counter */ + exchange->j2cnt+=1; + + return 0; } int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) { diff --git a/potentials/albe.h b/potentials/albe.h index 89cdae2..9687655 100644 --- a/potentials/albe.h +++ b/potentials/albe.h @@ -8,39 +8,36 @@ #ifndef ALBE_H #define ALBE_H -#define ALBE_MAXN 16*27 +#define ALBE_MAXN (4*27) /* albe exchange type */ typedef struct s_albe_exchange { - t_3dvec dist_ij; - double d_ij2; - double d_ij; + t_3dvec dist[ALBE_MAXN]; + double d2[ALBE_MAXN]; + double d[ALBE_MAXN]; - t_3dvec dist_ik[ALBE_MAXN]; - double d_ik2[ALBE_MAXN]; - double d_ik[ALBE_MAXN]; + double f_c[ALBE_MAXN]; + double df_c[ALBE_MAXN]; - double f_c_ik[ALBE_MAXN]; - double df_c_ik[ALBE_MAXN]; + double zeta[ALBE_MAXN]; + t_3dvec dzeta[ALBE_MAXN][ALBE_MAXN]; - double g[ALBE_MAXN]; - double dg[ALBE_MAXN]; - double cos_theta[ALBE_MAXN]; + u8 skip[ALBE_MAXN]; - double *gamma_i; - double *c_i; - double *d_i; - double *h_i; + double *gamma_[ALBE_MAXN]; + double *c_[ALBE_MAXN]; + double *d_[ALBE_MAXN]; + double *c2_[ALBE_MAXN]; + double *d2_[ALBE_MAXN]; + double *c2d2_[ALBE_MAXN]; + double *h_[ALBE_MAXN]; - double ci2; - double di2; - double ci2di2; - - double zeta_ij; double pre_dzeta; - int kcount; + int jcnt; + int j2cnt; + int kcnt; } t_albe_exchange; /* albe mult (2!) potential parameters */ @@ -65,10 +62,12 @@ typedef struct s_albe_mult_params { double gamma[2]; double gamma_m; double c[2]; + double c2[2]; double c_mixed; double c2[2]; double c2_mixed; double d[2]; + double d2[2]; double d_mixed; double d2[2]; double d2_mixed; diff --git a/potentials/albe_orig.c b/potentials/albe_orig.c new file mode 100644 index 0000000..908f6dd --- /dev/null +++ b/potentials/albe_orig.c @@ -0,0 +1,570 @@ +/* + * albe_orig.c - albe potential + * + * author: Frank Zirkelbach + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../moldyn.h" +#include "../math/math.h" +#include "albe_orig.h" + +/* create mixed terms from parameters and set them */ +int albe_orig_mult_set_params(t_moldyn *moldyn,int element1,int element2) { + + t_albe_orig_mult_params *p; + + // set cutoff before parameters (actually only necessary for some pots) + if(moldyn->cutoff==0.0) { + 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 orig] 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->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: + /* 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->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; + 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->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 orig] WARNING: element2"); + printf("\n"); + return -1; + } + break; + default: + printf("[albe orig] WARNING: element1\n"); + return -1; + } + + 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 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); + printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV); + 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 | %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_orig_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_albe_orig_mult_params *params; + t_albe_orig_exchange *exchange; + unsigned char brand; + double S2; + t_3dvec dist_ij; + double d_ij2,d_ij; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* reset zeta sum */ + exchange->zeta_ij=0.0; + + /* + * set ij depending values + */ + + brand=ai->brand; + if(brand==aj->brand) { + S2=params->S2[brand]; + } + else { + S2=params->S2mixed; + } + + /* dist_ij, d_ij2 */ + v3_sub(&dist_ij,&(aj->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&dist_ij); + d_ij2=v3_absolute_square(&dist_ij); + + /* if d_ij2 > S2 => no force & potential energy contribution */ + if(d_ij2>S2) { + moldyn->run3bp=0; + return 0; + } + + /* d_ij */ + d_ij=sqrt(d_ij2); + + /* store values */ + exchange->dist_ij=dist_ij; + exchange->d_ij2=d_ij2; + exchange->d_ij=d_ij; + + /* reset k counter for first k loop */ + exchange->kcount=0; + + return 0; +} + +/* albe 3 body potential function (first k loop) */ +int albe_orig_mult_3bp_k1(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { + + t_albe_orig_mult_params *params; + t_albe_orig_exchange *exchange; + unsigned char brand; + double R,S,S2; + t_3dvec dist_ij,dist_ik; + double d_ik2,d_ik,d_ij; + double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg; + 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_ORIG_MAXN) { + printf("FATAL: neighbours = %d\n",kcount); + printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag); + } + + /* ik constants */ + brand=ai->brand; + if(brand==ak->brand) { + R=params->R[brand]; + S=params->S[brand]; + S2=params->S2[brand]; + /* albe needs i,k depending c,d,h and gamma values */ + exchange->gamma_i=&(params->gamma[brand]); + 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; + S=params->Smixed; + S2=params->S2mixed; + /* albe needs i,k depending c,d,h and gamma values */ + exchange->gamma_i=&(params->gamma_m); + 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); + } + + /* dist_ik, d_ik2 */ + v3_sub(&dist_ik,&(ak->r),&(ai->r)); + if(bc) check_per_bound(moldyn,&dist_ik); + d_ik2=v3_absolute_square(&dist_ik); + + /* store data for second k loop */ + exchange->dist_ik[kcount]=dist_ik; + exchange->d_ik2[kcount]=d_ik2; + + /* return if not within cutoff */ + if(d_ik2>S2) { + exchange->kcount++; + return 0; + } + + /* d_ik */ + d_ik=sqrt(d_ik2); + + /* dist_ij, d_ij */ + dist_ij=exchange->dist_ij; + d_ij=exchange->d_ij; + + /* cos theta */ + cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); + + /* 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); + dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f.. + + /* zeta sum += f_c_ik * g_ijk */ + if(d_ik<=R) { + exchange->zeta_ij+=g; + f_c_ik=1.0; + df_c_ik=0.0; + } + else { + s_r=S-R; + arg=M_PI*(d_ik-R)/s_r; + f_c_ik=0.5+0.5*cos(arg); + df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik)); + exchange->zeta_ij+=f_c_ik*g; + } + + /* store even more data for second k loop */ + exchange->g[kcount]=g; + exchange->dg[kcount]=dg; + exchange->d_ik[kcount]=d_ik; + exchange->cos_theta[kcount]=cos_theta; + exchange->f_c_ik[kcount]=f_c_ik; + exchange->df_c_ik[kcount]=df_c_ik; + + /* increase k counter */ + exchange->kcount++; + + return 0; +} + +int albe_orig_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + 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; + double scale; + double mu,B; + double lambda,A; + double d_ij,r0; + unsigned char brand; + double S,R,s_r,arg; + double energy; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + brand=aj->brand; + if(brand==ai->brand) { + S=params->S[brand]; + R=params->R[brand]; + B=params->B[brand]; + A=params->A[brand]; + r0=params->r0[brand]; + mu=params->mu[brand]; + lambda=params->lambda[brand]; + } + else { + S=params->Smixed; + R=params->Rmixed; + B=params->Bmixed; + A=params->Amixed; + r0=params->r0_mixed; + mu=params->mu_m; + lambda=params->lambda_m; + } + + d_ij=exchange->d_ij; + + /* f_c, df_c */ + if(d_ijzeta_ij==0.0) { + b=1.0; + db=0.0; + } + else { + b=1.0/sqrt(1.0+exchange->zeta_ij); + db=-0.5*b/(1.0+exchange->zeta_ij); + } + + /* 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); + + /* 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])) + 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 + + /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ + exchange->pre_dzeta=0.5*f_a*f_c*db; + + /* energy contribution */ + energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism + moldyn->energy+=energy; + ai->e+=energy; + + /* reset k counter for second k loop */ + exchange->kcount=0; + + return 0; +} + +/* albe 3 body potential function (second k loop) */ +int albe_orig_mult_3bp_k2(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { + + t_albe_orig_mult_params *params; + t_albe_orig_exchange *exchange; + int kcount; + t_3dvec dist_ik,dist_ij; + double d_ik2,d_ik,d_ij2,d_ij; + unsigned char brand; + double S2; + double g,dg,cos_theta; + double pre_dzeta; + double f_c_ik,df_c_ik; + double dijdik_inv,fcdg,dfcg; + t_3dvec dcosdrj,dcosdrk; + t_3dvec force,tmp; + + /* continue if aj equals ak */ + if(aj==ak) + return 0; + + params=moldyn->pot_params; + exchange=&(params->exchange); + kcount=exchange->kcount; + + if(kcount>ALBE_ORIG_MAXN) + printf("FATAL: neighbours!\n"); + + /* d_ik2 */ + d_ik2=exchange->d_ik2[kcount]; + + brand=ak->brand; + if(brand==ai->brand) + S2=params->S2[brand]; + else + S2=params->S2mixed; + + /* return if d_ik > S */ + if(d_ik2>S2) { + exchange->kcount++; + return 0; + } + + /* prefactor dzeta */ + pre_dzeta=exchange->pre_dzeta; + + /* dist_ik, d_ik */ + dist_ik=exchange->dist_ik[kcount]; + d_ik=exchange->d_ik[kcount]; + + /* f_c_ik, df_c_ik */ + f_c_ik=exchange->f_c_ik[kcount]; + df_c_ik=exchange->df_c_ik[kcount]; + + /* dist_ij, d_ij2, d_ij */ + dist_ij=exchange->dist_ij; + d_ij2=exchange->d_ij2; + d_ij=exchange->d_ij; + + /* g, dg, cos_theta */ + g=exchange->g[kcount]; + dg=exchange->dg[kcount]; + cos_theta=exchange->cos_theta[kcount]; + + /* cos_theta derivatives wrt j,k */ + dijdik_inv=1.0/(d_ij*d_ik); + v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j + v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2); + v3_add(&dcosdrj,&dcosdrj,&tmp); + v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k + v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2); + v3_add(&dcosdrk,&dcosdrk,&tmp); + + /* f_c_ik * dg, df_c_ik * g */ + fcdg=f_c_ik*dg; + dfcg=df_c_ik*g; + + /* derivative wrt j */ + v3_scale(&force,&dcosdrj,fcdg*pre_dzeta); + + /* force contribution */ + 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); + printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); + printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI)); + printf(" d ij ik = %f %f\n",d_ij,d_ik); + } +} +#endif + + /* virial */ + virial_calc(ai,&force,&dist_ij); + + /* force contribution to atom i */ + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); + + /* derivative wrt k */ + 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); + + /* force contribution */ + v3_add(&(ak->f),&(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); + printf(" total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z); + printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI)); + printf(" d ij ik = %f %f\n",d_ij,d_ik); + } +} +#endif + + /* virial */ + virial_calc(ai,&force,&dist_ik); + + /* force contribution to atom i */ + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); + + /* increase k counter */ + exchange->kcount++; + + return 0; + +} + +int albe_orig_mult_check_2b_bond(t_moldyn *moldyn, + t_atom *itom,t_atom *jtom,u8 bc) { + + t_albe_orig_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; +} diff --git a/potentials/albe_orig.h b/potentials/albe_orig.h new file mode 100644 index 0000000..a52c7b2 --- /dev/null +++ b/potentials/albe_orig.h @@ -0,0 +1,100 @@ +/* + * albe_orig.h - albe potential header file + * + * author: Frank Zirkelbach + * + */ + +#ifndef ALBE_ORIG_H +#define ALBE_ORIG_H + +/* albe constants */ +#include "albe.h" + +#define ALBE_ORIG_MAXN 16*27 + +/* albe exchange type */ +typedef struct s_albe_orig_exchange { + + t_3dvec dist_ij; + double d_ij2; + double d_ij; + + t_3dvec dist_ik[ALBE_ORIG_MAXN]; + double d_ik2[ALBE_ORIG_MAXN]; + double d_ik[ALBE_ORIG_MAXN]; + + double f_c_ik[ALBE_ORIG_MAXN]; + double df_c_ik[ALBE_ORIG_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 zeta_ij; + double pre_dzeta; + + int kcount; +} t_albe_orig_exchange; + +/* albe mult (2!) potential parameters */ +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 */ + double Smixed; /* mixed S radius */ + double S2mixed; /* mixed S radius squared */ + double Rmixed; /* mixed R radius */ + double A[2]; /* factor of albe attractive part */ + double B[2]; /* factor of albe repulsive part */ + double r0[2]; /* r_0 */ + double Amixed; /* mixed A factor */ + double Bmixed; /* mixed B factor */ + double r0_mixed; /* mixed r_0 */ + double lambda[2]; /* albe lambda */ + double lambda_m; /* mixed lambda */ + double mu[2]; /* albe mu */ + double mu_m; /* mixed mu */ + + 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_orig_exchange exchange; /* exchange between 2bp and 3bp calc */ +} t_albe_orig_mult_params; + +/* function prototypes */ +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 */ +// -> albe.h + +#endif