completed crt implementation

[physik/posic.git] / mdrun.c

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;