t_atom *atom;
t_3dvec *dim,*vdim;
- double virial,scale;
+ double scale,v;
+ t_virial virial;
t_linkcell *lc;
int i;
vdim=&(moldyn->vis.dim);
lc=&(moldyn->lc);
- for(i=0;i<moldyn->count;i++)
- virial+=v3_norm(&(atom[i].virial));
+ memset(&virial,0,sizeof(t_virial));
+
+ for(i=0;i<moldyn->count;i++) {
+ virial.xx+=atom[i].virial.xx;
+ virial.yy+=atom[i].virial.yy;
+ virial.zz+=atom[i].virial.zz;
+ virial.xy+=atom[i].virial.xy;
+ virial.xz+=atom[i].virial.xz;
+ virial.yz+=atom[i].virial.yz;
+ }
-printf("%f\n",virial);
+ /* just a guess so far ... */
+ v=virial.xx+virial.yy+virial.zz;
+
+printf("%f\n",v);
/* get pressure from virial */
- moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial;
+ moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+ONE_THIRD*v;
moldyn->p/=moldyn->volume;
-printf("%f\n",moldyn->p/(ATM));
+printf("%f | %f\n",moldyn->p/(ATM),moldyn->p_ref/ATM);
/* scale factor */
if(moldyn->pt_scale&P_SCALE_BERENDSEN)
int i,j,k,count;
t_atom *itom,*jtom,*ktom;
+ t_virial *virial;
t_linkcell *lc;
t_list neighbour_i[27];
t_list neighbour_i2[27];
v3_zero(&(itom[i].f));
/* reset viral of atom i */
- v3_zero(&(itom[i].virial));
+ virial=&(itom[i].virial);
+ virial->xx=0.0;
+ virial->yy=0.0;
+ virial->zz=0.0;
+ virial->xy=0.0;
+ virial->xz=0.0;
+ virial->yz=0.0;
/* reset site energy */
itom[i].e=0.0;
int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
printf("[moldyn] tersoff parameter completion\n");
+ p->S2[0]=p->S[0]*p->S[0];
+ p->S2[1]=p->S[1]*p->S[1];
p->Smixed=sqrt(p->S[0]*p->S[1]);
+ p->S2mixed=p->Smixed*p->Smixed;
p->Rmixed=sqrt(p->R[0]*p->R[1]);
p->Amixed=sqrt(p->A[0]*p->A[1]);
p->Bmixed=sqrt(p->B[0]*p->B[1]);
t_tersoff_mult_params *params;
t_tersoff_exchange *exchange;
t_3dvec dist_ij,force;
- double d_ij;
- double A,B,R,S,lambda,mu;
+ double d_ij,d_ij2;
+ double A,B,R,S,S2,lambda,mu;
double f_r,df_r;
double f_c,df_c;
int brand;
*
*/
- /* dist_ij, d_ij */
- v3_sub(&dist_ij,&(aj->r),&(ai->r));
- if(bc) check_per_bound(moldyn,&dist_ij);
- d_ij=v3_norm(&dist_ij);
-
- /* save for use in 3bp */
- exchange->d_ij=d_ij;
- exchange->dist_ij=dist_ij;
-
/* constants */
if(brand==ai->brand) {
S=params->S[brand];
+ S2=params->S2[brand];
R=params->R[brand];
A=params->A[brand];
B=params->B[brand];
}
else {
S=params->Smixed;
+ S2=params->S2mixed;
R=params->Rmixed;
A=params->Amixed;
B=params->Bmixed;
params->exchange.chi=params->chi;
}
- /* if d_ij > S => no force & potential energy contribution */
- if(d_ij>S)
+ /* dist_ij, d_ij */
+ 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)
return 0;
+ /* now we will need the distance */
+ //d_ij=v3_norm(&dist_ij);
+ d_ij=sqrt(d_ij2);
+
+ /* save for use in 3bp */
+ exchange->d_ij=d_ij;
+ exchange->dist_ij=dist_ij;
+
/* more constants */
exchange->beta_j=&(params->beta[brand]);
exchange->n_j=&(params->n[brand]);
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)); /* MARK! */
df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
/* two body contribution (ij, ji) */
v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
/* add forces of 2bp (ij, ji) contribution
* dVij = dVji and we sum up both: no 1/2) */
v3_add(&(ai->f),&(ai->f),&force);
+
+ /* virial */
+ ai->virial.xx-=force.x*dist_ij.x;
+ ai->virial.yy-=force.y*dist_ij.y;
+ ai->virial.zz-=force.z*dist_ij.z;
+ ai->virial.xy-=force.x*dist_ij.y;
+ ai->virial.xz-=force.x*dist_ij.z;
+ ai->virial.yz-=force.y*dist_ij.z;
+
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVij, dVji (2bp) contrib:\n");
/* add force */
v3_add(&(ai->f),&(ai->f),&force);
+
+ /* virial */
+ ai->virial.xx-=force.x*dist_ij->x;
+ ai->virial.yy-=force.y*dist_ij->y;
+ ai->virial.zz-=force.z*dist_ij->z;
+ ai->virial.xy-=force.x*dist_ij->y;
+ ai->virial.xz-=force.x*dist_ij->z;
+ ai->virial.yz-=force.y*dist_ij->z;
+
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVij (3bp) contrib:\n");
/* add force */
v3_add(&(ai->f),&(ai->f),&force);
+
+ /* virial - plus sign, as dist_ij = - dist_ji - (really??) */
+ ai->virial.xx+=force.x*dist_ij->x;
+ ai->virial.yy+=force.y*dist_ij->y;
+ ai->virial.zz+=force.z*dist_ij->z;
+ ai->virial.xy+=force.x*dist_ij->y;
+ ai->virial.xz+=force.x*dist_ij->z;
+ ai->virial.yz+=force.y*dist_ij->z;
+
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVji (3bp) contrib:\n");
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)); /* MARK */
df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
/* zeta_ij */
v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
/* scaled with 0.5 ^ */
+
+ /* virial */
+ ai->virial.xx-=temp2.x*dist_jk.x;
+ ai->virial.yy-=temp2.y*dist_jk.y;
+ ai->virial.zz-=temp2.z*dist_jk.z;
+ ai->virial.xy-=temp2.x*dist_jk.y;
+ ai->virial.xz-=temp2.x*dist_jk.z;
+ ai->virial.yz-=temp2.y*dist_jk.z;
+
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
printf("dVjk (3bp) contrib:\n");