#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
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))
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,
#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
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:
/* 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;
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=(j<dnlc)?0:1;
if(jtom==&(itom[i]))
continue;
+ /* reset 3bp run */
+ moldyn->run3bp=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=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+ q=0;
+
+ while(neighbour_i[j][q]!=0) {
+
+ ktom=&(atom[neighbour_i[k][q]]);
+ q++;
+#else
+ that=&(neighbour_i2[k]);
+ list_reset_f(that);
+
+ if(that->start==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
}
}
}
- /* 3 body potential/force */
+ /* continued 3 body potential/force */
if(!(itom[i].attr&ATOM_ATTR_3BP))
continue;
/* 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++) {
if(!(ktom->attr&ATOM_ATTR_3BP))
continue;
- if(ktom==jtom)
- continue;
+ //if(ktom==jtom)
+ // continue;
if(ktom==&(itom[i]))
continue;
}
- 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++) {
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
}
}
- /* 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
}
#endif
}
-
+
#ifdef DEBUG
//printf("\n\n");
#endif
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;
#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
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:
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->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:
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;i<ALBE_MAXN;i++)
+ memset(exchange->dzeta[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];
}
/* 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(d<R) {
+ exchange->f_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);
+ /* k<j & check whether to run k */
+ j=exchange->jcnt;
+ 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
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];
lambda=params->lambda[brand];
}
else {
- S=params->Smixed;
- R=params->Rmixed;
B=params->Bmixed;
A=params->Amixed;
r0=params->r0_mixed;
lambda=params->lambda_m;
}
- d_ij=exchange->d_ij;
-
- /* f_c, df_c */
- if(d_ij<R) {
- f_c=1.0;
- df_c=0.0;
- }
- else {
- s_r=S-R;
- arg=M_PI*(d_ij-R)/s_r;
- f_c=0.5+0.5*cos(arg);
- df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
- }
-
/* f_a, df_a */
- f_a=-B*exp(-mu*(d_ij-r0));
- df_a=mu*f_a/d_ij;
+ f_a=-B*exp(-mu*(d-r0));
+ df_a=mu*f_a/d;
/* f_r, df_r */
- f_r=A*exp(-lambda*(d_ij-r0));
- df_r=lambda*f_r/d_ij;
+ f_r=A*exp(-lambda*(d-r0));
+ df_r=lambda*f_r/d;
/* b, db */
- if(exchange->zeta_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);
}
#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;
}
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])) {
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) {
#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 */
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;
--- /dev/null
+/*
+ * albe_orig.c - albe potential
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <math.h>
+
+#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_ij<R) {
+ f_c=1.0;
+ df_c=0.0;
+ }
+ else {
+ s_r=S-R;
+ arg=M_PI*(d_ij-R)/s_r;
+ f_c=0.5+0.5*cos(arg);
+ df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
+ }
+
+ /* f_a, df_a */
+ f_a=-B*exp(-mu*(d_ij-r0));
+ df_a=mu*f_a/d_ij;
+
+ /* f_r, df_r */
+ f_r=A*exp(-lambda*(d_ij-r0));
+ df_r=lambda*f_r/d_ij;
+
+ /* b, db */
+ if(exchange->zeta_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->time<DEND) {
+ if((ai==&(moldyn->atom[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->time<DEND) {
+ if(aj==&(moldyn->atom[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->time<DEND) {
+ if(ak==&(moldyn->atom[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;
+}
--- /dev/null
+/*
+ * albe_orig.h - albe potential header file
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#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
projects / physik / posic.git / commitdiff