+ /* close file */
+ close(fd);
+
+ return 0;
+}
+
+int check_pressure(t_moldyn *moldyn,t_mdrun *mdrun) {
+
+ double delta;
+
+ if(!(mdrun->sattr&SATTR_PRELAX))
+ return TRUE;
+
+ delta=moldyn->p_avg-moldyn->p_ref;
+
+ if(delta<0)
+ delta=-delta;
+
+ if(delta>mdrun->dp)
+ return FALSE;
+
+ return TRUE;
+}
+
+int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) {
+
+ double delta;
+
+ if(!(mdrun->sattr&SATTR_TRELAX))
+ return TRUE;
+
+ delta=moldyn->t_avg-moldyn->t_ref;
+
+ if(delta<0)
+ delta=-delta;
+
+ if(delta>mdrun->dt)
+ return FALSE;
+
+ return TRUE;
+}
+
+int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
+
+ insert_atoms_params *iap;
+ t_stage *stage;
+ t_atom *atom;
+ t_3dvec r,v,dist;
+ double d,dmin;
+ int cnt,i;
+ double x,y,z;
+ double x0,y0,z0;
+ u8 cr_check,run;
+
+ stage=mdrun->stage->current->data;
+ iap=stage->params;
+
+ cr_check=FALSE;
+
+ v.x=0.0; v.y=0.0; v.z=0.0;
+
+ switch(iap->type) {
+ case INS_TOTAL:
+ x=moldyn->dim.x;
+ x0=0.0;
+ y=moldyn->dim.y;
+ y0=0.0;
+ z=moldyn->dim.z;
+ z0=0.0;
+ cr_check=TRUE;
+ break;
+ case INS_REGION:
+ x=iap->x1-iap->x0;
+ x0=iap->x0;
+ y=iap->y1-iap->y0;
+ y0=iap->y0;
+ z=iap->z1-iap->z0;
+ z0=iap->z0;
+ cr_check=TRUE;
+ break;
+ default:
+ printf("%s unknown insertion mode\n");
+ return -1;
+ }
+
+ cnt=0;
+ while(cnt<iap->ins_atoms) {
+ run=1;
+ while(run) {
+ r.x=(rand_get_double(&(moldyn->random))-0.5)*x+x0;
+ r.y=(rand_get_double(&(moldyn->random))-0.5)*y+y0;
+ r.z=(rand_get_double(&(moldyn->random))-0.5)*z+z0;
+ run=0;
+ if(cr_check==TRUE) {
+ dmin=1000; // for sure too high!
+ for(i=0;i<moldyn->count;i++) {
+ atom=&(moldyn->atom[i]);
+ v3_sub(&dist,&(atom->r),&r);
+ check_per_bound(moldyn,&dist);
+ d=v3_absolute_square(&dist);
+ if(d<(iap->cr*iap->cr)) {
+ run=1;
+ break;
+ }
+ if(d<dmin)
+ dmin=d;
+ }
+ }
+ }
+ add_atom(moldyn,iap->element,pse_mass[iap->element],
+ iap->brand,iap->attr,&r,&v);
+ printf("%s atom inserted: %f %f %f | d squared = %f\n",
+ ME,r.x,r.y,r.z,dmin);
+ cnt+=1;
+ }
+
+ return 0;
+}
+
+int chaatrib(t_moldyn *moldyn,t_mdrun *mdrun) {
+
+ t_stage *stage;
+ t_chaattr_params *cap;
+ t_atom *atom;
+ u8 assigne;
+
+ stage=mdrun->stage->current->data;
+ cap=stage->params;
+
+ for(i=0;i<moldyn->count;i++) {
+ atom=&(moldyn->atom[i]);
+ if(cap->type&CHAATTR_ELEMENT) {
+ if(cap->element!=atom->element)
+ continue;
+ }
+ if(cap->type&CHAATTR_REGION) {
+ if(cap->x0<atom->r.x)
+ continue;
+ if(cap->y0<atom->r.y)
+ continue;
+ if(cap->z0<atom->r.z)
+ continue;
+ if(cap->x1>atom->r.x)
+ continue;
+ if(cap->y1>atom->r.y)
+ continue;
+ if(cap->z1>atom->r.z)
+ continue;
+ }
+ atom->attr=cap->attr;
+ }
+
+ return 0;
+}
+
+int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
+
+ t_stage *stage;
+ t_chsattr_params *csp;
+
+ stage=mdrun->stage->current->data;
+ csp=stage->params;
+
+ switch(csp->type) {
+ case CHSATTR_PCTRL:
+ set_p_scale(moldyn,csp->ctrl_type,csp->tau);
+ break;
+ case CHSATTR_TCTRL:
+ set_t_scale(moldyn,csp->ctrl_type,csp->tau);
+ break;
+ case CHSATTR_PRELAX:
+ mdrun->sattr&=(^(SATTR_PRELAX));
+ if(csp->ctrl==TRUE)
+ mdrun->sattr|=SATTR_PRELAX;
+ break;
+ case CHSATTR_TRELAX:
+ mdrun->sattr&=(^(SATTR_TRELAX));
+ if(csp->ctrl==TRUE)
+ mdrun->sattr|=SATTR_TRELAX;
+ break;
+ case CHSATTR_AVGRST:
+ mdrun->sattr&=(^(SATTR_AVGRST));
+ if(csp->ctrl==TRUE)
+ mdrun->sattr|=SATTR_AVGRST;
+ break;
+ default:
+ printf("%s unknown system attribute change\n",ME);
+ return -1;
+ }
+
+ return 0;
+}
+
+int mdrun_hook(t_moldyn *moldyn,void *ptr) {
+
+ t_mdrun *mdrun;
+ t_stage *stage;
+ t_list *sl;
+ int steps;
+ double tau;
+ u8 change_stage;
+
+ t_insert_atoms_params *iap;
+ t_continue_params *cp;
+ t_anneal_params *ap;
+
+ mdrun=ptr;
+ sl=mdrun->stage;
+
+ change_stage=FALSE;
+
+ /* return immediately if there are no more stage */
+ if(sl->current==NULL)
+ return 0;
+
+ /* get stage description */
+ stage=mdrun->sl->current->data;
+
+ /* default steps and tau values */
+ steps=mdrun->relax_steps;
+ tau=mdrun->timestep;
+
+ /* check whether relaxation steps are necessary */
+ if(!((check_pressure==FALSE)|(check_temperature==FALSE))) {
+
+ /* stage specific stuff */
+ switch(stage->type) {
+ case STAGE_INSERT_ATOMS:
+ iap=stage->params;
+ if(iap->cnt_steps==iap->ins_steps) {
+ change_stage=TRUE;
+ break;
+ }
+ insert_atoms(moldyn,mdrun);
+ iap->cnt_steps+=1;
+ break;
+ case STAGE_CONTINUE:
+ if(stage->executed==TRUE) {
+ change_stage=TRUE;
+ break;
+ }
+ cp=stage->params;
+ steps=cp->runs;
+ break;
+ case STAGE_ANNEAL:
+ ap=stage->params;
+ if(ap->count==ap->runs) {
+ change_stage=TRUE;
+ break;
+ }
+ set_temperature(moldyn,moldyn->t_ref+ap->dt);
+ ap->count+=1;
+ break;
+ case STAGE_CHAATTR:
+ chaatrib(moldyn,mdrun);
+ change_stage=TRUE;
+ break;
+ case STAGE_CHSATTR:
+ chsatrib(moldyn,mdrun);
+ change_stage=TRUE;
+ break;
+ default:
+ printf("%s unknwon stage type\n",ME);
+ break;
+ }
+
+ /* mark as executed */
+ stage->executed=TRUE;
+
+ /* change state */
+ if(change_stage==TRUE) {
+ printf("%s finished stage\n",ME);
+ if(list_next_f(sl)==L_NO_NEXT_ELEMENT) {
+ printf("%s no more stages\n",ME);
+ return 0;
+ }
+ steps=0;
+ tau=mdrun->timestep;
+ }
+
+ }
+ else {
+
+ /* averages */
+ if(mdrun->sattr&SATTR_AVGRST)
+ average_reset(moldyn);
+
+ }
+
+ /* continue simulation */
+ moldyn_add_schedule(moldyn,steps,tau);
+