NAMD
sortreplicas.c
Go to the documentation of this file.
1 
2 #include "largefiles.h" /* must be first! */
3 
4 #include "vmdplugin.h"
5 
6 extern int molfile_dcdplugin_init(void);
7 extern int molfile_dcdplugin_register(void *, vmdplugin_register_cb);
8 extern int molfile_dcdplugin_fini(void);
9 /*
10 extern int molfile_jsplugin_init(void);
11 extern int molfile_jsplugin_register(void *, vmdplugin_register_cb);
12 extern int molfile_jsplugin_fini(void);
13 */
14 
15 #include "molfile_plugin.h"
16 
17 static int register_cb(void *v, vmdplugin_t *p) {
18  *((vmdplugin_t **)v) = p;
19  return VMDPLUGIN_SUCCESS;
20 }
21 
22 #include <sys/types.h>
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <string.h>
26 #include <errno.h>
27 
28 int main(int argc, char **argv) {
29 
30  molfile_timestep_t frame;
31  molfile_plugin_t *plugin;
32  char *output_root;
33  char *filename;
34  int num_replicas;
35  int runs_per_frame;
36  long long int final_step = -1;
37  long long int checkstep = -1;
38  int colvars;
39  FILE **hist_in;
40  FILE **hist_out;
41  FILE **colv_in;
42  FILE **colv_out;
43  void **traj_in;
44  void **traj_out;
45  int natoms=MOLFILE_NUMATOMS_UNKNOWN;
46  int i, i_run;
47 
50 
51  if ( argc < 4 || argc > 5 ) {
52  fprintf(stderr, "args: <job_output_root> <num_replicas> <runs_per_frame> [final_step]\n");
53  exit(-1);
54  }
55  output_root = argv[1];
56  num_replicas = atoi(argv[2]);
57  runs_per_frame = atoi(argv[3]);
58  if ( argc > 4 ) {
59  sscanf(argv[4], "%lld", &final_step);
60  }
61 
62  filename = (char*) malloc(strlen(output_root)+100);
63  hist_in = (FILE**) malloc(num_replicas*sizeof(FILE*));
64  hist_out = (FILE**) malloc(num_replicas*sizeof(FILE*));
65  colv_in = (FILE**) malloc(num_replicas*sizeof(FILE*));
66  colv_out = (FILE**) malloc(num_replicas*sizeof(FILE*));
67  traj_in = (void**) malloc(num_replicas*sizeof(FILE*));
68  traj_out = (void**) malloc(num_replicas*sizeof(FILE*));
69 
70  for ( i=0; i<num_replicas; ++i ) {
71  char *root_end;
72  if ( strstr(output_root,"%s") ) {
73  char istr[10];
74  sprintf(istr,"%d",i);
75  sprintf(filename,output_root,istr);
76  } else {
77  sprintf(filename,output_root,i);
78  }
79  root_end = filename + strlen(filename);
80 
81  sprintf(root_end,".%d.history",i);
82  hist_in[i] = fopen(filename,"r");
83  if ( ! hist_in[i] ) {
84  fprintf(stderr, "error opening input file %s: %s\n",
85  filename, strerror(errno));
86  exit(-1);
87  }
88  sprintf(root_end,".%d.dcd",i);
89  traj_in[i] = plugin->open_file_read(filename,"dcd",&natoms);
90  if ( ! traj_in[i] ) {
91  fprintf(stderr, "error opening input file %s: %s\n",
92  filename, strerror(errno));
93  exit(-1);
94  }
95  sprintf(root_end,".%d.colvars.traj",i);
96  colv_in[i] = fopen(filename,"r");
97  if ( colv_in[i] ) {
98  if ( i == 0 ) {
99  printf("Found first input colvars trajectory file %s.\n", filename);
100  colvars = 1;
101  } else if ( ! colvars ) {
102  fprintf(stderr, "missing input colvars trajectory files before %s\n", filename);
103  exit(-1);
104  }
105  } else {
106  if ( i == 0 ) {
107  colvars = 0;
108  } else if ( colvars ) {
109  fprintf(stderr, "error opening input colvars trajectory file %s: %s\n",
110  filename, strerror(errno));
111  exit(-1);
112  }
113  }
114  }
115 
116  for ( i=0; i<num_replicas; ++i ) {
117  char *root_end;
118  if ( strstr(output_root,"%s") ) {
119  char istr[10];
120  sprintf(istr,"%d",i);
121  sprintf(filename,output_root,istr);
122  } else {
123  sprintf(filename,output_root,i);
124  }
125  root_end = filename + strlen(filename);
126 
127  sprintf(root_end,".%d.sort.history",i);
128  hist_out[i] = fopen(filename,"w");
129  if ( ! hist_out[i] ) {
130  fprintf(stderr, "error opening output file %s: %s\n",
131  filename, strerror(errno));
132  exit(-1);
133  }
134  sprintf(root_end,".%d.sort.dcd",i);
135  traj_out[i] = plugin->open_file_write(filename,"dcd",natoms);
136  if ( ! traj_out[i] ) {
137  fprintf(stderr, "error opening output file %s: %s\n",
138  filename, strerror(errno));
139  exit(-1);
140  }
141  if ( colvars ) {
142  sprintf(root_end,".%d.sort.colvars.traj",i);
143  colv_out[i] = fopen(filename,"w");
144  if ( ! colv_out[i] ) {
145  fprintf(stderr, "error opening output file %s: %s\n",
146  filename, strerror(errno));
147  exit(-1);
148  }
149  }
150  }
151 
152  frame.coords = (float*) malloc(3*natoms*sizeof(float));
153  frame.velocities = (float*) NULL;
154 
155 #define LINE_MAX 10000
156 
157  i_run = 0;
158  for ( ; 1; ++i_run ) { /* loop until read fails */
159  char line[LINE_MAX];
160  for ( i=0; i<num_replicas; ++i ) {
161  char *r;
162  char sav;
163  int f1,f2;
164  int rc;
165  int rep_id = -1;
166  long long int step;
167  r = fgets(line, LINE_MAX, hist_in[i]);
168  if ( ! r ) { break; }
169  rc = sscanf(line, "%lld %n%d%n", &step, &f1, &rep_id, &f2);
170  if ( rc != 2 ) {
171  fprintf(stderr,"Format error for replica %d at line %d: %s",
172  i, i_run, line);
173  exit(-1);
174  }
175  if ( i == 0 ) {
176  if ( step <= checkstep ) {
177  fprintf(stderr,"Step out of order for replica %d at line %d: %s",
178  i, i_run, line);
179  exit(-1);
180  }
181  checkstep = step;
182  if ( final_step >= 0 && checkstep > final_step ) {
183  printf("Stopping after final step %lld.\n", final_step);
184  break;
185  }
186  } else if ( step != checkstep ) {
187  fprintf(stderr,"Step mismatch for replica %d at line %d: %s",
188  i, i_run, line);
189  exit(-1);
190  }
191  if ( rep_id < 0 || rep_id >= num_replicas ) {
192  fprintf(stderr,"Invalid replica ID for replica %d at line %d: %s",
193  i, i_run, line);
194  exit(-1);
195  }
196  sav = line[f1];
197  line[f1] = 0;
198  fprintf(hist_out[rep_id],"%s%d%s",line,i,line+f2);
199  line[f1] = sav;
200  if ( colvars ) {
201  long long int oldcstep = -1;
202  while ( 1 ) {
203  long long int cstep;
204  char cline[LINE_MAX];
205 #ifdef WIN32
206  __int64 oldpos = _ftelli64(colv_in[i]);
207 #else
208  off_t oldpos = ftello(colv_in[i]);
209 #endif
210  r = fgets(cline, LINE_MAX, colv_in[i]);
211  if ( ! r ) { break; }
212  if ( cline[0] == '#' ) {
213  fprintf(colv_out[rep_id],"%s",cline);
214  continue;
215  }
216  rc = sscanf(cline, "%lld", &cstep);
217  if ( rc != 1 ) {
218  fprintf(stderr,"Format error in colvar trajectory for replica %d: %s",
219  i, cline);
220  exit(-1);
221  }
222  if ( cstep == oldcstep ) continue; /* filter out repeats */
223  if ( cstep < oldcstep ) {
224  fprintf(stderr,"Step out of order in colvar trajectory for replica %d: %s",
225  i, cline);
226  exit(-1);
227  }
228  if ( cstep > step ) {
229 #ifdef WIN32
230  _fseeki64(colv_in[i], oldpos, SEEK_SET);
231 #else
232  fseeko(colv_in[i], oldpos, SEEK_SET);
233 #endif
234  break;
235  }
236  if ( i_run != 0 || oldcstep != -1 ) { /* skip first entry */
237  fprintf(colv_out[rep_id],"%s",cline);
238  }
239  oldcstep = cstep;
240  }
241  }
242  if ( (i_run+1) % runs_per_frame ) continue;
243  rc = plugin->read_next_timestep(traj_in[i],natoms,&frame);
244  if ( rc == MOLFILE_SUCCESS ) {
245  plugin->write_timestep(traj_out[rep_id],&frame);
246  } else {
247  fprintf(stderr,"Unable to read frame for replica %d at line %d: %s",
248  i, i_run, line);
249  break;
250  }
251  }
252  if ( i < num_replicas ) {
253  printf("Processed %d runs.\n",i_run);
254  if ( i ) fprintf(stderr,"Uneven input lengths for replica %d at line %d: %s",
255  i, i_run, line);
256  break;
257  }
258  }
259 
260  free(frame.coords);
261 
262  for ( i=0; i<num_replicas; ++i ) {
263  if ( fclose(hist_in[i]) ) {
264  fprintf(stderr, "error closing history input file %d: %s\n", i, strerror(errno));
265  }
266  plugin->close_file_read(traj_in[i]);
267  if ( fclose(hist_out[i]) ) {
268  fprintf(stderr, "error closing history output file %d: %s\n", i, strerror(errno));
269  }
270  plugin->close_file_write(traj_out[i]);
271  if ( colvars ) {
272  if ( fclose(colv_in[i]) ) {
273  fprintf(stderr, "error closing colvars input file %d: %s\n", i, strerror(errno));
274  }
275  if ( fclose(colv_out[i]) ) {
276  fprintf(stderr, "error closing colvars output file %d: %s\n", i, strerror(errno));
277  }
278  }
279  }
280 
282  exit(0);
283 }
284 
int molfile_dcdplugin_init(void)
#define LINE_MAX
int register_cb(void *v, vmdplugin_t *p)
Definition: PluginIOMgr.C:6
int main(int argc, char *argv[])
Definition: diffbinpdb.c:15
int molfile_dcdplugin_fini(void)
int molfile_dcdplugin_register(void *, vmdplugin_register_cb)