completed crt implementation

[physik/posic.git] / mdrun.c

diff --git a/mdrun.c b/mdrun.c
index 29e5e48..3266dd3 100644 (file)
--- a/mdrun.c
+++ b/mdrun.c
@@ -129,6 +129,9 @@ int add_stage(t_mdrun *mdrun,u8 type,void *params) {
                case STAGE_THERMAL_INIT:
                        psize=0;
                        break;
+               case STAGE_CRT:
+                       psize=sizeof(t_crt_params);
+                       break;
                default:
                        printf("%s unknown stage type: %02x\n",ME,type);
                        return -1;
@@ -178,6 +181,7 @@ int mdrun_parse_config(t_mdrun *mdrun) {
        t_set_timestep_params stsp;
        t_fill_params fp;
        t_del_atoms_params delp;
+       t_crt_params crtp;
 
        /* open config file */
        fd=open(mdrun->cfile,O_RDONLY);
@@ -212,6 +216,7 @@ int mdrun_parse_config(t_mdrun *mdrun) {
                memset(&stsp,0,sizeof(t_set_timestep_params));
                memset(&fp,0,sizeof(t_fill_params));
                memset(&delp,0,sizeof(t_del_atoms_params));
+               memset(&crtp,0,sizeof(t_crt_params));
 
                // get command + args
                wcnt=0;
@@ -713,6 +718,12 @@ int mdrun_parse_config(t_mdrun *mdrun) {
                                stsp.tau=atof(word[2]);
                                add_stage(mdrun,STAGE_SET_TIMESTEP,&stsp);
                        }
+                       else if(!strncmp(word[1],"crt",3)) {
+                               crtp.type=atoi(word[2]);
+                               crtp.steps=atoi(word[3]);
+                               strncpy(crtp.file,word[4],127);
+                               add_stage(mdrun,STAGE_CRT,&crtp);
+                       }
                        else {
                                printf("%s unknown stage type: %s\n",
                                       ME,word[1]);
@@ -1198,6 +1209,114 @@ int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
        return 0;
 }
 
+int crt(t_moldyn *moldyn,t_mdrun *mdrun) {
+
+       t_stage *stage;
+       t_crt_params *crtp;
+
+       int fd;
+       char line[128];
+       char *wptr;
+       int acount;
+       int ret;
+       void *ptr;
+
+       t_atom *atom;
+       t_3dvec disp;
+       double frac;
+       int i;
+       
+       stage=mdrun->stage.current->data;
+       crtp=stage->params;
+
+       acount=0;
+
+       /* initial stuff */
+
+       if(crtp->count==0) {
+               printf("  crt init\n");
+               // read final positions, constraints and do the alloc
+               fd=open(crtp->file,O_RDONLY);
+               if(fd<0) {
+                       perror("[mdrun] FATAL reading constraints file");
+                       return fd;
+               }
+               while(1) {
+                       ret=get_line(fd,line,128);
+                       // check for end of file
+                       if(ret<=0) {
+                               printf("  read %d atom positions\n",acount);
+                               if(acount!=moldyn->count)
+                                       printf("  atom count mismatch!!!\n");
+                               printf("\n");
+                               break;
+                       }
+                       // ignore # lines and \n
+                       if((line[0]=='#')|(ret==1))
+                               continue;
+                       // allocate new memory
+                       ptr=realloc(crtp->r_fin,(acount+1)*sizeof(t_3dvec));
+                       if(ptr==NULL) {
+                               perror("[mdrun] FATAL realloc crt positions");
+                               return -1;
+                       }
+                       crtp->r_fin=ptr;
+                       ptr=realloc(constraints,(acount+1)*3*sizeof(u8));
+                       if(ptr==NULL) {
+                               perror("[mdrun] FATAL realloc crt constraints");
+                               return -1;
+                       }
+                       constraints=ptr;
+                       // ignore type
+                       wptr=strtok(line," \t");
+                       // read x y z
+                       wptr=strtok(NULL," \t");
+                       crtp->r_fin[acount].x=atof(wptr);
+                       wptr=strtok(NULL," \t");
+                       crtp->r_fin[acount].y=atof(wptr);
+                       wptr=strtok(NULL," \t");
+                       crtp->r_fin[acount].z=atof(wptr);
+                       // read constraints
+                       wptr=strtok(NULL," \t");
+                       constraints[3*acount]=atoi(wptr);
+                       wptr=strtok(NULL," \t");
+                       constraints[3*acount+1]=atoi(wptr);
+                       wptr=strtok(NULL," \t");
+                       constraints[3*acount+2]=atoi(wptr);
+                       // done reading
+                       acount+=1;
+               }
+               // allocate trafo angles
+               trafo_angle=malloc(acount*2*sizeof(double));
+               if(trafo_angle==NULL) {
+                       perror("[mdrun] FATAL alloc trafo angles");
+                       return -1;
+               }
+               // set crt mode
+               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;
+               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;
+}
+
 #define stage_print(m) if(!(stage->executed)) \
                                printf("%s",m)
 
@@ -1218,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;
@@ -1387,6 +1507,20 @@ int mdrun_hook(void *ptr1,void *ptr2) {
                                thermal_init(moldyn,TRUE);
                                change_stage=TRUE;
                                break;
+                       case STAGE_CRT:
+                               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;
+                               }
+                               crt(moldyn,mdrun);
+                               crtp->count+=1;
+                               break;
                        default:
                                printf("%s unknwon stage type\n",ME);
                                break;
@@ -1430,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;