X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=mdrun.c;h=6f7a2f02a0704491b5d8bf50762031e0149711ba;hb=d7f67c88195ab155f2737e57cc5e81973d3feb0c;hp=4a2822e86d10c5d967bcc55ce00de6964c26ba8b;hpb=61d24d027511c3e96b2f94558dc1b31c72725767;p=physik%2Fposic.git

diff --git a/mdrun.c b/mdrun.c
index 4a2822e..6f7a2f0 100644
--- a/mdrun.c
+++ b/mdrun.c
@@ -62,6 +62,75 @@ int mdrun_parse_argv(t_mdrun *mdrun,int argc,char **argv) {
 	return 0;
 }
 
+del_stages(t_mdrun *mdrun) {
+
+	t_list *sl;
+	t_stage *stage;
+
+	sl=mdrun->stage;
+
+	list_reset_f(sl);
+
+	if(sl->start==NULL)
+		return 0;
+
+	do {
+		stage=sl->current->data;
+		free(stage->params);
+		free(stage);
+	} while(list_next_f(sl)!=L_NO_NEXT_ELEMENT);
+
+	return 0;
+}
+
+add_stage(t_mdrun *mdrun,u8 type,void *params) {
+
+	int psize;
+
+	t_stage *stage;
+
+	switch(type) {
+		case STAGE_INSERT_ATOMS:
+			psize=sizeof(t_insert_atoms_params);
+			break;
+		case STAGE_CONTINUE:
+			psize=sizeof(t_continue_params);
+			break;
+		case STAGE_ANNEAL:
+			psize=sizeof(t_anneal_params);
+			break;
+		case STAGE_CHAATTR:
+			psize=sizeof(t_chaattr_params);
+			break;
+		case STAGE_CHSATTR:
+			psize=sizeof(t_chsattr_params);
+			break;
+		default:
+
+	}
+
+	stage=malloc(sizeof(t_stage));
+	if(stage==NULL) {
+		perror("[mdrun] malloc (add stage)");
+		return -1;
+	}
+
+	stage->type=type;
+	stage->executed=FALSE;
+
+	stage->params=malloc(psize);
+	if(stage->params==NULL) {
+		perror("[mdrun] malloc (stage params)");
+		return -1;
+	}
+
+	memcpy(stage->params,params,psize);
+
+	list_add_immediate_f(mdrun->stage,stage);
+
+	return 0;
+}
+
 int mdrun_parse_config(t_mdrun *mdrun) {
 
 	int fd,ret;
@@ -70,6 +139,13 @@ int mdrun_parse_config(t_mdrun *mdrun) {
 	char *wptr;
 	char word[16][32];
 	int wcnt;
+	int i;
+
+	t_insert_atoms_params iap;
+	t_continue_params cp;
+	t_anneal_params ap;
+	t_chaattr_params cap;
+	t_chsattr_params csp;
 
 	/* open config file */
 	fd=open(mdrun->cfile,O_RDONLY);
@@ -186,17 +262,51 @@ int mdrun_parse_config(t_mdrun *mdrun) {
 					return -1;
 			}
 		}
-		else if(!strncmp(word[0],"filllc",6)) {
-			mdrun->lx=atoi(word[1]);
-			mdrun->ly=atoi(word[2]);
-			mdrun->lz=atoi(word[3]);
+		else if(!strncmp(word[0],"fill",6)) {
+			// only lc mode by now
+			mdrun->lx=atoi(word[2]);
+			mdrun->ly=atoi(word[3]);
+			mdrun->lz=atoi(word[4]);
 		}
-		else if(!strncmp(word[0],"aattrib",7)) {
+		else if(!strncmp(word[0],"aattr",5)) {
 			// for aatrib line we need a special stage
 			// containing one schedule of 0 loops ...
-			printf("%s aattrib: %s\n",ME,word[1]);
+			if(!strncmp(word[1],"all",3)) {
+				cap.type=
+
+
+				HIER WEITER
+			}
+			for(i=0;i<strlen(wptr)) {
+				switch(word[2][i]) {
+					case 'b':
+						cap.attr|=ATOM_ATTR_VB;
+						break;
+					case 'h':
+						cap.attr|=ATOM_ATTR_HB;
+						break;
+					case 'v':
+						cap.attr|=ATOM_ATTR_VA;
+						break;
+					case 'f':
+						cap.attr|=ATOM_ATTR_FP;
+						break;
+					case '1':
+						cap.attr|=ATOM_ATTR_1BP;
+						break;
+					case '2':
+						cap.attr|=ATOM_ATTR_2BP;
+						break;
+					case '3':
+						cap.attr|=ATOM_ATTR_3BP;
+						break;
+					default:
+						break;
+				}
+			}
+			add_stage(mdrun,STAGE_CHAATTR,&cap);
 		}
-		else if(!strncmp(word[0],"eattrib",7)) {
+		else if(!strncmp(word[0],"ectrl",5)) {
 			mdrun->p_tau=atof(word[1]);
 			mdrun->t_tau=atof(word[2]);
 		}
@@ -221,11 +331,303 @@ int mdrun_parse_config(t_mdrun *mdrun) {
 			printf("%s stage: %s\n",ME,word[1]);
 		}
 		else {
-			printf("%s unknown command %s, skipped!\n",ME,wptr);
+			printf("%s unknown command %s, skipped!\n",ME,word[0]);
 			continue;
 		}
 	}
 
+	/* 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);
+
 	return 0;
 }
 
@@ -242,9 +644,15 @@ int main(int argc,char **argv) {
 	if(mdrun_parse_argv(&mdrun,argc,argv)<0)
 		return -1;
 
+	/* initialize list system */
+	list_init_f(&(mdrun.stage));
+
 	/* parse config file */
 	mdrun_parse_config(&mdrun);
 
+	/* reset the stage list */
+	list_reset_f(&(mdrun.stage));
+
 	/* sanity check! */
 
 	/* prepare simulation */
@@ -326,13 +734,22 @@ int main(int argc,char **argv) {
 	moldyn_set_log(&moldyn,CREATE_REPORT,0);
 	set_avg_skip(&moldyn,mdrun.avgskip);
 
+	/* pt scaling */
+	if(mdrun.p_tau!=0)
+		set_p_scale(&moldyn,P_SCALE_BERENDSEN,mdrun.p_tau);
+	if(mdrun.t_tau!=0)
+		set_t_scale(&moldyn,T_SCALE_BERENDSEN,mdrun.t_tau);
+
 	/* prepare the hook function */
+	moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);
 
 	/* run the simulation */
 	moldyn_integrate(&moldyn);
 
 	/* shutdown */
 	moldyn_shutdown(&moldyn);
+	del_stages(&mdrun);
+	list_destroy_f(&(mdrun.stage));
 
 	return 0;
 }