switch(type) {
case MOLDYN_POTENTIAL_TM:
- moldyn->func1b=tersoff_mult_1bp;
- moldyn->func3b_j1=tersoff_mult_3bp_j1;
- moldyn->func3b_k1=tersoff_mult_3bp_k1;
- moldyn->func3b_j2=tersoff_mult_3bp_j2;
- moldyn->func3b_k2=tersoff_mult_3bp_k2;
+ moldyn->func_i0=tersoff_mult_1bp;
+ moldyn->func_j1=tersoff_mult_3bp_j1;
+ moldyn->func_j1_k0=tersoff_mult_3bp_k1;
+ moldyn->func_j1c=tersoff_mult_3bp_j2;
+ moldyn->func_j1_k1=tersoff_mult_3bp_k2;
moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
break;
/*
break;
*/
case MOLDYN_POTENTIAL_AM:
- moldyn->func1b=albe_mult_i0;
- moldyn->func2b=albe_mult_i0_j0;
- moldyn->func3b_0=albe_mult_i0_j0_k0;
- moldyn->func3b_1=albe_mult_i0_j1;
- moldyn->func3b_j1=albe_mult_i0_j2;
- moldyn->func3b_k2=albe_mult_i0_j2_k0;
- moldyn->func3b_j2=albe_mult_i0_j3;
+ 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;
jtom=&(atom[neighbour_i[j][p]]);
p++;
-
- if(jtom==&(itom[i]))
- continue;
-
- if((jtom->attr&ATOM_ATTR_2BP)&
- (itom[i].attr&ATOM_ATTR_2BP)) {
- moldyn->func2b(moldyn,
- &(itom[i]),
- jtom,
- bc_ij);
- }
- }
#else
this=&(neighbour_i[j]);
list_reset_f(this);
do {
jtom=this->current->data;
+#endif
if(jtom==&(itom[i]))
continue;
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);
}
-// REWRITE ONCE WORKING!!!
/* 3 body potential/force */
/* in j loop, 3bp run can be skipped */
if(!(jtom->attr&ATOM_ATTR_3BP))
continue;
- if(moldyn->func3b_0==NULL)
+ 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;
if(!(ktom->attr&ATOM_ATTR_3BP))
continue;
- if(ktom==jtom)
- continue;
+ //if(ktom==jtom)
+ // continue;
if(ktom==&(itom[i]))
continue;
- moldyn->func3b_0(moldyn,
- &(itom[i]),
- jtom,
- ktom,
- bc_ik|bc_ij);
+ moldyn->func_j0_k0(moldyn,
+ &(itom[i]),
+ jtom,
+ ktom,
+ bc_ik|bc_ij);
#ifdef STATIC_LISTS
}
#else
}
-
- if(moldyn->func3b_1)
- moldyn->func3b_1(moldyn,
+ /* finish of first j loop after first k loop */
+ if(moldyn->func_j0e)
+ moldyn->func_j0e(moldyn,
&(itom[i]),
jtom,
bc_ij);
-// UNTIL HERE + change function names!
-
-
-
+#ifdef STATIC_LISTS
+ }
+#else
} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
#endif
}
}
- /* 3 body potential/force */
+ /* continued 3 body potential/force */
if(!(itom[i].attr&ATOM_ATTR_3BP))
continue;
- /* copy the neighbour lists */
-#ifdef STATIC_LISTS
- /* no copy needed for static lists */
-#else
- memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
-#endif
-
/* second loop over atoms j */
for(j=0;j<27;j++) {
/* 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;
- moldyn->func3b_k1(moldyn,
- &(itom[i]),
- jtom,
- ktom,
- bc_ik|bc_ij);
+ moldyn->func_j1_k0(moldyn,
+ &(itom[i]),
+ jtom,
+ ktom,
+ bc_ik|bc_ij);
#ifdef STATIC_LISTS
}
#else
}
- 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_0)(struct s_moldyn *moldyn,
- t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
- int (*func3b_1)(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;
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;
}
/* first i loop */
int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
- int i;
t_albe_mult_params *params;
t_albe_exchange *exchange;
+ int i;
+
params=moldyn->pot_params;
exchange=&(params->exchange);
/* zero exchange values */
- memset(exchange->dist,0,ALBE_MAXN*sizeof(t_3dvec));
- memset(exchange->d2,0,ALBE_MAXN*sizeof(double));
- memset(exchange->d,0,ALBE_MAXN*sizeof(double));
memset(exchange->zeta,0,ALBE_MAXN*sizeof(double));
for(i=0;i<ALBE_MAXN;i++)
- memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(double));
- memset(exchange->skip,0,ALBE_MAXN*sizeof(u8));
+ memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(t_3dvec));
exchange->jcnt=0;
exchange->j2cnt=0;
t_albe_mult_params *params;
t_albe_exchange *exchange;
+
double S2,S,R,d2,d,s_r,arg;
t_3dvec dist;
int j;
/* dist_ij, d_ij2 */
v3_sub(&dist,&(aj->r),&(ai->r));
- exchange->dist[j]=dist;
if(bc) check_per_bound(moldyn,&dist);
+ exchange->dist[j]=dist;
d2=v3_absolute_square(&dist);
exchange->d2[j]=d2;
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];
- /* set albe needs i,(j/k) depending c,d,h and gamma values */
+ /* 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;
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);
}
- exchange->c2_[j]=*exchange->c_[j]**exchange->c_[j];
- exchange->d2_[j]=*exchange->d_[j]**exchange->d_[j];
- exchange->c2d2_[j]=exchange->c2_[j]/exchange->d2_[j];
/* d_ij */
- d=sqrt(exchange->d2[j]);
+ d=sqrt(d2);
exchange->d[j]=d;
/* f_c, df_c */
t_albe_mult_params *params;
t_albe_exchange *exchange;
+
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,**dzeta;
+ t_3dvec dcosdrj,dcosdrk,tmp;
+ t_3dvec *dzjj,*dzkk,*dzjk,*dzkj;
params=moldyn->pot_params;
exchange=&(params->exchange);
+ if(aj==ak) {
+ exchange->kcnt+=1;
+ return 0;
+ }
+
/* k<j & check whether to run k */
j=exchange->jcnt;
k=exchange->kcnt;
exchange->kcnt+=1;
return 0;
}
-
+
/* distances */
distj=exchange->dist[j];
distk=exchange->dist[k];
/* cos theta */
cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
- /* g(cos(theta)) - in albe: ik-depending values! */
+ /* 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;
+ 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) {
}
else {
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;
+ 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;
}
- /* zeta */
+ /* zeta - for albe: ik depending g function */
+//if(ai->tag==0) {
+// printf("------> %.15f %.15f\n",dj,dk);
+// printf("------> %.15f %.15f\n",dj,dk);
+//}
+
exchange->zeta[j]+=(exchange->f_c[k]*gk);
exchange->zeta[k]+=(exchange->f_c[j]*gj);
v3_add(&dcosdrk,&dcosdrk,&tmp);
/* zeta derivatives */
- dzeta=exchange->dzeta;
+ dzjj=&(exchange->dzeta[j][j]);
+ dzkk=&(exchange->dzeta[k][k]);
+ dzjk=&(exchange->dzeta[j][k]);
+ dzkj=&(exchange->dzeta[k][j]);
v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
- v3_add(&dzeta[j][j],&dzeta[j][j],&tmp); // j j
+ v3_add(dzjj,dzjj,&tmp); // j j
v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
- v3_add(&dzeta[k][k],&dzeta[k][k],&tmp); // k k
- v3_scale(&tmp,&distk,exchange->df_c[k]*gk/dk);
- v3_add(&dzeta[j][k],&dzeta[j][k],&tmp);
+ v3_add(dzkk,dzkk,&tmp); // k k
+ v3_scale(&tmp,&distk,-exchange->df_c[k]*gk); // dk rik = - di rik
+ v3_add(dzjk,dzjk,&tmp);
v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
- v3_add(&dzeta[j][k],&dzeta[j][k],&tmp); // j k
- v3_scale(&tmp,&distj,exchange->df_c[j]*gj/dj);
- v3_add(&dzeta[k][j],&dzeta[k][j],&tmp);
+ v3_add(dzjk,dzjk,&tmp); // j k
+ v3_scale(&tmp,&distj,-exchange->df_c[j]*gj); // dj rij = - di rij
+ v3_add(dzkj,dzkj,&tmp);
v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
- v3_add(&dzeta[k][j],&dzeta[k][j],&tmp); // k j
+ v3_add(dzkj,dzkj,&tmp); // k j
/* increase k counter */
exchange->kcnt+=1;
v3_scale(&force,dist,scale);
v3_add(&(ai->f),&(ai->f),&force);
+#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);
+#ifdef NDEBUG
+if(aj->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (ji) %.15f\n",force.x,force.y,force.z,aj->tag,ai->tag,exchange->zeta[j]);
+printf(" t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
+}
+#endif
+
/* virial */
virial_calc(aj,&force,dist);
/* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
- exchange->pre_dzeta=0.5*f_a*exchange->f_c[j]*db;
+ exchange->pre_dzeta=0.5*f_a*f_c*db;
/* 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
+
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
+
/* virial */
virial_calc(ai,&force,dist);
params=moldyn->pot_params;
exchange=&(params->exchange);
+ if(aj==ak) {
+ exchange->kcnt+=1;
+ return 0;
+ }
+
/* k!=j & check whether to run k */
k=exchange->kcnt;
j=exchange->j2cnt;
/* force contribution (drk derivative) */
v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta);
v3_add(&(ak->f),&(ak->f),&force);
+
+#ifdef NDEBUG
+if(ak->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d %d (k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag);
+printf(" t: %.15f %.15f %.15f\n",ak->f.x,ak->f.y,ak->f.z);
+}
+#endif
+
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 %d -- %d(i contr k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag,k);
+printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
+printf(" ## %f\n",exchange->d[k]);
+}
+#endif
+
/* virial */
virial_calc(ai,&force,&(exchange->dist[k]));
projects / physik / posic.git / commitdiff