From ff44bcba9df0d021568abcb553d1d490c442f696 Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 19 May 2010 17:17:40 +0200 Subject: [PATCH] completed crt implementation --- Makefile | 10 ++- mdrun.c | 61 +++++++++++------ mdrun.h | 9 +++ moldyn.c | 202 ++++++++++++++++++++----------------------------------- 4 files changed, 128 insertions(+), 154 deletions(-) diff --git a/Makefile b/Makefile index fc5705a..fa40909 100644 --- a/Makefile +++ b/Makefile @@ -28,7 +28,6 @@ CFLAGS += -DLOWMEM_LISTS #CFLAGS += -DDSTART=50 -DDEND=60 -DDATOM=0 #CFLAGS += -DVDEBUG -#CFLAGS += -DCONSTRAINT_110_5832 #CFLAGS += -DQUENCH LDFLAGS = -lm @@ -46,9 +45,10 @@ SRC += potentials/lennard_jones.c potentials/harmonic_oscillator.c SRC += potentials/tersoff.c potentials/albe.c SRC += potentials/albe_fast.c -ALL = mdrun sic fluctuation_calc postproc pair_correlation_calc diffusion_calc +#ALL = mdrun sic fluctuation_calc postproc pair_correlation_calc diffusion_calc +ALL = mdrun fluctuation_calc postproc pair_correlation_calc diffusion_calc ALL += diffusion_calc_ver2 bond_analyze search_bonds visual_atoms -ALL += display_atom_data atom_match +ALL += display_atom_data atom_match msd_calc s2xyz all: $(ALL) @@ -76,6 +76,10 @@ display_atom_data: $(DEPS) atom_match: $(DEPS) +msd_calc: $(DEPS) + +s2xyz: $(DEPS) + .PHONY:clean clean: rm -vf $(ALL) *.o */*.o diff --git a/mdrun.c b/mdrun.c index 130725d..3266dd3 100644 --- a/mdrun.c +++ b/mdrun.c @@ -129,8 +129,8 @@ int add_stage(t_mdrun *mdrun,u8 type,void *params) { case STAGE_THERMAL_INIT: psize=0; break; - case STAGE_CONSTRAINT_RELAXATION_TECHNIQUE: - psize=sizeof(t_constraint_relaxation_technique); + case STAGE_CRT: + psize=sizeof(t_crt_params); break; default: printf("%s unknown stage type: %02x\n",ME,type); @@ -1221,12 +1221,10 @@ int crt(t_moldyn *moldyn,t_mdrun *mdrun) { int ret; void *ptr; - extern u8 crt; - extern u8 *constraints; - extern double *trafo_angles; - t_atom *atom; - double dx,dy,dz; + t_3dvec disp; + double frac; + int i; stage=mdrun->stage.current->data; crtp=stage->params; @@ -1236,7 +1234,7 @@ int crt(t_moldyn *moldyn,t_mdrun *mdrun) { /* initial stuff */ if(crtp->count==0) { - printf(" crt init\n",acount); + printf(" crt init\n"); // read final positions, constraints and do the alloc fd=open(crtp->file,O_RDONLY); if(fd<0) { @@ -1247,8 +1245,10 @@ int crt(t_moldyn *moldyn,t_mdrun *mdrun) { ret=get_line(fd,line,128); // check for end of file if(ret<=0) { - printf(" -> read %d atom positions\n",acount); - crtp->acnt=acount; + printf(" read %d atom positions\n",acount); + if(acount!=moldyn->count) + printf(" atom count mismatch!!!\n"); + printf("\n"); break; } // ignore # lines and \n @@ -1271,18 +1271,18 @@ int crt(t_moldyn *moldyn,t_mdrun *mdrun) { wptr=strtok(line," \t"); // read x y z wptr=strtok(NULL," \t"); - crtp->r_fin.x=atof(wptr); + crtp->r_fin[acount].x=atof(wptr); wptr=strtok(NULL," \t"); - crtp->r_fin.y=atof(wptr); + crtp->r_fin[acount].y=atof(wptr); wptr=strtok(NULL," \t"); - crtp->r_fin.z=atof(wptr); + crtp->r_fin[acount].z=atof(wptr); // read constraints wptr=strtok(NULL," \t"); - constraints[acount]=atoi(wptr); + constraints[3*acount]=atoi(wptr); wptr=strtok(NULL," \t"); - constraints[acount+1]=atoi(wptr); + constraints[3*acount+1]=atoi(wptr); wptr=strtok(NULL," \t"); - constraints[acount+2]=atoi(wptr); + constraints[3*acount+2]=atoi(wptr); // done reading acount+=1; } @@ -1293,17 +1293,25 @@ int crt(t_moldyn *moldyn,t_mdrun *mdrun) { return -1; } // set crt mode - crt=crtp->type; + crtt=crtp->type; } /* crt routines: calculate displacement + set individual constraints */ + printf(" crt step %d of %d in total\n\n",crtp->count+1,crtp->steps); + for(i=0;icount;i++) { + // calc displacements atom=moldyn->atom; - dx=atom[i].r.x-crtp->r_fin[i].x; - dy=atom[i].r.y-crtp->r_fin[i].y; - dz=atom[i].r.z-crtp->r_fin[i].z; - // HIER WEITER + v3_sub(&disp,&(crtp->r_fin[i]),&(atom[i].r)); + // angles + trafo_angle[2*i]=atan2(disp.x,disp.y); + trafo_angle[2*i+1]=-atan2(disp.z, + sqrt(disp.x*disp.x+disp.y*disp.y)); + // move atoms + frac=1.0/(crtp->steps-crtp->count); + v3_scale(&disp,&disp,frac); + v3_add(&(atom[i].r),&(atom[i].r),&disp); } return 0; @@ -1329,6 +1337,7 @@ int mdrun_hook(void *ptr1,void *ptr2) { t_set_temp_params *stp; t_set_timestep_params *stsp; t_fill_params *fp; + t_crt_params *crtp; moldyn=ptr1; mdrun=ptr2; @@ -1499,10 +1508,13 @@ int mdrun_hook(void *ptr1,void *ptr2) { change_stage=TRUE; break; case STAGE_CRT: - stage_print(" -> constraint relaxation") + stage_print(" -> constraint relaxation"); stage_print(" technique\n\n"); crtp=stage->params; if(crtp->count==crtp->steps) { + free(constraints); + free(trafo_angle); + free(crtp->r_fin); change_stage=TRUE; break; } @@ -1552,6 +1564,11 @@ int main(int argc,char **argv) { memset(&mdrun,0,sizeof(t_mdrun)); memset(&moldyn,0,sizeof(t_moldyn)); + /* init crt variables */ + crtt=0; + constraints=NULL; + trafo_angle=NULL; + /* parse arguments */ if(mdrun_parse_argv(&mdrun,argc,argv)<0) return -1; diff --git a/mdrun.h b/mdrun.h index 102c3d2..b8ed70c 100644 --- a/mdrun.h +++ b/mdrun.h @@ -228,6 +228,15 @@ typedef struct s_crt_params { int count; } t_crt_params; +/* + * extern variables + */ + +// constraint relaxation technique +extern u8 crtt; +extern u8 *constraints; +extern double *trafo_angle; + /* * function prototypes */ diff --git a/moldyn.c b/moldyn.c index ccc0b56..a6f5c45 100644 --- a/moldyn.c +++ b/moldyn.c @@ -55,12 +55,10 @@ pthread_mutex_t *amutex; pthread_mutex_t emutex; #endif -#ifdef CRT /* fully constrained relaxation technique - global pointers */ -u8 crt; +u8 crtt; u8 *constraints; double *trafo_angle; -#endif /* * the moldyn functions @@ -85,14 +83,9 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { pthread_mutex_init(&emutex,NULL); #endif -#ifdef CONSTRAINT_110_5832 - printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n"); - printf("\n\n\n!! -- constraints enabled -- !!\n\n\n"); -#endif -#ifdef CONSTRAINT_11X_5832 - printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n"); - printf("\n\n\n!! -- constraints enabled -- !!\n\n\n"); -#endif + if(crtt) + printf("USING CRT\n"); + return 0; } @@ -2267,6 +2260,43 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n", return 0; } +/* basis trafo */ + +#define FORWARD 0 +#define BACKWARD 1 + +int basis_trafo(t_3dvec *r,u8 dir,double z,double x) { + + t_3dvec tmp; + + if(dir==FORWARD) { + if(z!=0.0) { + v3_copy(&tmp,r); + r->x=cos(z)*tmp.x-sin(z)*tmp.y; + r->y=sin(z)*tmp.x+cos(z)*tmp.y; + } + if(x!=0.0) { + v3_copy(&tmp,r); + r->y=cos(x)*tmp.y-sin(x)*tmp.z; + r->z=sin(x)*tmp.y+cos(x)*tmp.z; + } + } + else { + if(x!=0.0) { + v3_copy(&tmp,r); + r->y=cos(-x)*tmp.y-sin(-x)*tmp.z; + r->z=sin(-x)*tmp.y+cos(-x)*tmp.z; + } + if(z!=0.0) { + v3_copy(&tmp,r); + r->x=cos(-z)*tmp.x-sin(-z)*tmp.y; + r->y=sin(-z)*tmp.x+cos(-z)*tmp.y; + } + } + + return 0; +} + /* velocity verlet */ int velocity_verlet(t_moldyn *moldyn) { @@ -2275,113 +2305,40 @@ int velocity_verlet(t_moldyn *moldyn) { double tau,tau_square,h; t_3dvec delta; t_atom *atom; -#ifdef CONSTRAINT_11X_5832 - double xt,yt,zt; - double xtt,ytt,ztt; -#endif atom=moldyn->atom; count=moldyn->count; tau=moldyn->tau; tau_square=moldyn->tau_square; -#ifdef CONSTRAINT_110_5832 - if(count==5833) { - atom[5832].f.x=0.5*(atom[5832].f.x-atom[5832].f.y); - atom[5832].f.y=-atom[5832].f.x; - } -#endif -#ifdef CONSTRAINT_11X_5832 - if(count==5833) { - // second trafo - xt=atom[5832].f.x; - yt=atom[5832].f.y*cos(-0.16935129)-atom[5832].f.z*sin(-0.16935129); - zt=atom[5832].f.y*sin(-0.16935129)+atom[5832].f.z*cos(-0.16935129); - // first trafo - xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0); - ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0); - ztt=zt; - // apply constraints - ytt=0.0; - // first trafo backwards - xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0); - yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0); - zt=ztt; - // second trafo backwards - atom[5832].f.x=xt; - atom[5832].f.y=yt*cos(-0.16935129)+zt*sin(-0.16935129); - atom[5832].f.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129); - } -#endif for(i=0;i