started with crt ...

diff --git a/mdrun.c b/mdrun.c
index 29e5e487fcf9058048e30e522a5159ce1a20a993..130725d48fc70c74fcaf47146e963ff578928b15 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_CONSTRAINT_RELAXATION_TECHNIQUE:
+                       psize=sizeof(t_constraint_relaxation_technique);
+                       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,106 @@ 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;
+
+       extern u8 crt;
+       extern u8 *constraints;
+       extern double *trafo_angles;
+
+       t_atom *atom;
+       double dx,dy,dz;
+       
+       stage=mdrun->stage.current->data;
+       crtp=stage->params;
+
+       acount=0;
+
+       /* initial stuff */
+
+       if(crtp->count==0) {
+               printf("  crt init\n",acount);
+               // 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);
+                               crtp->acnt=acount;
+                               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.x=atof(wptr);
+                       wptr=strtok(NULL," \t");
+                       crtp->r_fin.y=atof(wptr);
+                       wptr=strtok(NULL," \t");
+                       crtp->r_fin.z=atof(wptr);
+                       // read constraints
+                       wptr=strtok(NULL," \t");
+                       constraints[acount]=atoi(wptr);
+                       wptr=strtok(NULL," \t");
+                       constraints[acount+1]=atoi(wptr);
+                       wptr=strtok(NULL," \t");
+                       constraints[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
+               crt=crtp->type;
+       }
+
+       /* crt routines: calculate displacement + set individual constraints */
+
+       for(i=0;i<moldyn->count;i++) {
+               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
+       }
+
+       return 0;
+}
+
 #define stage_print(m) if(!(stage->executed)) \
                                printf("%s",m)
 
@@ -1387,6 +1498,17 @@ 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) {
+                                       change_stage=TRUE;
+                                       break;
+                               }
+                               crt(moldyn,mdrun);
+                               crtp->count+=1;
+                               break;
                        default:
                                printf("%s unknwon stage type\n",ME);
                                break;

diff --git a/mdrun.h b/mdrun.h
index 0a5ceabaff4abcca73df7b7ecc4536c1fdc64fa7..102c3d2b6f3d098a216ca58758139a52a7bdeeee 100644 (file)
--- a/mdrun.h
+++ b/mdrun.h
@@ -54,6 +54,7 @@ typedef struct s_stage {
 #define STAGE_THERMAL_INIT                     0x10
 #define STAGE_DEL_ATOMS                                0x11
 #define STAGE_MODIFY_ATOMS                     0x12
+#define STAGE_CRT                              0x13
 
 typedef struct s_mdrun {
        char cfile[128];                        // config file
@@ -218,6 +219,15 @@ typedef struct s_fill_params {
        t_offset_params o_params;
 } t_fill_params;
 
+typedef struct s_crt_params {
+       u8 type;
+       char file[128];
+       t_3dvec *r_fin;
+       u8 *constraints;
+       int steps;
+       int count;
+} t_crt_params;
+
 /*
  * function prototypes
  */

diff --git a/moldyn.c b/moldyn.c
index bd29afc9a3e6853bded38acae91dfafdf15a662a..ccc0b56f99159f0bab42b18051116241355eb0f8 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -55,6 +55,13 @@ pthread_mutex_t *amutex;
 pthread_mutex_t emutex;
 #endif
 
+#ifdef CRT
+/* fully constrained relaxation technique - global pointers */
+u8 crt;
+u8 *constraints;
+double *trafo_angle;
+#endif
+
 /*
  * the moldyn functions
  */