completed crt implementation

diff --git a/Makefile b/Makefile
index fc5705a..fa40909 100644 (file)
--- 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 (file)
--- 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;i<moldyn->count;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 (file)
--- 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 (file)
--- 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<count;i++) {
+
                /* check whether fixed atom */
                if(atom[i].attr&ATOM_ATTR_FP)
                        continue;
+
                /* new positions */
                h=0.5/atom[i].mass;
-               v3_scale(&delta,&(atom[i].v),tau);
-#ifdef CONSTRAINT_110_5832
-               if(i==5832) {
-                       delta.y=-delta.x;
+
+               /* constraint relaxation */
+               if(crtt) {
+                       basis_trafo(&(atom[i].f),FORWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
+                       if(constraints[3*i])
+                               atom[i].f.x=0;
+                       if(constraints[3*i+1])
+                               atom[i].f.y=0;
+                       if(constraints[3*i+2])
+                               atom[i].f.z=0;
+                       basis_trafo(&(atom[i].f),BACKWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
                }
-#endif
-#ifdef JHDVHJDV
-#ifdef CONSTRAINT_11X_5832
-       if(i==5833) {
-               // second trafo
-               xt=delta.x;
-               yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
-               zt=delta.y*sin(-0.16935129)+delta.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
-               delta.x=xt;
-               delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
-               delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
-       }
-#endif
-#endif
+
 #ifndef QUENCH
+               v3_scale(&delta,&(atom[i].v),tau);
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
 #endif
                v3_scale(&delta,&(atom[i].f),h*tau_square);
-#ifdef CONSTRAINT_110_5832
-               if(i==5832) {
-                       delta.y=-delta.x;
-               }
-#endif
-#ifdef GHDJHBSJBSD
-#ifdef CONSTRAINT_11X_5832
-       if(i==5833) {
-               // second trafo
-               xt=delta.x;
-               yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
-               zt=delta.y*sin(-0.16935129)+delta.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
-               delta.x=xt;
-               delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
-               delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
-       }
-#endif
-#endif
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
                //check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
                check_per_bound(moldyn,&(atom[i].r));
@@ -2408,38 +2365,25 @@ int velocity_verlet(t_moldyn *moldyn) {
        albe_potential_force_calc(moldyn);
 #endif
 
-#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<count;i++) {
                /* check whether fixed atom */
                if(atom[i].attr&ATOM_ATTR_FP)
                        continue;
+
+               /* constraint relaxation */
+               if(crtt) {
+                       basis_trafo(&(atom[i].f),FORWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
+                       if(constraints[3*i])
+                               atom[i].f.x=0;
+                       if(constraints[3*i+1])
+                               atom[i].f.y=0;
+                       if(constraints[3*i+2])
+                               atom[i].f.z=0;
+                       basis_trafo(&(atom[i].f),BACKWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
+               }
+
                /* again velocities [actually v(t+tau)] */
                v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
                v3_add(&(atom[i].v),&(atom[i].v),&delta);