+ }
+
+ /* check for hooks */
+ if(schedule->hook)
+ schedule->hook(moldyn,schedule->hook_params);
+
+ /* get a new info line */
+ printf("\n");
+
+ }
+
+ return 0;
+}
+
+/* velocity verlet */
+
+int velocity_verlet(t_moldyn *moldyn) {
+
+ int i,count;
+ double tau,tau_square;
+ t_3dvec delta;
+ t_atom *atom;
+
+ atom=moldyn->atom;
+ count=moldyn->count;
+ tau=moldyn->tau;
+ tau_square=moldyn->tau_square;
+
+ for(i=0;i<count;i++) {
+ /* new positions */
+ v3_scale(&delta,&(atom[i].v),tau);
+ v3_add(&(atom[i].r),&(atom[i].r),&delta);
+ v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
+ v3_add(&(atom[i].r),&(atom[i].r),&delta);
+ check_per_bound(moldyn,&(atom[i].r));
+
+ /* velocities */
+ v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
+ v3_add(&(atom[i].v),&(atom[i].v),&delta);
+ }
+
+ /* neighbour list update */
+ link_cell_update(moldyn);
+
+ /* forces depending on chosen potential */
+ potential_force_calc(moldyn);
+
+ for(i=0;i<count;i++) {
+ /* again velocities */
+ v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
+ v3_add(&(atom[i].v),&(atom[i].v),&delta);
+ }
+
+ return 0;
+}
+
+
+/*
+ *
+ * potentials & corresponding forces
+ *
+ */
+
+/* generic potential and force calculation */
+
+int potential_force_calc(t_moldyn *moldyn) {
+printf("start pot force calc\n");
+
+ int i,j,k,count;
+ t_atom *itom,*jtom,*ktom;
+ t_linkcell *lc;
+ t_list neighbour_i[27];
+ t_list neighbour_i2[27];
+ //t_list neighbour_j[27];
+ t_list *this,*that;
+ u8 bc_ij,bc_ik;
+ int dnlc;
+
+ count=moldyn->count;
+ itom=moldyn->atom;
+ lc=&(moldyn->lc);
+
+ /* reset energy */
+ moldyn->energy=0.0;
+
+ /* get energy and force of every atom */
+ for(i=0;i<count;i++) {
+printf("atom %d: %f\n",i,itom[i].r.x);
+printf("atom %d: %f\n",i,itom[i].v.x);
+printf("atom %d: %f\n",i,itom[i].f.x);
+
+ /* reset force */
+ v3_zero(&(itom[i].f));
+
+ /* single particle potential/force */
+ if(itom[i].attr&ATOM_ATTR_1BP)
+ moldyn->func1b(moldyn,&(itom[i]));
+
+ /* 2 body pair potential/force */
+ if(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
+
+ link_cell_neighbour_index(moldyn,
+ (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+ (itom[i].r.y+moldyn->dim.y/2)/lc->y,
+ (itom[i].r.z+moldyn->dim.z/2)/lc->z,
+ neighbour_i);
+
+ dnlc=lc->dnlc;
+
+ for(j=0;j<27;j++) {
+
+ this=&(neighbour_i[j]);
+ list_reset(this);
+
+ if(this->start==NULL)
+ continue;
+
+ bc_ij=(j<dnlc)?0:1;
+
+ do {
+ jtom=this->current->data;
+
+ 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);
+
+ /* 3 body potential/force */
+
+ if(!(itom[i].attr&ATOM_ATTR_3BP)||
+ !(jtom->attr&ATOM_ATTR_3BP))
+ continue;
+
+ /* copy the neighbour lists */
+ memcpy(neighbour_i2,neighbour_i,
+ 27*sizeof(t_list));
+
+ /* get neighbours of i */
+ for(k=0;k<27;k++) {
+
+ that=&(neighbour_i2[k]);
+ list_reset(that);
+
+ if(that->start==NULL)
+ continue;
+
+ bc_ik=(k<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_ik|bc_ij);
+
+ } while(list_next(that)!=\
+ L_NO_NEXT_ELEMENT);
+
+ }
+
+ } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+
+ /* 2bp post function */
+ if(moldyn->func2b_post) {
+ moldyn->func2b_post(moldyn,
+ &(itom[i]),
+ jtom,bc_ij);
+ }
+
+ }
+ }
+ }
+
+printf("end pot force calc\n");
+ return 0;
+}
+
+/*
+ * periodic boundayr checking
+ */
+
+int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
+
+ double x,y,z;
+ t_3dvec *dim;
+
+ dim=&(moldyn->dim);
+
+ x=0.5*dim->x;
+ y=0.5*dim->y;
+ z=0.5*dim->z;
+
+ if(moldyn->status&MOLDYN_STAT_PBX) {
+ if(a->x>=x) a->x-=dim->x;
+ else if(-a->x>x) a->x+=dim->x;
+ }
+ if(moldyn->status&MOLDYN_STAT_PBY) {
+ if(a->y>=y) a->y-=dim->y;
+ else if(-a->y>y) a->y+=dim->y;
+ }
+ if(moldyn->status&MOLDYN_STAT_PBZ) {
+ if(a->z>=z) a->z-=dim->z;
+ else if(-a->z>z) a->z+=dim->z;
+ }
+
+ return 0;
+}
+
+
+/*
+ * example potentials
+ */
+
+/* harmonic oscillator potential and force */
+
+int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+ t_ho_params *params;
+ t_3dvec force,distance;
+ double d;
+ double sc,equi_dist;
+
+ params=moldyn->pot2b_params;
+ sc=params->spring_constant;
+ equi_dist=params->equilibrium_distance;
+
+ v3_sub(&distance,&(aj->r),&(ai->r));
+
+ if(bc) check_per_bound(moldyn,&distance);
+ d=v3_norm(&distance);
+ if(d<=moldyn->cutoff) {
+ /* energy is 1/2 (d-d0)^2, but we will add this twice ... */
+ moldyn->energy+=(0.25*sc*(d-equi_dist)*(d-equi_dist));
+ /* f = -grad E; grad r_ij = -1 1/r_ij distance */
+ v3_scale(&force,&distance,sc*(1.0-(equi_dist/d)));
+ v3_add(&(ai->f),&(ai->f),&force);
+ }
+
+ return 0;
+}
+
+/* lennard jones potential & force for one sort of atoms */
+
+int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+ t_lj_params *params;
+ t_3dvec force,distance;
+ double d,h1,h2;
+ double eps,sig6,sig12;
+
+ params=moldyn->pot2b_params;
+ eps=params->epsilon4;
+ sig6=params->sigma6;
+ sig12=params->sigma12;
+
+ v3_sub(&distance,&(aj->r),&(ai->r));
+ if(bc) check_per_bound(moldyn,&distance);
+ d=v3_absolute_square(&distance); /* 1/r^2 */
+ if(d<=moldyn->cutoff_square) {
+ d=1.0/d; /* 1/r^2 */
+ h2=d*d; /* 1/r^4 */
+ h2*=d; /* 1/r^6 */
+ h1=h2*h2; /* 1/r^12 */
+ /* energy is eps*..., but we will add this twice ... */
+ moldyn->energy+=0.5*eps*(sig12*h1-sig6*h2);
+ h2*=d; /* 1/r^8 */
+ h1*=d; /* 1/r^14 */
+ h2*=6*sig6;
+ h1*=12*sig12;
+ d=+h1-h2;
+ d*=eps;
+ v3_scale(&force,&distance,-1.0*d); /* f = - grad E */
+ v3_add(&(ai->f),&(ai->f),&force);
+ }
+
+ return 0;
+}
+
+/*
+ * tersoff potential & force for 2 sorts of atoms
+ */
+
+/* create mixed terms from parameters and set them */
+int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
+
+ printf("[moldyn] tersoff parameter completion\n");
+ p->Smixed=sqrt(p->S[0]*p->S[1]);
+ 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]);
+ 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 (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(" 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;
+}
+
+/* tersoff 1 body part */
+int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
+
+ int num;
+ t_tersoff_mult_params *params;
+ t_tersoff_exchange *exchange;
+
+ num=ai->bnum;
+ params=moldyn->pot1b_params;
+ exchange=&(params->exchange);
+
+ /*
+ * simple: point constant parameters only depending on atom i to
+ * their right values
+ */
+
+ exchange->beta_i=&(params->beta[num]);
+ exchange->n_i=&(params->n[num]);
+ exchange->c_i=&(params->c[num]);
+ exchange->d_i=&(params->d[num]);
+ exchange->h_i=&(params->h[num]);
+
+ exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
+ exchange->ci2=params->c[num]*params->c[num];
+ exchange->di2=params->d[num]*params->d[num];
+ exchange->ci2di2=exchange->ci2/exchange->di2;
+
+ return 0;
+}
+
+/* tersoff 2 body part */
+int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+ 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 f_r,df_r;
+ double f_c,df_c;
+ int num;
+ double s_r;
+ double arg;
+
+ params=moldyn->pot2b_params;
+ num=aj->bnum;
+ exchange=&(params->exchange);
+
+ /* clear 3bp and 2bp post run */
+ exchange->run3bp=0;
+ exchange->run2bp_post=0;
+
+ /* reset S > r > R mark */
+ exchange->d_ij_between_rs=0;
+
+ /*
+ * calc of 2bp contribution of V_ij and dV_ij/ji
+ *
+ * for Vij and dV_ij we need:
+ * - f_c_ij, df_c_ij
+ * - f_r_ij, df_r_ij
+ *
+ * for dV_ji we need:
+ * - f_c_ji = f_c_ij, df_c_ji = df_c_ij
+ * - f_r_ji = f_r_ij; df_r_ji = df_r_ij
+ *
+ */
+
+ /* 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(num==ai->bnum) {
+ S=params->S[num];
+ R=params->R[num];
+ A=params->A[num];
+ B=params->B[num];
+ lambda=params->lambda[num];
+ mu=params->mu[num];
+ exchange->chi=1.0;
+ }
+ else {
+ S=params->Smixed;
+ R=params->Rmixed;
+ A=params->Amixed;
+ B=params->Bmixed;
+ lambda=params->lambda_m;
+ mu=params->mu_m;
+ params->exchange.chi=params->chi;
+ }
+
+ /* if d_ij > S => no force & potential energy contribution */
+ if(d_ij>S)
+ return 0;
+
+ /* more constants */
+ exchange->beta_j=&(params->beta[num]);
+ exchange->n_j=&(params->n[num]);
+ exchange->c_j=&(params->c[num]);
+ exchange->d_j=&(params->d[num]);
+ exchange->h_j=&(params->h[num]);
+ if(num==ai->bnum) {
+ exchange->betajnj=exchange->betaini;
+ exchange->cj2=exchange->ci2;
+ exchange->dj2=exchange->di2;
+ exchange->cj2dj2=exchange->ci2di2;
+ }
+ else {
+ exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
+ exchange->cj2=params->c[num]*params->c[num];
+ exchange->dj2=params->d[num]*params->d[num];
+ exchange->cj2dj2=exchange->cj2/exchange->dj2;
+ }
+
+ /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */
+ f_r=A*exp(-lambda*d_ij);
+ df_r=lambda*f_r/d_ij;
+
+ /* f_a, df_a calc (again, same for ij and ji) | save for later use! */
+ exchange->f_a=-B*exp(-mu*d_ij);
+ exchange->df_a=-mu*exchange->f_a/d_ij;
+
+ /* f_c, df_c calc (again, same for ij and ji) */
+ if(d_ij<R) {
+ /* f_c = 1, df_c = 0 */
+ f_c=1.0;
+ df_c=0.0;
+ /* two body contribution (ij, ji) */
+ v3_scale(&force,&dist_ij,-df_r);
+ }
+ 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));
+ /* two body contribution (ij, ji) */
+ v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
+ /* tell 3bp that S > r > R */
+ exchange->d_ij_between_rs=1;
+ }
+
+ /* add forces of 2bp (ij, ji) contribution
+ * dVij = dVji and we sum up both: no 1/2) */
+ v3_add(&(ai->f),&(ai->f),&force);
+
+ /* energy 2bp contribution (ij, ji) 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;
+ exchange->df_c=df_c;
+
+ /* enable the run of 3bp function and 2bp post processing */
+ exchange->run3bp=1;
+ exchange->run2bp_post=1;
+
+ /* reset 3bp sums */
+ exchange->zeta_ij=0.0;
+ exchange->zeta_ji=0.0;
+ v3_zero(&(exchange->dzeta_ij));
+ v3_zero(&(exchange->dzeta_ji));
+
+ return 0;
+}
+
+/* tersoff 2 body post part */