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;
+ }
+
+ /* just a guess so far ... */
+ v=virial.xx+virial.yy+virial.zz;
-printf("%f\n",virial);
+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)
t_moldyn_schedule *schedule;
schedule=&(moldyn->schedule);
- count=++(schedule->content_count);
+ count=++(schedule->total_sched);
- ptr=realloc(moldyn->schedule.runs,count*sizeof(int));
+ ptr=realloc(schedule->runs,count*sizeof(int));
if(!ptr) {
perror("[moldyn] realloc (runs)");
return -1;
}
- moldyn->schedule.runs=ptr;
- moldyn->schedule.runs[count-1]=runs;
+ schedule->runs=ptr;
+ schedule->runs[count-1]=runs;
ptr=realloc(schedule->tau,count*sizeof(double));
if(!ptr) {
perror("[moldyn] realloc (tau)");
return -1;
}
- moldyn->schedule.tau=ptr;
- moldyn->schedule.tau[count-1]=tau;
+ schedule->tau=ptr;
+ schedule->tau[count-1]=tau;
return 0;
}
int moldyn_integrate(t_moldyn *moldyn) {
- int i,sched;
+ int i;
unsigned int e,m,s,v;
t_3dvec p;
- t_moldyn_schedule *schedule;
+ t_moldyn_schedule *sched;
t_atom *atom;
int fd;
char dir[128];
double ds;
- schedule=&(moldyn->schedule);
+ sched=&(moldyn->schedule);
atom=moldyn->atom;
/* initialize linked cell method */
moldyn->debug=0;
/* executing the schedule */
- for(sched=0;sched<moldyn->schedule.content_count;sched++) {
+ for(sched->count=0;sched->count<sched->total_sched;sched->count++) {
/* setting amount of runs and finite time step size */
- moldyn->tau=schedule->tau[sched];
+ moldyn->tau=sched->tau[sched->count];
moldyn->tau_square=moldyn->tau*moldyn->tau;
- moldyn->time_steps=schedule->runs[sched];
+ moldyn->time_steps=sched->runs[sched->count];
/* integration according to schedule */
visual_atoms(&(moldyn->vis),moldyn->time,
moldyn->atom,moldyn->count);
printf("\rsched: %d, steps: %d, debug: %d",
- sched,i,moldyn->debug);
+ sched->count,i,moldyn->debug);
fflush(stdout);
}
}
}
/* check for hooks */
- if(schedule->hook)
- schedule->hook(moldyn,schedule->hook_params);
+ if(sched->hook)
+ sched->hook(moldyn,sched->hook_params);
/* get a new info line */
printf("\n");
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;
}
}
+#ifdef DEBUG
+printf("\n\n");
+#endif
+#ifdef VDEBUG
+printf("\n\n");
+#endif
return 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->d_ij2=d_ij2;
+ 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);
* 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");
+ printf("%f | %f\n",force.x,ai->f.x);
+ printf("%f | %f\n",force.y,ai->f.y);
+ printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVij, dVji (2bp) contrib:\n");
+ printf("%f | %f\n",force.x*dist_ij.x,ai->virial.xx);
+ printf("%f | %f\n",force.y*dist_ij.y,ai->virial.yy);
+ printf("%f | %f\n",force.z*dist_ij.z,ai->virial.zz);
+}
+#endif
+
/* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
moldyn->energy+=(0.5*f_r*f_c);
/* 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");
+ printf("%f | %f\n",force.x,ai->f.x);
+ printf("%f | %f\n",force.y,ai->f.y);
+ printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVij (3bp) contrib:\n");
+ printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx);
+ printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy);
+ printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz);
+}
+#endif
+
/* add energy of 3bp sum */
moldyn->energy+=(0.5*f_c*b*f_a);
/* 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");
+ printf("%f | %f\n",force.x,ai->f.x);
+ printf("%f | %f\n",force.y,ai->f.y);
+ printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVji (3bp) contrib:\n");
+ printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx);
+ printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy);
+ printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz);
+}
+#endif
+
return 0;
}
t_3dvec dist_ij,dist_ik,dist_jk;
t_3dvec temp1,temp2;
t_3dvec *dzeta;
- double R,S,s_r;
+ double R,S,S2,s_r;
double B,mu;
- double d_ij,d_ik,d_jk;
+ double d_ij,d_ik,d_jk,d_ij2,d_ik2,d_jk2;
double rr,dd;
double f_c,df_c;
double f_c_ik,df_c_ik,arg;
/* dist_ij, d_ij - this is < S_ij ! */
dist_ij=exchange->dist_ij;
d_ij=exchange->d_ij;
+ d_ij2=exchange->d_ij2;
/* f_c_ij, df_c_ij (same for ji) */
f_c=exchange->f_c;
/* dist_ik, d_ik */
v3_sub(&dist_ik,&(ak->r),&(ai->r));
if(bc) check_per_bound(moldyn,&dist_ik);
- d_ik=v3_norm(&dist_ik);
+ d_ik2=v3_absolute_square(&dist_ik);
/* ik constants */
brand=ai->brand;
if(brand==ak->brand) {
R=params->R[brand];
S=params->S[brand];
+ S2=params->S2[brand];
}
else {
R=params->Rmixed;
S=params->Smixed;
+ S2=params->S2mixed;
}
/* zeta_ij/dzeta_ij contribution only for d_ik < S */
- if(d_ik<S) {
+ if(d_ik2<S2) {
+
+ /* now we need d_ik */
+ d_ik=sqrt(d_ik2);
/* get constants_i from exchange data */
n=*(exchange->n_i);
/* d_costheta */
tmp=1.0/dd;
- d_costheta1=cos_theta/(d_ij*d_ij)-tmp;
- d_costheta2=cos_theta/(d_ik*d_ik)-tmp;
+ d_costheta1=cos_theta/d_ij2-tmp;
+ d_costheta2=cos_theta/d_ik2-tmp;
/* some usefull values */
h_cos=(h-cos_theta);
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 */
/* dist_jk, d_jk */
v3_sub(&dist_jk,&(ak->r),&(aj->r));
if(bc) check_per_bound(moldyn,&dist_jk);
- d_jk=v3_norm(&dist_jk);
+ d_jk2=v3_absolute_square(&dist_jk);
/* jk constants */
brand=aj->brand;
if(brand==ak->brand) {
R=params->R[brand];
S=params->S[brand];
+ S2=params->S2[brand];
B=params->B[brand];
mu=params->mu[brand];
chi=1.0;
else {
R=params->Rmixed;
S=params->Smixed;
+ S2=params->S2mixed;
B=params->Bmixed;
mu=params->mu_m;
chi=params->chi;
}
/* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */
- if(d_jk<S) {
+ if(d_jk2<S2) {
+
+ /* now we need d_ik */
+ d_jk=sqrt(d_jk2);
/* constants_j from exchange data */
n=*(exchange->n_j);
/* d_costheta */
d_costheta1=1.0/dd;
- d_costheta2=cos_theta/(d_ij*d_ij);
+ d_costheta2=cos_theta/d_ij2;
/* some usefull values */
h_cos=(h-cos_theta);
/* g(cos_theta) */
g=1.0+c2d2-frac;
- /* d_costheta_ij and dg(cos_theta) - needed in any case! */
+ /* d_costheta_jik and dg(cos_theta) - needed in any case! */
v3_scale(&temp1,&dist_jk,d_costheta1);
v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */
- v3_add(&temp1,&temp1,&temp2);
+ //v3_add(&temp1,&temp1,&temp2);
+ v3_sub(&temp1,&temp1,&temp2); /* there is a minus! */
v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
/* store dg in temp2 and use it for dVjk later */
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");
+ printf("%f | %f\n",temp2.x,ai->f.x);
+ printf("%f | %f\n",temp2.y,ai->f.y);
+ printf("%f | %f\n",temp2.z,ai->f.z);
+}
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVjk (3bp) contrib:\n");
+ printf("%f | %f\n",temp2.x*dist_jk.x,ai->virial.xx);
+ printf("%f | %f\n",temp2.y*dist_jk.y,ai->virial.yy);
+ printf("%f | %f\n",temp2.z*dist_jk.z,ai->virial.zz);
+}
+#endif
+
}
return 0;
projects / physik / posic.git / blobdiff