return 0;
}
-int set_temperature(t_moldyn *moldyn,double t) {
+int set_temperature(t_moldyn *moldyn,double t_ref) {
- moldyn->t=t;
+ moldyn->t_ref=t_ref;
+
+ return 0;
+}
+
+int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
+
+ moldyn->pt_scale=(ptype|ttype);
+ moldyn->t_tc=ttc;
+ moldyn->p_tc=ptc;
return 0;
}
return 0;
}
-int thermal_init(t_moldyn *moldyn) {
+int thermal_init(t_moldyn *moldyn,u8 equi_init) {
/*
* - gaussian distribution of velocities
/* gaussian distribution of velocities */
v3_zero(&p_total);
for(i=0;i<moldyn->count;i++) {
- sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
+ sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
/* x direction */
v=sigma*rand_get_gauss(random);
atom[i].v.x=v;
}
/* velocity scaling */
- scale_velocity(moldyn,VSCALE_INIT_EQUI);
+ scale_velocity(moldyn,equi_init);
return 0;
}
-int scale_velocity(t_moldyn *moldyn,u8 type) {
+int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
int i;
double e,scale;
* - velocity scaling (E = 3/2 N k T), E: kinetic energy
*/
+ /* get kinetic energy / temperature & count involved atoms */
e=0.0;
count=0;
for(i=0;i<moldyn->count;i++) {
- if(atom[i].attr&ATOM_ATTR_HB) {
+ if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
count+=1;
}
}
-
- /* temporary hack for e,t = 0 */
+ if(count!=0) moldyn->t=(2.0*e)/(3.0*count*K_BOLTZMANN);
+ else return 0; /* no atoms involved in scaling! */
+
+ /* (temporary) hack for e,t = 0 */
if(e==0.0) {
- if(moldyn->t!=0.0)
- thermal_init(moldyn);
+ moldyn->t=0.0;
+ if(moldyn->t_ref!=0.0)
+ thermal_init(moldyn,equi_init);
else
- return 0;
+ return 0; /* no scaling needed */
}
- /* direct scaling */
- scale=(1.5*count*K_BOLTZMANN*moldyn->t)/e;
- if(type&VSCALE_INIT_EQUI) scale*=2.0; /* equipartition theorem */
+
+ /* get scaling factor */
+ scale=moldyn->t_ref/moldyn->t;
+ if(equi_init&TRUE)
+ scale*=2.0;
+ else
+ if(moldyn->pt_scale&T_SCALE_BERENDSEN)
+ scale=1.0+moldyn->tau*(scale-1.0)/moldyn->t_tc;
scale=sqrt(scale);
+
+ /* velocity scaling */
for(i=0;i<moldyn->count;i++)
- if(atom[i].attr&ATOM_ATTR_HB)
+ if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
v3_scale(&(atom[i].v),&(atom[i].v),scale);
return 0;
/* integration step */
moldyn->integrate(moldyn);
+ /* p/t scaling */
+ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
+ scale_velocity(moldyn,FALSE);
+
/* increase absolute time */
moldyn->time+=moldyn->tau;
int i,j,k,count;
t_atom *itom,*jtom,*ktom;
t_linkcell *lc;
- t_list neighbour_i[27],neighbour_j[27];
+ t_list neighbour_i[27];
+ t_list neighbour_i2[27];
+ //t_list neighbour_j[27];
t_list *this,*that;
u8 bc_ij,bc_ijk;
int countn,dnlc;
!(jtom->attr&ATOM_ATTR_3BP))
continue;
- link_cell_neighbour_index(moldyn,
- (jtom->r.x+moldyn->dim.x/2)/lc->x,
- (jtom->r.y+moldyn->dim.y/2)/lc->y,
- (jtom->r.z+moldyn->dim.z/2)/lc->z,
- neighbour_j);
-
- /* neighbours of j */
- for(k=0;k<lc->countn;k++) {
-
- that=&(neighbour_j[k]);
- list_reset(that);
-
- if(that->start==NULL)
- continue;
-
- bc_ijk=(k<lc->dnlc)?0:1;
-
- do {
-
- ktom=that->current->data;
-
- if(!(ktom->attr&ATOM_ATTR_3BP))
- continue;
-
- if(ktom==jtom)
- continue;
-
- if(ktom==&(itom[i]))
- continue;
-
- moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
-
- } while(list_next(that)!=\
- L_NO_NEXT_ELEMENT);
-
- }
-
- /* neighbours of i */
+ /*
+ * according to mr. nordlund, we dont need to take the
+ * sum over all atoms now, as 'this is centered' around
+ * atom i ...
+ * i am not quite sure though! there is a not vanishing
+ * part even if f_c_ik is zero ...
+ * this analytical potentials suck!
+ * switching from mc to md to dft soon!
+ */
+
+ // link_cell_neighbour_index(moldyn,
+ // (jtom->r.x+moldyn->dim.x/2)/lc->x,
+ // (jtom->r.y+moldyn->dim.y/2)/lc->y,
+ // (jtom->r.z+moldyn->dim.z/2)/lc->z,
+ // neighbour_j);
+
+// /* neighbours of j */
+// for(k=0;k<lc->countn;k++) {
+//
+// that=&(neighbour_j[k]);
+// list_reset(that);
+//
+// if(that->start==NULL)
+// continue;
+//
+// bc_ijk=(k<lc->dnlc)?0:1;
+//
+// do {
+//
+// ktom=that->current->data;
+//
+// if(!(ktom->attr&ATOM_ATTR_3BP))
+// continue;
+//
+// if(ktom==jtom)
+// continue;
+//
+// if(ktom==&(itom[i]))
+// continue;
+//
+// moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+//
+/* } while(list_next(that)!=\ */
+// L_NO_NEXT_ELEMENT);
+//
+// }
+
+ /* copy the neighbour lists */
+ memcpy(neighbour_i2,neighbour_i,
+ 27*sizeof(t_list));
+
+ /* get neighbours of i */
for(k=0;k<countn;k++) {
- that=&(neighbour_i[k]);
+ that=&(neighbour_i2[k]);
list_reset(that);
if(that->start==NULL)
if(ktom==&(itom[i]))
continue;
+printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
} while(list_next(that)!=\
L_NO_NEXT_ELEMENT);
}
} while(list_next(this)!=L_NO_NEXT_ELEMENT);
+
+ /* 2bp post function */
+ if(moldyn->func2b_post)
+ mlodyn->func2b_post(moldyn,
+ &(itom[i]),
+ jtom,bc_ij);
+
}
}
}
p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
p->mu_m=0.5*(p->mu[0]+p->mu[1]);
+ printf("[moldyn] tersoff mult parameter info:\n");
+ printf(" S (m) | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed);
+ printf(" R (m) | %.12f | %.12f | %.12f\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(" beta | %.10f | %.10f\n",p->beta[0],p->beta[1]);
+ printf(" n | %f | %f\n",p->n[0],p->n[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(" chi | %f \n",p->chi);
+
return 0;
}
if(bc) check_per_bound(moldyn,&dist_ij);
- /* save for use in 3bp */ /* REALLY ?!?!?! */
+ d_ij=v3_norm(&dist_ij);
+
+ /* save for use in 3bp */
+ exchange->d_ij=d_ij;
exchange->dist_ij=dist_ij;
+ exchange->d_ij2=d_ij*d_ij;
/* constants */
if(num==aj->bnum) {
S=params->S[num];
R=params->R[num];
A=params->A[num];
+ B=params->B[num];
lambda=params->lambda[num];
- /* more constants depending of atoms i and j, needed in 3bp */
- params->exchange.B=&(params->B[num]);
- params->exchange.mu=&(params->mu[num]);
mu=params->mu[num];
params->exchange.chi=1.0;
}
S=params->Smixed;
R=params->Rmixed;
A=params->Amixed;
+ B=params->Bmixed;
lambda=params->lambda_m;
- /* more constants depending of atoms i and j, needed in 3bp */
- params->exchange.B=&(params->Bmixed);
- params->exchange.mu=&(params->mu_m);
mu=params->mu_m;
params->exchange.chi=params->chi;
}
- d_ij=v3_norm(&dist_ij);
-
- /* save for use in 3bp */
- exchange->d_ij=d_ij;
-
if(d_ij>S)
return 0;
/* add forces */
v3_add(&(ai->f),&(ai->f),&force);
- /* energy is 0.5 f_r f_c, but we will sum it up twice ... */
- moldyn->energy+=(0.25*f_r*f_c);
+ /* energy is 0.5 f_r f_c ... */
+ moldyn->energy+=(0.5*f_r*f_c);
/* save for use in 3bp */
exchange->f_c=f_c;
/* enable the run of 3bp function */
exchange->run3bp=1;
+ /* reset 3bp sums */
+ exchange->3bp_sum1=0.0;
+ exchange->3bp_sum2=0.0;
+
+ return 0;
+}
+
+/* tersoff 2 body post part */
+
+int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+
+ /* here we have to allow for the 3bp sums */
+
+ t_tersoff_mult_params *params;
+ t_tersoff_exchange *exchange;
+
+ t_3dvec force,temp,*db_ij;
+ double db_ij_scale1,db_ij_scale2;
+ double b_ij;
+ double f_c,df_c,f_a,df_a;
+
+ params=moldyn->pot2b_params;
+ exchange=&(moldyn->exchange);
+
+ db_ij=&(exchange->db_ij);
+ f_c=exchange->f_c;
+ df_c=exchange->df_c;
+ f_a=exchange->f_a;
+ df_a=exchange->df_a;
+
+ db_ij_scale1=(1+betan*3bp_sum1);
+ db_ij_scale2=(n*betan*3bp_sum2);
+ help=pow(db_ij_scale1,-1.0/(2*n)-1);
+ b_ij=chi*db_ij_scale1*help;
+ db_ij_scale1=-chi/(2*n)*help;
+
+ v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
+ v3_scale(db_ij,db_ij,f_a);
+
+ v3_scale(&temp,dist_ij,b_ij*df_a);
+
+ v3_add(&force,&temp,db_ij);
+ v3_scale(&force,&force,f_c);
+
+ v3_scale(&temp,&dist_ij,f_a*b_ij*df_c);
+
+ /* add energy of 3bp sum */
+ moldyn->energy+=(0.5*f_c*b_ij*f_a);
+ /* add force of 3bp calculation */
+ v3_add(&(ai->f),&temp,&force);
+
return 0;
}
d_ij=exchange->d_ij;
d_ij2=exchange->d_ij2;
+ dist_ij=exchange->dist_ij;
f_a=params->exchange.f_a;
df_a=params->exchange.df_a;
+
+ f_c=exchange->f_c;
+ df_c=exchange->df_c;
/* d_ij is <= S, as we didn't return so far! */
*
*/
-
v3_sub(&dist_ik,&(ai->r),&(ak->r));
if(bc) check_per_bound(moldyn,&dist_ik);
d_ik=v3_norm(&dist_ik);
}
/* calc of f_c_ik */
- if(d_ik>S)
- return 0;
-
- if(d_ik<R) {
- /* f_c_ik = 1, df_c_ik = 0 */
+ if(d_ik>S) {
+ f_c_ik=0.0;
+ df_c_ik=0.0;
+ }
+ else if(d_ik<R) {
f_c_ik=1.0;
df_c_ik=0.0;
}
c=*(exchange->c);
d=*(exchange->d);
h=*(exchange->h);
+ chi=exchange->chi;
c2=exchange->c2;
d2=exchange->d2;
c2d2=exchange->c2d2;
numer=d_ij2+d_ik*d_ik-d_jk*d_jk;
denom=2*d_ij*d_ik;
cos_theta=numer/denom;
+ //cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
sin_theta=sqrt(1.0-(cos_theta*cos_theta));
theta=acos(cos_theta);
d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom);
h_cos=(h-cos_theta);
h_cos2=h_cos*h_cos;
- d2_h_cos2=d2-h_cos2;
+ d2_h_cos2=d2+h_cos2;
/* some usefull expressions */
- frac1=c2/(d2-h_cos2);
+ frac1=c2/(d2_h_cos2);
bracket1=1+c2d2-frac1;
- bracket2=f_c_ik*bracket1;
- bracket2_n_1=pow(bracket2,n-1.0);
- bracket2_n=bracket2_n_1*bracket2;
- bracket3=1+betan*bracket2_n;
+ if(f_c_ik==0.0) {
+ bracket2=0.0;
+ bracket2_n_1=0.0;
+ bracket2_n=0.0;
+ bracket3=1.0;
+ printf("Foo -> 0: ");
+ }
+ else {
+ bracket2=f_c_ik*bracket1;
+ bracket2_n_1=pow(bracket2,n-1.0);
+ bracket2_n=bracket2_n_1*bracket2;
+ bracket3=1.0+betan*bracket2_n;
+ printf("Foo -> 1: ");
+ }
bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0);
bracket3_pow=bracket3_pow_1*bracket3;
+printf("%.15f %.15f %.15f\n",bracket2_n_1,bracket2_n);
/* now go on with calc of b_ij and derivation of b_ij */
b_ij=chi*bracket3_pow;
/* add forces */
v3_add(&(ai->f),&(ai->f),&force);
- /* energy is 0.5 f_r f_c, but we will sum it up twice ... */
- moldyn->energy+=(0.25*f_a*b_ij*f_c);
+ /* energy is 0.5 f_r f_c */
+ moldyn->energy+=(0.5*f_a*b_ij*f_c);
return 0;
}