NAMD
Output.C
Go to the documentation of this file.
1 
7 /*
8  This object outputs the data collected on the master node
9 */
10 
11 #include "largefiles.h" // must be first!
12 
13 #include <string.h>
14 #include <stdlib.h>
15 
16 #include "InfoStream.h"
17 #include "IMDOutput.h"
18 #include "Output.h"
19 #include "dcdlib.h"
20 #include "strlib.h"
21 #include "Molecule.h"
22 #include "Node.h"
23 #include "Parameters.h"
24 #include "PDB.h"
25 #include "SimParameters.h"
26 #include "Vector.h"
27 #include "structures.h"
28 #include "MStream.h"
29 #include "Communicate.h"
30 #include "PatchMap.h"
31 #include "PatchMap.inl"
32 #include "ScriptTcl.h"
33 #include "Lattice.h"
34 #include "DataExchanger.h"
35 #include <fcntl.h>
36 #include <sys/stat.h>
37 #include "PatchData.h"
38 #ifdef WIN32
39 #include <io.h>
40 #define access(PATH,MODE) _access(PATH,00)
41 #endif
42 static const short maxDCD=257; //1 byte of userdcd plus standard
43 #if defined(WIN32) && !defined(__CYGWIN__)
44 #define PATHSEPSTR "\\"
45 #define MKDIR(X) mkdir(X)
46 #else
47 #define PATHSEPSTR "/"
48 #define MKDIR(X) mkdir(X,0777)
49 #endif
50 
51 #define NAMD_open NAMD_open64
52 #define NAMD_write NAMD_write64
53 #define NAMD_close NAMD_close64
54 
55 #ifndef O_LARGEFILE
56 #define O_LARGEFILE 0x0
57 #endif
58 
59 
60 // same as open, only does error checking internally
61 int NAMD_open(const char *fname) {
62  int fd;
63 
64  // open the file and die if the open fails
65 #ifdef WIN32
66  while ( (fd = _open(fname, O_WRONLY|O_CREAT|O_EXCL|O_BINARY|O_LARGEFILE,_S_IREAD|_S_IWRITE)) < 0) {
67 #else
68 #ifdef NAMD_NO_O_EXCL
69  while ( (fd = open(fname, O_WRONLY|O_CREAT|O_TRUNC|O_LARGEFILE,
70 #else
71  while ( (fd = open(fname, O_WRONLY|O_CREAT|O_EXCL|O_LARGEFILE,
72 #endif
73  S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH)) < 0) {
74 #endif
75  if ( errno != EINTR ) {
76  char errmsg[1024];
77  sprintf(errmsg, "Unable to open binary file %s", fname);
78  NAMD_err(errmsg);
79  }
80  }
81 
82  return fd;
83 }
84 
85 // same as write, only does error checking internally
86 void NAMD_write(int fd, const char *buf, size_t count, const char *errmsg="NAMD_write64") {
87  double firsttime = 0.;
88  while ( count ) {
89 #if defined(WIN32) && !defined(__CYGWIN__)
90  long retval = _write(fd,buf,count);
91 #else
92  ssize_t retval = write(fd,buf,count);
93 #endif
94  if ( retval < 0 && errno == EINTR ) retval = 0;
95  if ( retval < 0 && errno == ENOMEM ) {
96  if ( firsttime == 0. ) firsttime = CmiWallTimer();
97  if ( (CmiWallTimer() - firsttime) < 300. ) retval = 0;
98  }
99  if ( retval < 0 ) NAMD_err(errmsg);
100  if ( retval > count ) NAMD_bug("extra bytes written in NAMD_write64()");
101  buf += retval;
102  count -= retval;
103  }
104  if ( firsttime != 0. ) {
105  iout << iWARN << errmsg << ": NAMD_write64() retried for " << (CmiWallTimer() - firsttime) << " seconds.\n" << endi;
106  }
107 }
108 
109 // same as close, only does error checking internally
110 void NAMD_close(int fd, const char *fname) {
111 #ifdef WIN32
112  while ( _close(fd) ) {
113 #else
114  while ( close(fd) ) {
115 #endif
116  if ( errno != EINTR ) {
117  char errmsg[1024];
118  sprintf(errmsg, "Error on closing file %s", fname);
119  NAMD_err(errmsg);
120  }
121  }
122 }
123 
124 
125 #define seek_dcdfile NAMD_seek
126 
127 // These make the NAMD 1 names work in NAMD 2
128 #define namdMyNode Node::Object()
129 #define simParams simParameters
130 #define pdbData pdb
131 
132 static void lattice_to_unitcell(const Lattice *lattice, double *unitcell) {
133 
134  unitcell[0] = unitcell[2] = unitcell[5] = 0.0;
135  unitcell[1] = unitcell[3] = unitcell[4] = 0.0;
136 
137  if (lattice) {
138  const Vector &a=lattice->a();
139  const Vector &b=lattice->b();
140  const Vector &c=lattice->c();
141  unitcell[0] = (lattice->a_p()) ? a.length() : 0.0;
142  unitcell[2] = (lattice->b_p()) ? b.length() : 0.0;
143  unitcell[5] = (lattice->c_p()) ? c.length() : 0.0;
144  double cosAB = (lattice->a_p() && lattice->b_p() ) ?
145  (a*b)/(unitcell[0]*unitcell[2]) : 0.0;
146  double cosAC = (lattice->a_p() && lattice->c_p() ) ?
147  (a*c)/(unitcell[0]*unitcell[5]) : 0.0;
148  double cosBC = (lattice->b_p() && lattice->c_p() ) ?
149  (b*c)/(unitcell[2]*unitcell[5]) : 0.0;
150  if (cosAB > 1.0) cosAB = 1.0; else if (cosAB < -1.0) cosAB = -1.0;
151  if (cosAC > 1.0) cosAC = 1.0; else if (cosAC < -1.0) cosAC = -1.0;
152  if (cosBC > 1.0) cosBC = 1.0; else if (cosBC < -1.0) cosBC = -1.0;
153  unitcell[1] = cosAB;
154  unitcell[3] = cosAC;
155  unitcell[4] = cosBC;
156  }
157 }
158 
159 
160 /************************************************************************/
161 /* */
162 /* FUNCTION Output */
163 /* */
164 /* This is the constructor for the Ouput class. It just sets */
165 /* up some values for the VMD connection. */
166 /* */
167 /************************************************************************/
168 
169 Output::Output() : replicaDcdActive(0) { }
170 
171 /* END OF FUNCTION Output */
172 
173 /************************************************************************/
174 /* */
175 /* FUNCTION ~Output */
176 /* */
177 /************************************************************************/
178 
179 Output::~Output() { }
180 
181 /* END OF FUNCTION ~Output */
182 
183 /************************************************************************/
184 /* */
185 /* FUNCTION coordinate */
186 /* */
187 /* INPUTS: */
188 /* timestep - Timestep coordinates were accumulated for */
189 /* n - This is the number of coordinates accumulated. */
190 /* vel - Array of Vectors containing the velocities */
191 /* */
192 /* This function receives the coordinates accumulated for a given */
193 /* timestep from the Collect object and calls the appropriate output */
194 /* functions. ALL routines used to output coordinates information */
195 /* should be called from here. */
196 /* */
197 /************************************************************************/
198 
199 std::pair<int,int> Output::coordinateNeeded(int timestep)
200 {
201  SimParameters *simParams = Node::Object()->simParameters;
202 
203  if(simParams->benchTimestep) return std::make_pair(0,0);
204 
205  int positionsNeeded = 0;
206  uint16 dcdSelectionTag = 0;
207  if ( timestep >= 0 ) {
208 
209  // Output a DCD trajectory
210  if ( simParams->dcdFrequency &&
211  ((timestep % simParams->dcdFrequency) == 0) )
212  { positionsNeeded |= 1; }
213 
214  // Output a DCD Selection trajectory
215 
216  if (simParams->dcdSelectionOn)
217  {
218  Molecule *mol = Node::Object()->molecule;
219  for(int index=0; index<16;++index)
220  {
221  int frequency=mol->dcdSelectionParams[index].frequency;
222  uint16 tag=mol->dcdSelectionParams[index].tag;
223  if(frequency && (timestep % frequency) == 0)
224  {
225  if(dcdSelectionTag>0)
226  positionsNeeded |=1; //multi selection means all positions
227  positionsNeeded |= 4; dcdSelectionTag = tag;
228  }
229  }
230  }
231  // Output a restart file
232  if ( simParams->restartFrequency &&
233  ((timestep % simParams->restartFrequency) == 0) )
234  { positionsNeeded |= 2; }
235 
236  // Iteractive MD
237  if ( simParams->IMDon &&
238  ( ((timestep % simParams->IMDfreq) == 0) ||
239  (timestep == simParams->firstTimestep) ) )
240  { positionsNeeded |= 1; }
241 
242  }
243 
244  // Output final coordinates
245  if (timestep == FILE_OUTPUT || timestep == END_OF_RUN ||
246  timestep == EVAL_MEASURE)
247  {
248  positionsNeeded |= 2;
249  }
250  return std::make_pair(positionsNeeded,dcdSelectionTag);
251 }
252 
253 template <class xVector, class xDone>
254 void wrap_coor_int(xVector *coor, Lattice &lattice, xDone *done) {
255  SimParameters *simParams = Node::Object()->simParameters;
256  if ( *done ) return;
257  *done = 1;
258  if ( ! ( simParams->wrapAll || simParams->wrapWater ) ) return;
259  const int wrapNearest = simParams->wrapNearest;
260  const int wrapAll = simParams->wrapAll;
261  Molecule *molecule = Node::Object()->molecule;
262  int n = molecule->numAtoms;
263  int i;
264 #ifndef MEM_OPT_VERSION
265  Position *con = new Position[n];
266  for ( i = 0; i < n; ++i ) {
267  con[i] = 0;
268  int ci = molecule->get_cluster(i);
269  con[ci] += coor[i];
270  }
271  for ( i = 0; i < n; ++i ) {
272  if ( ! wrapAll && ! molecule->is_water(i) ) continue;
273  int ci = molecule->get_cluster(i);
274  if ( ci == i ) {
275  Vector coni = con[i] / molecule->get_clusterSize(i);
276  Vector trans = ( wrapNearest ?
277  lattice.wrap_nearest_delta(coni) : lattice.wrap_delta(coni) );
278  con[i] = trans;
279  }
280  coor[i] = coor[i] + con[ci];
281  }
282  delete [] con;
283 #endif
284 }
285 
286 // handle index translation between global and selection spaces
287 template <class xVector, class xDone>
288  void wrap_coor_int_dcd_selection(xVector *coor, Lattice &lattice, xDone *done, int index) {
289  SimParameters *simParams = Node::Object()->simParameters;
290  if ( *done ) return;
291  *done = 1;
292  if ( ! ( simParams->wrapAll || simParams->wrapWater ) ) return;
293  const int wrapNearest = simParams->wrapNearest;
294  const int wrapAll = simParams->wrapAll;
295  Molecule *molecule = Node::Object()->molecule;
296  const int n = molecule->get_dcd_selection_size(index);
297  int i;
298 #ifndef MEM_OPT_VERSION
299  Position *con = new Position[n];
300  for ( i = 0; i < n; ++i ) {
301  int molIndex = molecule->get_atom_index_from_dcd_selection(index, i);
302  con[i] = 0;
303  int ci = molecule->get_cluster(molIndex);
304  int dcdci = molecule->get_dcd_selection_index_from_atom_id(index, ci);
305  // the cluster index might not be in the selection
306  if(dcdci>=0)
307  con[dcdci] += coor[i];
308  }
309  for ( i = 0; i < n; ++i ) {
310  int molIndex = molecule->get_atom_index_from_dcd_selection(index, i);
311  if ( ! wrapAll && ! molecule->is_water(molIndex) ) continue;
312  int ci = molecule->get_cluster(molIndex);
313  int dcdci = molecule->get_dcd_selection_index_from_atom_id(index, ci);
314  if ( ci == molIndex ) {
315  Vector coni = con[i] / molecule->get_clusterSize(molIndex);
316  Vector trans = ( wrapNearest ?
317  lattice.wrap_nearest_delta(coni) : lattice.wrap_delta(coni) );
318  con[i] = trans;
319  }
320  if(dcdci>0)
321  coor[i] = coor[i] + con[dcdci];
322  else
323  coor[i] = coor[i];
324  }
325  delete [] con;
326 #endif
327 }
328 
329  void wrap_coor_dcd_selection(Vector *coor, Lattice &lattice, double *done, int index) {
330  wrap_coor_int_dcd_selection(coor,lattice,done, index);
331 };
332 
333  void wrap_coor_dcd_selection(FloatVector *coor, Lattice &lattice, float *done, int index) {
334  wrap_coor_int_dcd_selection(coor,lattice,done, index);
335 };
336 
337 void wrap_coor(Vector *coor, Lattice &lattice, double *done) {
338  wrap_coor_int(coor,lattice,done);
339 };
340 
341 void wrap_coor(FloatVector *coor, Lattice &lattice, float *done) {
342  wrap_coor_int(coor,lattice,done);
343 };
344 
345 void Output::coordinate(int timestep, int n, Vector *coor, FloatVector *fcoor,
346  Lattice &lattice)
347 {
348  SimParameters *simParams = Node::Object()->simParameters;
349  Molecule *molecule = Node::Object()->molecule;
350  double coor_wrapped = 0;
351  float fcoor_wrapped = 0;
352 
353  if ( timestep >= 0 ) {
354 
355  // Output a DCD trajectory
356  if ( simParams->dcdFrequency &&
357  ((timestep % simParams->dcdFrequency) == 0) )
358  {
359  wrap_coor(fcoor,lattice,&fcoor_wrapped);
360  output_dcdfile(timestep, n, fcoor,
361  simParams->dcdUnitCell ? &lattice : NULL, 0);
362  }
363 
364  // Output any DCD selection trajectories that this timestep mods
365  if ( simParams->dcdSelectionOn)
366  {
367  bool needCleanup = false;
368 
369  for(int index=0; index<16;++index)
370  {
371  int frequency=molecule->dcdSelectionParams[index].frequency;
372  uint16 tag=molecule->dcdSelectionParams[index].tag;
373  int size=molecule->dcdSelectionParams[index].size;
374  if(frequency>0 && (timestep % frequency) == 0)
375  {
376  FloatVector *dcdSelectionFcoor = fcoor;
377  int dcdSelectionNpoints = n;
378  if(dcdSelectionNpoints != molecule->get_dcd_selection_size(index))
379  {
380  // on intervals when some other output is also happening, we have
381  // to filter out a copy of the selection from the full
382  // collected output then wrap and output accordingly
383  dcdSelectionNpoints = size;
384  dcdSelectionFcoor = new FloatVector[dcdSelectionNpoints];
385  needCleanup = true;
386  for(int i=0; i<n ;i++)
387  {
388  int offset = molecule->get_dcd_selection_index_from_atom_id(index, i);
389  if(offset >= 0)
390  dcdSelectionFcoor[offset] = fcoor[i];
391  }
392  }
393  else
394  {
395  dcdSelectionFcoor = fcoor;
396  }
397  wrap_coor_dcd_selection(dcdSelectionFcoor, lattice, &fcoor_wrapped, index);
398  output_dcdfile(timestep, dcdSelectionNpoints, dcdSelectionFcoor,
399  simParams->dcdUnitCell ? &lattice : NULL, tag);
400  if(needCleanup)
401  delete [] dcdSelectionFcoor;
402  }
403 
404  }
405  }
406 
407  // Output a restart file
408  if ( simParams->restartFrequency &&
409  ((timestep % simParams->restartFrequency) == 0) )
410  {
411  iout << "WRITING COORDINATES TO RESTART FILE AT STEP "
412  << timestep << "\n" << endi;
413  wrap_coor(coor,lattice,&coor_wrapped);
414  output_restart_coordinates(coor, n, timestep);
415  iout << "FINISHED WRITING RESTART COORDINATES\n" <<endi;
416  fflush(stdout);
417  }
418 
419  // Interactive MD
420  if ( simParams->IMDon &&
421  ( ((timestep % simParams->IMDfreq) == 0) ||
422  (timestep == simParams->firstTimestep) ) )
423  {
424  IMDOutput *imd = NULL;
425 #ifdef NODEGROUP_FORCE_REGISTER
426  if (simParams->CUDASOAintegrate) {
427  CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
428  imd = cpdata.ckLocalBranch()->imd;
429  }
430  else
431 #endif
432  {
433  imd = Node::Object()->imd;
434  }
435  wrap_coor(fcoor,lattice,&fcoor_wrapped);
436  if (imd != NULL) imd->gather_coordinates(timestep, n, fcoor);
437  }
438 
439  }
440 
441  if (timestep == EVAL_MEASURE)
442  {
443 #ifdef NAMD_TCL
444  wrap_coor(coor,lattice,&coor_wrapped);
445  Node::Object()->getScript()->measure(coor);
446 #endif
447  }
448 
449  // Output final coordinates
450  if (timestep == FILE_OUTPUT || timestep == END_OF_RUN)
451  {
452  int realstep = ( timestep == FILE_OUTPUT ?
453  simParams->firstTimestep : simParams->N );
454  iout << "WRITING COORDINATES TO OUTPUT FILE AT STEP "
455  << realstep << "\n" << endi;
456  fflush(stdout);
457  wrap_coor(coor,lattice,&coor_wrapped);
458  output_final_coordinates(coor, n, realstep);
459  }
460 
461  // Close trajectory files
462  if (timestep == END_OF_RUN)
463  {
464  if (simParams->dcdFrequency) output_dcdfile(END_OF_RUN,0,0,
465  simParams->dcdUnitCell
466  ? &lattice : NULL, 0);
467  if ( simParams->dcdSelectionOn)
468  {
469  for(int dcdSelectionIndex=0; dcdSelectionIndex<16;++dcdSelectionIndex)
470  {
471  int dcdSelectionFrequency = molecule->dcdSelectionParams[dcdSelectionIndex].frequency;
472  int dcdSelectionTag = molecule->dcdSelectionParams[dcdSelectionIndex].tag;
473  if(dcdSelectionFrequency)
474  output_dcdfile(END_OF_RUN,0,0,
475  simParams->dcdUnitCell
476  ? &lattice : NULL, dcdSelectionTag);
477  }
478  }
479  }
480 
481 }
482 /* END OF FUNCTION coordinate */
483 
484 /************************************************************************/
485 /* */
486 /* FUNCTION velocity */
487 /* */
488 /* INPUTS: */
489 /* timestep - Timestep velocities were accumulated for */
490 /* n - This is the number of velocities accumulated. */
491 /* vel - Array of Vectors containing the velocities */
492 /* */
493 /* This function receives the velocities accumulated for a given */
494 /* timestep from the Collect object and calls the appropriate output */
495 /* functions. ALL routines used to output velocity information should*/
496 /* be called from here. */
497 /* */
498 /************************************************************************/
499 
500 int Output::velocityNeeded(int timestep)
501 {
502  SimParameters *simParams = Node::Object()->simParameters;
503 
504  if(simParams->benchTimestep) return 0;
505 
506  int velocitiesNeeded = 0;
507 
508  if ( timestep >= 0 ) {
509 
510  // Output a velocity DCD trajectory
511  if ( simParams->velDcdFrequency &&
512  ((timestep % simParams->velDcdFrequency) == 0) )
513  { velocitiesNeeded |= 1; }
514 
515  // Output a restart file
516  if ( simParams->restartFrequency &&
517  ((timestep % simParams->restartFrequency) == 0) )
518  { velocitiesNeeded |= 2; }
519 
520  }
521 
522  // Output final velocities
523  if (timestep == FILE_OUTPUT || timestep == END_OF_RUN)
524  {
525  velocitiesNeeded |= 2;
526  }
527 
528  return velocitiesNeeded;
529 }
530 
531 void Output::velocity(int timestep, int n, Vector *vel)
532 {
533  SimParameters *simParams = Node::Object()->simParameters;
534 
535  if ( timestep >= 0 ) {
536 
537  // Output velocity DCD trajectory
538  if ( simParams->velDcdFrequency &&
539  ((timestep % simParams->velDcdFrequency) == 0) )
540  {
541  output_veldcdfile(timestep, n, vel);
542  }
543 
544  // Output restart file
545  if ( simParams->restartFrequency &&
546  ((timestep % simParams->restartFrequency) == 0) )
547  {
548  iout << "WRITING VELOCITIES TO RESTART FILE AT STEP "
549  << timestep << "\n" << endi;
550  output_restart_velocities(timestep, n, vel);
551  iout << "FINISHED WRITING RESTART VELOCITIES\n" <<endi;
552  fflush(stdout);
553  }
554 
555  }
556 
557  // Output final velocities
558  if (timestep == FILE_OUTPUT || timestep == END_OF_RUN)
559  {
560  int realstep = ( timestep == FILE_OUTPUT ?
561  simParams->firstTimestep : simParams->N );
562  iout << "WRITING VELOCITIES TO OUTPUT FILE AT STEP "
563  << realstep << "\n" << endi;
564  fflush(stdout);
565  output_final_velocities(realstep, n, vel);
566  }
567 
568  // Close trajectory files
569  if (timestep == END_OF_RUN)
570  {
571  if (simParams->velDcdFrequency) output_veldcdfile(END_OF_RUN,0,0);
572  // close force dcd file here since no final force output below
573  if (simParams->forceDcdFrequency) output_forcedcdfile(END_OF_RUN,0,0);
574  }
575 
576 }
577 /* END OF FUNCTION velocity */
578 
579 /************************************************************************/
580 /* */
581 /* FUNCTION force */
582 /* */
583 /* INPUTS: */
584 /* timestep - Timestep forces were accumulated for */
585 /* n - This is the number of forces accumulated. */
586 /* frc - Array of Vectors containing the forces */
587 /* */
588 /* This function receives the forces accumulated for a given */
589 /* timestep from the Collect object and calls the appropriate output */
590 /* functions. ALL routines used to output force information should*/
591 /* be called from here. */
592 /* */
593 /************************************************************************/
594 
595 int Output::forceNeeded(int timestep)
596 {
597  SimParameters *simParams = Node::Object()->simParameters;
598 
599  if(simParams->benchTimestep) return 0;
600 
601  int forcesNeeded = 0;
602 
603  if ( timestep >= 0 ) {
604 
605  // Output a force DCD trajectory
606  if ( simParams->forceDcdFrequency &&
607  ((timestep % simParams->forceDcdFrequency) == 0) )
608  { forcesNeeded |= 1; }
609 
610  }
611 
612  // Output forces
613  if (timestep == FORCE_OUTPUT)
614  {
615  forcesNeeded |= 2;
616  }
617 
618  return forcesNeeded;
619 }
620 
621 void Output::force(int timestep, int n, Vector *frc)
622 {
623  SimParameters *simParams = Node::Object()->simParameters;
624 
625  if ( timestep >= 0 ) {
626 
627  // Output force DCD trajectory
628  if ( simParams->forceDcdFrequency &&
629  ((timestep % simParams->forceDcdFrequency) == 0) )
630  {
631  output_forcedcdfile(timestep, n, frc);
632  }
633 
634  }
635 
636  // Output forces
637  if (timestep == FORCE_OUTPUT)
638  {
639  int realstep = simParams->firstTimestep;
640  iout << "WRITING FORCES TO OUTPUT FILE AT STEP "
641  << realstep << "\n" << endi;
642  fflush(stdout);
643  output_forces(realstep, n, frc);
644  }
645 
646  // Trajectory file closed by velocity() above
647 
648 }
649 /* END OF FUNCTION force */
650 
651 /************************************************************************/
652 /* */
653 /* FUNCTION output_restart_coordinates */
654 /* */
655 /* INPUTS: */
656 /* coor - Array of vectors containing current positions */
657 /* n - Number of coordinates to output */
658 /* timestep - Timestep for which the coordinates are being written */
659 /* */
660 /* This function writes out the current positions of all the atoms */
661 /* in PDB format to the restart file specified by the user in the */
662 /* configuration file. */
663 /* */
664 /************************************************************************/
665 
666 void Output::output_restart_coordinates(Vector *coor, int n, int timestep)
667 
668 {
669  char comment[128]; // Comment for header of PDB file
670  char timestepstr[20];
671 
672  int baselen = strlen(namdMyNode->simParams->restartFilename);
673  char *restart_name = new char[baselen+26];
674  const char *bsuffix = ".old";
675 
676  strcpy(restart_name, namdMyNode->simParams->restartFilename);
677  if ( namdMyNode->simParams->restartSave ) {
678  sprintf(timestepstr,".%d",timestep);
679  strcat(restart_name, timestepstr);
680  bsuffix = ".BAK";
681  }
682  strcat(restart_name, ".coor");
683 
684  NAMD_backup_file(restart_name,bsuffix);
685 
686  // Check to see if we should generate a binary or PDB file
687  if (!namdMyNode->simParams->binaryRestart)
688  {
689  // Generate a PDB restart file
690  sprintf(comment, "RESTART COORDINATES WRITTEN BY NAMD AT TIMESTEP %d", timestep);
691 
692 #ifdef NODEGROUP_FORCE_REGISTER
693  if(namdMyNode->simParams->CUDASOAintegrate){
694  CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
695  PatchData* pData = cpdata.ckLocalBranch();
696  pData->pdbData->set_all_positions(coor);
697  pData->pdbData->write(restart_name, comment);
698  }
699  else
700 #endif
701  {
702  namdMyNode->pdbData->set_all_positions(coor);
703  namdMyNode->pdbData->write(restart_name, comment);
704  }
705  }
706  else
707  {
708  // Generate a binary restart file
709  write_binary_file(restart_name, n, coor);
710  }
711 
712  delete [] restart_name;
713 
714  if ( namdMyNode->simParams->restartSaveDcd ) {
715  if ( ! output_dcdfile(END_OF_RUN, 0, 0, 0, 0) ) { // close old file
716  const char *old_name = namdMyNode->simParams->dcdFilename;
717  int old_len = strlen(old_name);
718  char *new_name = new char[old_len+26];
719  strcpy(new_name, old_name);
720  if ( old_len >= 4 && ! strcmp(new_name+old_len-4,".dcd") ) {
721  old_len -= 4;
722  new_name[old_len] = 0;
723  }
724  sprintf(timestepstr,".%d",timestep);
725  strcat(new_name, timestepstr);
726  strcat(new_name, ".dcd");
727  iout << "RENAMING COORDINATE DCD FILE " << old_name << " TO " << new_name << "\n" << endi;
728  NAMD_backup_file(new_name,".BAK");
729  while ( rename(old_name, new_name) ) {
730  if ( errno == EINTR || errno == EXDEV ) continue;
731  char err_msg[257];
732  sprintf(err_msg, "Unable to rename DCD file %s to %s", old_name, new_name);
733  NAMD_err(err_msg);
734  }
735  delete [] new_name;
736  }
737  }
738 
739 }
740 /* END OF FUNCTION output_restart_coordinates */
741 
742 /************************************************************************/
743 /* */
744 /* FUNCTION output_restart_velocities */
745 /* */
746 /* INPUTS: */
747 /* vel - Array of vectors containing current velocities */
748 /* timestep - Timestep for which the velocities are being written */
749 /* */
750 /* This function writes out the current velocites of all the atoms */
751 /* in PDB format to the restart file specified by the user in the */
752 /* configuration file. */
753 /* */
754 /************************************************************************/
755 
756 void Output::output_restart_velocities(int timestep, int n, Vector *vel)
757 
758 {
759  char comment[128]; // comment for the header of PDB file
760  char timestepstr[20];
761 
762  int baselen = strlen(namdMyNode->simParams->restartFilename);
763  char *restart_name = new char[baselen+26];
764  const char *bsuffix = ".old";
765 
766  strcpy(restart_name, namdMyNode->simParams->restartFilename);
767  if ( namdMyNode->simParams->restartSave ) {
768  sprintf(timestepstr,".%d",timestep);
769  strcat(restart_name, timestepstr);
770  bsuffix = ".BAK";
771  }
772  strcat(restart_name, ".vel");
773 
774  NAMD_backup_file(restart_name,bsuffix);
775 
776  // Check to see if we should write out a PDB or a binary file
777  if (!namdMyNode->simParams->binaryRestart)
778  {
779  // Write the coordinates to a PDB file. Multiple them by 20
780  // first to make the numbers bigger
781  sprintf(comment, "RESTART VELOCITIES WRITTEN BY NAMD AT TIMESTEP %d", timestep);
782 
783  scale_vels(vel, n, PDBVELFACTOR);
784 
785 #ifdef NODEGROUP_FORCE_REGISTER
786  if(namdMyNode->simParams->CUDASOAintegrate){
787  CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
788  PatchData* pData = cpdata.ckLocalBranch();
789  pData->pdbData->set_all_positions(vel);
790  pData->pdbData->write(restart_name, comment);
791  }
792  else
793 #endif
794  {
795  namdMyNode->pdbData->set_all_positions(vel);
796  namdMyNode->pdbData->write(restart_name, comment);
797  }
798  scale_vels(vel, n, PDBVELINVFACTOR);
799  }
800  else
801  {
802  // Write the velocities to a binary file
803  write_binary_file(restart_name, n, vel);
804  }
805 
806  delete [] restart_name;
807 }
808 /* END OF FUNCTION output_restart_velocities */
809 
810 
811 
812 // This is here so it can access Output.h and Node.h
813 void SimParameters::close_dcdfile() {
814 
815  Output *output = Node::Object()->output;
816  if ( ! output ) return;
817 
818  output->output_dcdfile(END_OF_RUN, 0, 0, 0, 0);
819  Molecule *molecule = Node::Object()->molecule;
820  for(int dcdSelectionIndex=0; dcdSelectionIndex<16;++dcdSelectionIndex)
821  {
822  int dcdSelectionFrequency = molecule->dcdSelectionParams[dcdSelectionIndex].frequency;
823  int dcdSelectionTag = molecule->dcdSelectionParams[dcdSelectionIndex].tag;
824  if(dcdSelectionFrequency)
825  output->output_dcdfile(END_OF_RUN,0,0,0,dcdSelectionTag);
826  }
827 }
828 
829 void SimParameters::close_veldcdfile() {
830 
831  Output *output = Node::Object()->output;
832  if ( ! output ) return;
833 
834  output->output_veldcdfile(END_OF_RUN, 0, 0);
835 
836 }
837 
838 void Output::setReplicaDcdIndex(int index) {
839  replicaDcdActive = 1;
840  replicaDcdIndex = index;
841 }
842 
843 void Output::replicaDcdInit(int index, const char *filename) {
844  replicaDcdActive = 1;
845  replicaDcdIndex = index;
846  int msgsize = sizeof(ReplicaDcdInitMsg) + strlen(filename);
847  ReplicaDcdInitMsg *msg = (ReplicaDcdInitMsg *) CmiAlloc(msgsize);
848  msg->srcPart = CmiMyPartition();
849  msg->dcdIndex = replicaDcdIndex;
850  msg->dcdSelectIndex = 16;
851  strcpy(msg->data, filename);
852  sendReplicaDcdInit(abs(replicaDcdIndex) % CmiNumPartitions(), msg, msgsize);
853 }
854 
855 void Output::replicaDcdSelectInit(int index, const char *tag, const char *filename) {
856  replicaDcdActive = 1;
857  replicaDcdIndex = index;
858  int msgsize = sizeof(ReplicaDcdInitMsg) + strlen(filename);
859  ReplicaDcdInitMsg *msg = (ReplicaDcdInitMsg *) CmiAlloc(msgsize);
860  msg->srcPart = CmiMyPartition();
861  msg->dcdIndex = replicaDcdIndex;
862  Molecule *molecule = namdMyNode->molecule;
863  msg->dcdSelectIndex = molecule->find_or_create_dcd_selection_index(tag);
864  strcpy(msg->data, filename);
865  sendReplicaDcdInit(abs(replicaDcdIndex) % CmiNumPartitions(), msg, msgsize);
866 }
867 
868 
869 void Output::recvReplicaDcdInit(ReplicaDcdInitMsg *msg) {
870  replicaDcdFile &f = (msg->dcdSelectIndex<16) ? replicaDcdSelectFiles[std::make_pair(msg->dcdIndex,msg->dcdSelectIndex)]: replicaDcdFiles[msg->dcdIndex];
871  if ( f.fileid ) {
872  iout << "CLOSING REPLICA DCD FILE " << msg->dcdIndex << " " << f.filename.c_str() << "\n" << endi;
873  close_dcd_write(f.fileid);
874  f.fileid = 0;
875  }
876  f.filename = (const char*) msg->data;
877  sendReplicaDcdAck(msg->srcPart, (ReplicaDcdAckMsg*) CmiAlloc(sizeof(ReplicaDcdAckMsg)));
878 }
879 
880 void Output::recvReplicaDcdData(ReplicaDcdDataMsg *msg) {
881  Molecule *molecule = namdMyNode->molecule;
882  if ( msg->dcdSelectIndex>=16 && ! replicaDcdFiles.count(msg->dcdIndex) ||
883  ( ! replicaDcdFiles.count(msg->dcdIndex) )) {
884  char err_msg[257];
885  sprintf(err_msg, "Unknown replicaDcdFile identifier %d\n", msg->dcdIndex);
886  NAMD_die(err_msg);
887  }
888  replicaDcdFile &f = (msg->dcdSelectIndex<16) ? replicaDcdSelectFiles[std::make_pair(msg->dcdIndex,msg->dcdSelectIndex)]: replicaDcdFiles[msg->dcdIndex];
889 
890  if ( ! f.fileid ) {
891  // Open the DCD file
892  iout << "OPENING REPLICA DCD FILE " << msg->dcdIndex << " " << f.filename.c_str() << "\n" << endi;
893 
894  f.fileid=open_dcd_write(f.filename.c_str());
895 
896  if (f.fileid == DCD_FILEEXISTS) {
897  char err_msg[257];
898  sprintf(err_msg, "DCD file %s already exists!!", f.filename.c_str());
899  NAMD_err(err_msg);
900  } else if (f.fileid < 0) {
901  char err_msg[257];
902  if(msg->dcdSelectIndex<16)
903  sprintf(err_msg, "Couldn't open DCD Index %d Select Index %d file %s", msg->dcdIndex, msg->dcdSelectIndex, f.filename.c_str());
904  else
905  sprintf(err_msg, "Couldn't open DCD Index %d file %s", msg->dcdIndex, f.filename.c_str());
906  NAMD_err(err_msg);
907  } else if (! f.fileid) {
908  NAMD_bug("Output::recvReplicaDcdData open_dcd_write returned fileid of zero");
909  }
910 
911  // Write out the header
912  int ret_code = write_dcdheader(f.fileid, f.filename.c_str(),
913  msg->numAtoms, msg->NFILE, msg->NPRIV, msg->NSAVC, msg->NSTEP,
914  msg->DELTA, msg->with_unitcell);
915 
916  if (ret_code<0) {
917  NAMD_err("Writing of DCD header failed!!");
918  }
919  }
920 
921  // Write out the values for this timestep
922  iout << "WRITING TO REPLICA DCD FILE " << msg->dcdIndex << " " << f.filename.c_str() << "\n" << endi;
923  float *msgx = (float*) msg->data;
924  float *msgy = msgx + msg->numAtoms;
925  float *msgz = msgy + msg->numAtoms;
926  int ret_code = write_dcdstep(f.fileid, msg->numAtoms, msgx, msgy, msgz,
927  msg->with_unitcell ? msg->unitcell : 0);
928  if (ret_code < 0) NAMD_err("Writing of DCD step failed!!");
929 
930  sendReplicaDcdAck(msg->srcPart, (ReplicaDcdAckMsg*) CmiAlloc(sizeof(ReplicaDcdAckMsg)));
931 }
932 
933 
934 /************************************************************************/
935 /* */
936 /* FUNCTION output_dcdfile */
937 /* */
938 /* INPUTS: */
939 /* timestep - Current timestep */
940 /* n - Number of atoms in simulation */
941 /* coor - Coordinate vectors for all atoms */
942 /* lattice - periodic cell data; NULL if not to be written */
943 /* */
944 /* This function maintains the interface between the Output object */
945 /* and the dcd writing routines contained in dcdlib. */
946 /* */
947 /************************************************************************/
948 
949 #define RAD2DEG 180.0/3.14159265359
950 
951 int Output::output_dcdfile(int timestep, int n, FloatVector *coor,
952  const Lattice *lattice, int dcdSelectTag)
953 
954 {
955  static Bool first[maxDCD]={true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true}; // Flag indicating first call
956  static int fileid[maxDCD]; // File id for the dcd file
957 
958  static float *x, *y, *z; // Arrays to hold x, y, and z arrays
959  static int n_alloc; // allocated size
960 
961  int i; // Loop counter
962  int ret_code; // Return code from DCD calls
963  SimParameters *simParams = namdMyNode->simParams;
964  Molecule *molecule = namdMyNode->molecule;
965  // If this is the last time we will be writing coordinates,
966  // close the file before exiting
967  char *dcdFilename = simParams->dcdFilename;
968  int dcdSelectIndex=dcdSelectTag-1;
969  if(dcdSelectTag >0 && dcdSelectTag < 16)
970  {// look it up
971  dcdFilename = molecule->dcdSelectionParams[dcdSelectIndex].outFilename;
972  }
973  if ( timestep == END_OF_RUN ) {
974  for ( std::map<int,replicaDcdFile>::iterator it = replicaDcdFiles.begin();
975  it != replicaDcdFiles.end(); ++it ) {
976  replicaDcdFile &f = it->second;
977  if ( f.fileid ) {
978  iout << "CLOSING REPLICA DCD FILE " << it->first << " " << f.filename.c_str() << "\n" << endi;
979  close_dcd_write(f.fileid);
980  f.fileid = 0;
981  }
982  }
983  for ( std::map<std::pair<int,uint16_t>,replicaDcdFile>::iterator it = replicaDcdSelectFiles.begin();
984  it != replicaDcdSelectFiles.end(); ++it ) {
985  replicaDcdFile &f = it->second;
986  if ( f.fileid ) {
987  iout << "CLOSING REPLICA DCD SELECT FILE " << it->first.first << " " << f.filename.c_str() << "\n" << endi;
988  close_dcd_write(f.fileid);
989  f.fileid = 0;
990  }
991  }
992 
993  int rval = 0;
994  if ( ! first[dcdSelectTag] ) {
995  iout << "CLOSING COORDINATE DCD FILE " << dcdFilename << "\n" << endi;
996  close_dcd_write(fileid[dcdSelectTag]);
997  } else {
998  iout << "COORDINATE DCD FILE " << dcdFilename << " WAS NOT CREATED\n" << endi;
999  rval = -1;
1000  }
1001  first[dcdSelectTag] = 1;
1002  fileid[dcdSelectTag] = 0;
1003  return rval;
1004  }
1005 
1006  if ( replicaDcdActive ) {
1007  int msgsize = sizeof(ReplicaDcdDataMsg) + 3*n*sizeof(float);
1008  ReplicaDcdDataMsg *msg = (ReplicaDcdDataMsg *) CmiAlloc(msgsize);
1009  float *msgx = (float*) msg->data;
1010  float *msgy = msgx + n;
1011  float *msgz = msgy + n;
1012  for (i=0; i<n; i++) { msgx[i] = coor[i].x; }
1013  for (i=0; i<n; i++) { msgy[i] = coor[i].y; }
1014  for (i=0; i<n; i++) { msgz[i] = coor[i].z; }
1015  msg->numAtoms = n;
1016  lattice_to_unitcell(lattice,msg->unitcell);
1017  msg->with_unitcell = lattice ? 1 : 0;
1018  if(dcdSelectTag==0)
1019  msg->NSAVC = simParams->dcdFrequency;
1020  else
1021  msg->NSAVC = molecule->dcdSelectionParams[dcdSelectIndex].frequency;
1022  msg->NPRIV = timestep;
1023  msg->NSTEP = msg->NPRIV - msg->NSAVC;
1024  msg->NFILE = 0;
1025  msg->DELTA = simParams->dt/TIMEFACTOR;
1026  msg->srcPart = CmiMyPartition();
1027  msg->dcdIndex = replicaDcdIndex;
1028  msg->dcdSelectIndex = dcdSelectIndex;
1029  sendReplicaDcdData(abs(replicaDcdIndex) % CmiNumPartitions(), msg, msgsize);
1030  return 0;
1031  }
1032  if (first[dcdSelectTag])
1033  {
1034  // Allocate x, y, and z arrays since the DCD file routines
1035  // need them passed as three independant arrays to be
1036  // efficient
1037  if ( n > n_alloc ) {
1038  delete [] x; x = new float[3*n];
1039  y = x + n;
1040  z = x + 2*n;
1041  n_alloc = n;
1042  }
1043 
1044  // Open the DCD file
1045  iout << "OPENING COORDINATE DCD FILE\n" << endi;
1046  fileid[dcdSelectTag] = open_dcd_write(dcdFilename);
1047  if (fileid[dcdSelectTag] == DCD_FILEEXISTS)
1048  {
1049  char err_msg[257];
1050 
1051  sprintf(err_msg, "DCD file %s already exists!!",
1052  dcdFilename);
1053 
1054  NAMD_err(err_msg);
1055  }
1056  else if (fileid[dcdSelectTag] < 0)
1057  {
1058  char err_msg[257];
1059 
1060  sprintf(err_msg, "Couldn't open DCD file %s",
1061  dcdFilename);
1062 
1063  NAMD_err(err_msg);
1064  }
1065 
1066  int NSAVC, NFILE, NPRIV, NSTEP;
1067  if(dcdSelectTag==0)
1068  NSAVC = simParams->dcdFrequency;
1069  else
1070  NSAVC = molecule->dcdSelectionParams[dcdSelectIndex].frequency;
1071  NPRIV = timestep;
1072  NSTEP = NPRIV - NSAVC;
1073  NFILE = 0;
1074 
1075  // Write out the header
1076  ret_code = write_dcdheader(fileid[dcdSelectTag],
1077  dcdFilename,
1078  n, NFILE, NPRIV, NSAVC, NSTEP,
1079  simParams->dt/TIMEFACTOR, lattice != NULL);
1080 
1081  if (ret_code<0)
1082  {
1083  NAMD_err("Writing of DCD header failed!!");
1084  }
1085 
1086  first[dcdSelectTag] = FALSE;
1087  }
1088 
1089  // Copy the coordinates for output
1090  for (i=0; i<n; i++)
1091  {
1092  x[i] = coor[i].x;
1093  y[i] = coor[i].y;
1094  z[i] = coor[i].z;
1095  }
1096 
1097  // Write out the values for this timestep
1098  iout << "WRITING COORDINATES TO DCD FILE " << dcdFilename << " AT STEP "
1099  << timestep << "\n" << endi;
1100  fflush(stdout);
1101  if (lattice) {
1102  double unitcell[6];
1103  lattice_to_unitcell(lattice,unitcell);
1104  ret_code = write_dcdstep(fileid[dcdSelectTag], n, x, y, z, unitcell);
1105  } else {
1106  ret_code = write_dcdstep(fileid[dcdSelectTag], n, x, y, z, NULL);
1107  }
1108  if (ret_code < 0)
1109  {
1110  NAMD_err("Writing of DCD step failed!!");
1111  }
1112 
1113  return 0;
1114 }
1115 /* END OF FUNCTION output_dcdfile */
1116 
1117 /************************************************************************/
1118 /* */
1119 /* FUNCTION output_final_coordinates */
1120 /* */
1121 /* INPUTS: */
1122 /* coor - Array of vectors containing final coordinates */
1123 /* n - Number of coordinates to output */
1124 /* timestep - Timestep that coordinates are being written in */
1125 /* */
1126 /* This function writes out the final coordinates for the */
1127 /* simulation in PDB format to the file specified in the config */
1128 /* file. */
1129 /* */
1130 /************************************************************************/
1131 
1132 void Output::output_final_coordinates(Vector *coor, int n, int timestep)
1133 
1134 {
1135  char output_name[140]; // Output filename
1136  char comment[128]; // comment for PDB header
1137 
1138  // Built the output filename
1139  strcpy(output_name, namdMyNode->simParams->outputFilename);
1140  strcat(output_name, ".coor");
1141 
1142  NAMD_backup_file(output_name);
1143 
1144  // Check to see if we should write out a binary file or a
1145  // PDB file
1146  if (!namdMyNode->simParams->binaryOutput)
1147  {
1148  sprintf(comment, "FINAL COORDINATES WRITTEN BY NAMD AT TIMESTEP %d", timestep);
1149 
1150 #ifdef NODEGROUP_FORCE_REGISTER
1151  if(namdMyNode->simParams->CUDASOAintegrate){
1152  CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
1153  PatchData* pData = cpdata.ckLocalBranch();
1154  pData->pdbData->set_all_positions(coor);
1155  pData->pdbData->write(output_name, comment);
1156  }
1157  else
1158 #endif
1159  {
1160  namdMyNode->pdbData->set_all_positions(coor);
1161  namdMyNode->pdbData->write(output_name, comment);
1162  }
1163  }
1164  else
1165  {
1166  // Write the velocities to a binary file
1167  write_binary_file(output_name, n, coor);
1168  }
1169 }
1170 /* END OF FUNCTION output_final_coordinates */
1171 
1172 /************************************************************************/
1173 /* */
1174 /* FUNCTION output_final_velocities */
1175 /* */
1176 /* INPUTS: */
1177 /* vel - Array of vectors containing final velocities */
1178 /* timestep - Timestep that vleocities are being written in */
1179 /* */
1180 /* This function writes out the final vleocities for the */
1181 /* simulation in PDB format to the file specified in the config */
1182 /* file. */
1183 /* */
1184 /************************************************************************/
1185 
1186 void Output::output_final_velocities(int timestep, int n, Vector *vel)
1187 
1188 {
1189  char output_name[140]; // Output filename
1190  char comment[128]; // Comment for PDB header
1191 
1192  // Build the output filename
1193  strcpy(output_name, namdMyNode->simParams->outputFilename);
1194  strcat(output_name, ".vel");
1195 
1196  NAMD_backup_file(output_name);
1197 
1198  // Check to see if we should write a PDB or binary file
1199  if (!(namdMyNode->simParams->binaryOutput))
1200  {
1201  // Write the final velocities to a PDB file
1202  sprintf(comment, "FINAL VELOCITIES WRITTEN BY NAMD AT TIMESTEP %d", timestep);
1203 
1204  scale_vels(vel, n, PDBVELFACTOR);
1205 #ifdef NODEGROUP_FORCE_REGISTER
1206  if(namdMyNode->simParams->CUDASOAintegrate){
1207  CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
1208  PatchData* pData = cpdata.ckLocalBranch();
1209  pData->pdbData->set_all_positions(vel);
1210  pData->pdbData->write(output_name, comment);
1211  }
1212  else
1213 #endif
1214  {
1215  namdMyNode->pdbData->set_all_positions(vel);
1216  namdMyNode->pdbData->write(output_name, comment);
1217  }
1218  scale_vels(vel, n, PDBVELINVFACTOR);
1219  }
1220  else
1221  {
1222  // Write the coordinates to a binary file
1223  write_binary_file(output_name, n, vel);
1224  }
1225 
1226 }
1227 /* END OF FUNCTION output_final_velocities */
1228 
1229 /************************************************************************/
1230 /* */
1231 /* FUNCTION output_veldcdfile */
1232 /* */
1233 /* INPUTS: */
1234 /* timestep - Current timestep */
1235 /* n - Number of atoms in simulation */
1236 /* coor - velocity vectors for all atoms */
1237 /* */
1238 /* This function maintains the interface between the Output object */
1239 /* and the dcd writing routines contained in dcdlib. This fucntion */
1240 /* writes out the velocity vectors in DCD format. */
1241 /* */
1242 /************************************************************************/
1243 
1244 void Output::output_veldcdfile(int timestep, int n, Vector *vel)
1245 
1246 {
1247  static Bool first=TRUE; // Flag indicating first call
1248  static int fileid; // File id for the dcd file
1249  static float *x, *y, *z; // Arrays to hold x, y, and z arrays
1250  static int n_alloc; // allocated size
1251  int i; // Loop counter
1252  int ret_code; // Return code from DCD calls
1253  SimParameters *simParams = Node::Object()->simParameters;
1254 
1255  // If this is the last time we will be writing coordinates,
1256  // close the file before exiting
1257  if ( timestep == END_OF_RUN ) {
1258  if ( ! first ) {
1259  iout << "CLOSING VELOCITY DCD FILE\n" << endi;
1260  close_dcd_write(fileid);
1261  } else {
1262  iout << "VELOCITY DCD FILE WAS NOT CREATED\n" << endi;
1263  }
1264  first = TRUE;
1265  fileid = 0;
1266  return;
1267  }
1268 
1269  if (first)
1270  {
1271  // Allocate x, y, and z arrays since the DCD file routines
1272  // need them passed as three independant arrays to be
1273  // efficient
1274  if ( n > n_alloc ) {
1275  delete [] x; x = new float[3*n];
1276  y = x + n;
1277  z = x + 2*n;
1278  n_alloc = n;
1279  }
1280 
1281  // Open the DCD file
1282  iout << "OPENING VELOCITY DCD FILE\n" << endi;
1283 
1284  fileid=open_dcd_write(namdMyNode->simParams->velDcdFilename);
1285 
1286  if (fileid == DCD_FILEEXISTS)
1287  {
1288  char err_msg[257];
1289 
1290  sprintf(err_msg, "Velocity DCD file %s already exists!!",
1291  namdMyNode->simParams->velDcdFilename);
1292 
1293  NAMD_err(err_msg);
1294  }
1295  else if (fileid < 0)
1296  {
1297  char err_msg[257];
1298 
1299  sprintf(err_msg, "Couldn't open velocity DCD file %s",
1300  namdMyNode->simParams->velDcdFilename);
1301 
1302  NAMD_err(err_msg);
1303  }
1304 
1305  int NSAVC, NFILE, NPRIV, NSTEP;
1306  NSAVC = simParams->velDcdFrequency;
1307  NPRIV = timestep;
1308  NSTEP = NPRIV - NSAVC;
1309  NFILE = 0;
1310 
1311  // Write out the header
1312  const int with_unitcell = 0;
1313  ret_code = write_dcdheader(fileid,
1314  simParams->velDcdFilename,
1315  n, NFILE, NPRIV, NSAVC, NSTEP,
1316  simParams->dt/TIMEFACTOR, with_unitcell);
1317 
1318 
1319  if (ret_code<0)
1320  {
1321  NAMD_err("Writing of velocity DCD header failed!!");
1322  }
1323 
1324  first = FALSE;
1325  }
1326 
1327  // Copy the coordinates for output
1328  for (i=0; i<n; i++)
1329  {
1330  x[i] = vel[i].x;
1331  y[i] = vel[i].y;
1332  z[i] = vel[i].z;
1333  }
1334 
1335  // Write out the values for this timestep
1336  iout << "WRITING VELOCITIES TO DCD FILE AT STEP "
1337  << timestep << "\n" << endi;
1338  fflush(stdout);
1339  ret_code = write_dcdstep(fileid, n, x, y, z, NULL);
1340 
1341  if (ret_code < 0)
1342  {
1343  NAMD_err("Writing of velocity DCD step failed!!");
1344  }
1345 
1346 }
1347 /* END OF FUNCTION output_veldcdfile */
1348 
1349 /************************************************************************/
1350 /* */
1351 /* FUNCTION output_forces */
1352 /* */
1353 /* INPUTS: */
1354 /* frc - Array of vectors containing final forces */
1355 /* timestep - Timestep that vleocities are being written in */
1356 /* */
1357 /* This function writes out the final vleocities for the */
1358 /* simulation in PDB format to the file specified in the config */
1359 /* file. */
1360 /* */
1361 /************************************************************************/
1362 
1363 void Output::output_forces(int timestep, int n, Vector *frc)
1364 
1365 {
1366  char output_name[140]; // Output filename
1367  char comment[128]; // Comment for PDB header
1368 
1369  // Build the output filename
1370  strcpy(output_name, namdMyNode->simParams->outputFilename);
1371  strcat(output_name, ".force");
1372 
1373  NAMD_backup_file(output_name);
1374 
1375  // Check to see if we should write a PDB or binary file
1376  if (!(namdMyNode->simParams->binaryOutput))
1377  {
1378  // Write the forces to a PDB file
1379  sprintf(comment, "FORCES WRITTEN BY NAMD AT TIMESTEP %d", timestep);
1380 #ifdef NODEGROUP_FORCE_REGISTER
1381  if(namdMyNode->simParams->CUDASOAintegrate){
1382  CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
1383  PatchData* pData = cpdata.ckLocalBranch();
1384  pData->pdbData->set_all_positions(frc);
1385  pData->pdbData->write(output_name, comment);
1386  }
1387  else
1388 #endif
1389  {
1390  namdMyNode->pdbData->set_all_positions(frc);
1391  namdMyNode->pdbData->write(output_name, comment);
1392  }
1393  }
1394  else
1395  {
1396  // Write the coordinates to a binary file
1397  write_binary_file(output_name, n, frc);
1398  }
1399 
1400 }
1401 /* END OF FUNCTION output_forces */
1402 
1403 /************************************************************************/
1404 /* */
1405 /* FUNCTION output_forcedcdfile */
1406 /* */
1407 /* INPUTS: */
1408 /* timestep - Current timestep */
1409 /* n - Number of atoms in simulation */
1410 /* frc - force vectors for all atoms */
1411 /* */
1412 /* This function maintains the interface between the Output object */
1413 /* and the dcd writing routines contained in dcdlib. This fucntion */
1414 /* writes out the force vectors in DCD format. */
1415 /* */
1416 /************************************************************************/
1417 
1418 void Output::output_forcedcdfile(int timestep, int n, Vector *frc)
1419 
1420 {
1421  static Bool first=TRUE; // Flag indicating first call
1422  static int fileid; // File id for the dcd file
1423  static float *x, *y, *z; // Arrays to hold x, y, and z arrays
1424  static int n_alloc; // allocated size
1425  int i; // Loop counter
1426  int ret_code; // Return code from DCD calls
1427  SimParameters *simParams = Node::Object()->simParameters;
1428 
1429  // If this is the last time we will be writing coordinates,
1430  // close the file before exiting
1431  if ( timestep == END_OF_RUN ) {
1432  if ( ! first ) {
1433  iout << "CLOSING FORCE DCD FILE\n" << endi;
1434  close_dcd_write(fileid);
1435  } else {
1436  iout << "FORCE DCD FILE WAS NOT CREATED\n" << endi;
1437  }
1438  return;
1439  }
1440 
1441  if (first)
1442  {
1443  // Allocate x, y, and z arrays since the DCD file routines
1444  // need them passed as three independant arrays to be
1445  // efficient
1446  if ( n > n_alloc ) {
1447  delete [] x; x = new float[3*n];
1448  y = x + n;
1449  z = x + 2*n;
1450  n_alloc = n;
1451  }
1452 
1453  // Open the DCD file
1454  iout << "OPENING FORCE DCD FILE\n" << endi;
1455 
1456  fileid=open_dcd_write(namdMyNode->simParams->forceDcdFilename);
1457 
1458  if (fileid == DCD_FILEEXISTS)
1459  {
1460  char err_msg[257];
1461 
1462  sprintf(err_msg, "Force DCD file %s already exists!!",
1463  namdMyNode->simParams->forceDcdFilename);
1464 
1465  NAMD_err(err_msg);
1466  }
1467  else if (fileid < 0)
1468  {
1469  char err_msg[257];
1470 
1471  sprintf(err_msg, "Couldn't open force DCD file %s",
1472  namdMyNode->simParams->forceDcdFilename);
1473 
1474  NAMD_err(err_msg);
1475  }
1476 
1477  int NSAVC, NFILE, NPRIV, NSTEP;
1478  NSAVC = simParams->forceDcdFrequency;
1479  NPRIV = timestep;
1480  NSTEP = NPRIV - NSAVC;
1481  NFILE = 0;
1482 
1483  // Write out the header
1484  const int with_unitcell = 0;
1485  ret_code = write_dcdheader(fileid,
1486  simParams->forceDcdFilename,
1487  n, NFILE, NPRIV, NSAVC, NSTEP,
1488  simParams->dt/TIMEFACTOR, with_unitcell);
1489 
1490 
1491  if (ret_code<0)
1492  {
1493  NAMD_err("Writing of force DCD header failed!!");
1494  }
1495 
1496  first = FALSE;
1497  }
1498 
1499  // Copy the coordinates for output
1500  for (i=0; i<n; i++)
1501  {
1502  x[i] = frc[i].x;
1503  y[i] = frc[i].y;
1504  z[i] = frc[i].z;
1505  }
1506 
1507  // Write out the values for this timestep
1508  iout << "WRITING FORCES TO DCD FILE AT STEP "
1509  << timestep << "\n" << endi;
1510  fflush(stdout);
1511  ret_code = write_dcdstep(fileid, n, x, y, z, NULL);
1512 
1513  if (ret_code < 0)
1514  {
1515  NAMD_err("Writing of force DCD step failed!!");
1516  }
1517 
1518 }
1519 /* END OF FUNCTION output_forcedcdfile */
1520 
1521 /************************************************************************/
1522 /* */
1523 /* FUNCTION write_binary_file */
1524 /* */
1525 /* INPUTS: */
1526 /* fname - file name to write velocities to */
1527 /* n - Number of atoms in system */
1528 /* vels - Array of vectors */
1529 /* */
1530 /* This function writes out vectors in binary format to */
1531 /* the specified file. */
1532 /* */
1533 /************************************************************************/
1534 
1535 void Output::write_binary_file(char *fname, int n, Vector *vecs)
1536 
1537 {
1538  char errmsg[256];
1539  int fd; // File descriptor
1540  int32 n32 = n;
1541 
1542  fd = NAMD_open(fname);
1543 
1544  sprintf(errmsg, "Error on write to binary file %s", fname);
1545 
1546  // Write out the number of atoms and the vectors
1547  NAMD_write(fd, (char *) &n32, sizeof(int32), errmsg);
1548  NAMD_write(fd, (char *) vecs, sizeof(Vector)*n, errmsg);
1549 
1550  NAMD_close(fd, fname);
1551 }
1552 /* END OF FUNCTION write_binary_file */
1553 
1554 /************************************************************************/
1555 /* */
1556 /* FUNCTION scale_vels */
1557 /* */
1558 /* INPUTS: */
1559 /* v - Array of velocity vectors */
1560 /* n - Number of atoms in system */
1561 /* fact - Scaling factor */
1562 /* */
1563 /* This function scales all the vectors passed in by a constant */
1564 /* factor. This is used before writing out velocity vectors that */
1565 /* need to be resized. */
1566 /* */
1567 /************************************************************************/
1568 
1569 void Output::scale_vels(Vector *v, int n, Real fact)
1570 
1571 {
1572  int i;
1573 
1574  for (i=0; i<n; i++)
1575  {
1576  v[i].x *= fact;
1577  v[i].y *= fact;
1578  v[i].z *= fact;
1579  }
1580 }
1581 /* END OF FUNCTION scale_vels */
1582 
1583 
1584 #ifdef MEM_OPT_VERSION
1585 void ParOutput::velocityMaster(int timestep, int n){
1587  SimParameters *simParams = Node::Object()->simParameters;
1588 
1589  if ( timestep >= 0 ) {
1590 
1591  // Output velocity DCD trajectory
1592  if ( simParams->velDcdFrequency &&
1593  ((timestep % simParams->velDcdFrequency) == 0) )
1594  {
1595  output_veldcdfile_master(timestep, n);
1596  }
1597 
1598  // Output restart file
1599  if ( simParams->restartFrequency &&
1600  ((timestep % simParams->restartFrequency) == 0) )
1601  {
1602  iout << "WRITING VELOCITIES TO RESTART FILE AT STEP "
1603  << timestep << "\n" << endi;
1604  output_restart_velocities_master(timestep, n);
1605  iout << "FINISHED WRITING RESTART VELOCITIES\n" <<endi;
1606  fflush(stdout);
1607  }
1608 
1609  }
1610 
1611  // Output final velocities
1612  if (timestep == FILE_OUTPUT || timestep == END_OF_RUN)
1613  {
1614  int realstep = ( timestep == FILE_OUTPUT ?
1615  simParams->firstTimestep : simParams->N );
1616  iout << "WRITING VELOCITIES TO OUTPUT FILE AT STEP "
1617  << realstep << "\n" << endi;
1618  fflush(stdout);
1619  output_final_velocities_master(n);
1620  }
1621 
1622  // Close trajectory files
1623  if (timestep == END_OF_RUN)
1624  {
1625  if (simParams->velDcdFrequency) output_veldcdfile_master(END_OF_RUN,0);
1626  if (simParams->forceDcdFrequency) output_forcedcdfile_master(END_OF_RUN,0);
1627  }
1628 
1629 }
1630 
1631 //output atoms' velocities from id fID to tID.
1632 void ParOutput::velocitySlave(int timestep, int fID, int tID, Vector *vecs){
1633  SimParameters *simParams = Node::Object()->simParameters;
1634 
1635  if ( timestep >= 0 ) {
1636 
1637  // Output velocity DCD trajectory
1638  if ( simParams->velDcdFrequency &&
1639  ((timestep % simParams->velDcdFrequency) == 0) )
1640  {
1641  output_veldcdfile_slave(timestep, fID, tID, vecs);
1642  }
1643 
1644  // Output restart file
1645  if ( simParams->restartFrequency &&
1646  ((timestep % simParams->restartFrequency) == 0) )
1647  {
1648  int64 offset = sizeof(int)+sizeof(Vector)*((int64)fID);
1649  output_restart_velocities_slave(timestep, fID, tID, vecs, offset);
1650  }
1651 
1652  }
1653 
1654  // Output final velocities
1655  if (timestep == FILE_OUTPUT || timestep == END_OF_RUN)
1656  {
1657  int64 offset = sizeof(int)+sizeof(Vector)*((int64)fID);
1658  output_final_velocities_slave(fID, tID, vecs, offset);
1659  }
1660 
1661  // Close trajectory files
1662  if (timestep == END_OF_RUN)
1663  {
1664  if (simParams->velDcdFrequency) output_veldcdfile_slave(END_OF_RUN, 0, 0, NULL);
1665  if (simParams->forceDcdFrequency) output_forcedcdfile_slave(END_OF_RUN, 0, 0, NULL);
1666  }
1667 }
1668 
1669 void ParOutput::output_veldcdfile_master(int timestep, int n){
1670  int ret_code; // Return code from DCD calls
1671  SimParameters *simParams = Node::Object()->simParameters;
1672 
1673  // If this is the last time we will be writing coordinates,
1674  // close the file before exiting
1675  if ( timestep == END_OF_RUN ) {
1676  if ( ! veldcdFirst ) {
1677  iout << "CLOSING VELOCITY DCD FILE\n" << endi;
1678  close_dcd_write(veldcdFileID);
1679  } else {
1680  iout << "VELOCITY DCD FILE WAS NOT CREATED\n" << endi;
1681  }
1682  return;
1683  }
1684 
1685  if (veldcdFirst)
1686  {
1687  // Open the DCD file
1688  iout << "OPENING VELOCITY DCD FILE\n" << endi;
1689 
1690 #ifndef OUTPUT_SINGLE_FILE
1691 #error OUTPUT_SINGLE_FILE not defined!
1692 #endif
1693 
1694  #if OUTPUT_SINGLE_FILE
1695  char *veldcdFilename = simParams->velDcdFilename;
1696  #else
1697  char *veldcdFilename = buildFileName(veldcdType);
1698  #endif
1699 
1700  veldcdFileID=open_dcd_write(veldcdFilename);
1701 
1702  if (veldcdFileID == DCD_FILEEXISTS)
1703  {
1704  char err_msg[257];
1705  sprintf(err_msg, "Velocity DCD file %s already exists!",veldcdFilename);
1706  NAMD_err(err_msg);
1707  }
1708  else if (veldcdFileID < 0)
1709  {
1710  char err_msg[257];
1711  sprintf(err_msg, "Couldn't open velocity DCD file %s",veldcdFilename);
1712  NAMD_err(err_msg);
1713  }
1714 
1715  #if !OUTPUT_SINGLE_FILE
1716  // Write out extra fields as MAGIC number etc.
1717  int32 tmpInt = OUTPUT_MAGIC_NUMBER;
1718  float tmpFlt = OUTPUT_FILE_VERSION;
1719  NAMD_write(veldcdFileID, (char *) &tmpInt, sizeof(int32));
1720  NAMD_write(veldcdFileID, (char *) &tmpFlt, sizeof(float));
1721  tmpInt = simParams->numoutputprocs;
1722  NAMD_write(veldcdFileID, (char *) &tmpInt, sizeof(int32));
1723  #endif
1724 
1725  int NSAVC, NFILE, NPRIV, NSTEP;
1726  NSAVC = simParams->velDcdFrequency;
1727  NPRIV = timestep;
1728  NSTEP = NPRIV - NSAVC;
1729  NFILE = 0;
1730 
1731  // Write out the header
1732  const int with_unitcell = 0;
1733  ret_code = write_dcdheader(veldcdFileID,
1734  veldcdFilename,
1735  n, NFILE, NPRIV, NSAVC, NSTEP,
1736  simParams->dt/TIMEFACTOR, with_unitcell);
1737 
1738 
1739  if (ret_code<0)
1740  {
1741  NAMD_err("Writing of velocity DCD header failed!!");
1742  }
1743 
1744  #if !OUTPUT_SINGLE_FILE
1745  //In this case, the file name is dynamically allocated
1746  delete [] veldcdFilename;
1747  #endif
1748 
1749  veldcdFirst = FALSE;
1750  }
1751 
1752  // Write out the values for this timestep
1753  iout << "WRITING VELOCITIES TO DCD FILE AT STEP "
1754  << timestep << "\n" << endi;
1755  fflush(stdout);
1756 
1757  //In the case of writing to multiple files, only the header
1758  //of the dcd file needs to be updated. Note that the format of
1759  //the new dcd file has also changed! -Chao Mei
1760 
1761 #if OUTPUT_SINGLE_FILE
1762  //write X,Y,Z headers
1763  int totalAtoms = namdMyNode->molecule->numAtoms;
1764  write_dcdstep_par_XYZUnits(veldcdFileID, totalAtoms);
1765 #endif
1766 
1767  //update the header
1768  update_dcdstep_par_header(veldcdFileID);
1769 
1770 }
1771 
1772 void ParOutput::output_veldcdfile_slave(int timestep, int fID, int tID, Vector *vecs){
1773  int ret_code; // Return code from DCD calls
1774  SimParameters *simParams = Node::Object()->simParameters;
1775 
1776  // If this is the last time we will be writing coordinates,
1777  // close the file before exiting
1778  if ( timestep == END_OF_RUN ) {
1779  if ( ! veldcdFirst ) {
1780  close_dcd_write(veldcdFileID);
1781  }
1782 #if OUTPUT_SINGLE_FILE
1783  delete [] veldcdX;
1784  delete [] veldcdY;
1785  delete [] veldcdZ;
1786 #endif
1787  return;
1788  }
1789 
1790  int parN = tID-fID+1;
1791 
1792  if (veldcdFirst)
1793  {
1794 
1795  #if OUTPUT_SINGLE_FILE
1796  char *veldcdFilename = namdMyNode->simParams->velDcdFilename;
1797  #else
1798  char *veldcdFilename = buildFileName(veldcdType);
1799  #endif
1800  veldcdFileID=open_dcd_write_par_slave(veldcdFilename);
1801  if(veldcdFileID < 0)
1802  {
1803  char err_msg[257];
1804  sprintf(err_msg, "Couldn't open velocity DCD file %s",veldcdFilename);
1805  NAMD_err(err_msg);
1806  }
1807  #if OUTPUT_SINGLE_FILE
1808  //If outputting to a single file, dcd files conforms to the old format
1809  //as data are organized as 3 seperate arrays of X,Y,Z, while in the new
1810  //format used in outputing multiple files, the data are organized as an
1811  //array of XYZs.
1812  veldcdX = new float[parN];
1813  veldcdY = new float[parN];
1814  veldcdZ = new float[parN];
1815  //seek to beginning of X,Y,Z sections which means skipping header.
1816  //Cell data is not needed because this is velocity trajectory
1817  int skipbytes = get_dcdheader_size();
1818  seek_dcdfile(veldcdFileID, skipbytes, SEEK_SET);
1819  #endif
1820 
1821  #if !OUTPUT_SINGLE_FILE
1822  delete [] veldcdFilename;
1823  // Write out extra fields as MAGIC number etc.
1824  int32 tmpInt = OUTPUT_MAGIC_NUMBER;
1825  float tmpFlt = OUTPUT_FILE_VERSION;
1826  NAMD_write(veldcdFileID, (char *) &tmpInt, sizeof(int32));
1827  NAMD_write(veldcdFileID, (char *) &tmpFlt, sizeof(float));
1828  NAMD_write(veldcdFileID, (char *) &outputID, sizeof(int));
1829  NAMD_write(veldcdFileID, (char *) &fID, sizeof(int));
1830  NAMD_write(veldcdFileID, (char *) &tID, sizeof(int));
1831  #endif
1832 
1833  veldcdFirst = FALSE;
1834  }
1835 
1836 #if OUTPUT_SINGLE_FILE
1837  //The following seek will set the stream position to the
1838  //beginning of the place where a new timestep output should
1839  //be performed.
1840  CmiAssert(sizeof(off_t)==8);
1841  int totalAtoms = namdMyNode->molecule->numAtoms;
1842 
1843  for(int i=0; i<parN; i++){
1844  veldcdX[i] = vecs[i].x;
1845  veldcdY[i] = vecs[i].y;
1846  veldcdZ[i] = vecs[i].z;
1847  }
1848 
1849  write_dcdstep_par_slave(veldcdFileID, fID, tID, totalAtoms, veldcdX, veldcdY, veldcdZ);
1850 
1851  //same with the slave output for coordiantes trajectory file
1852  //but cell data is not needed because this is velocity trajectory
1853  int atomsRemains = (totalAtoms-1)-(tID+1)+1;
1854  off_t offset = ((off_t)atomsRemains)*sizeof(float)+1*sizeof(int);
1855  seek_dcdfile(veldcdFileID, offset, SEEK_CUR);
1856 
1857 #else
1858  //write the timestep
1859  NAMD_write(veldcdFileID, (char *)&timestep, sizeof(int));
1860  //write the values for this timestep
1861  NAMD_write(veldcdFileID, (char *)vecs, sizeof(Vector)*parN);
1862 #endif
1863 }
1864 
1865 void ParOutput::output_restart_velocities_master(int timestep, int n){
1866 #if OUTPUT_SINGLE_FILE
1867  char timestepstr[20];
1868 
1869  int baselen = strlen(namdMyNode->simParams->restartFilename);
1870  char *restart_name = new char[baselen+26];
1871 
1872  strcpy(restart_name, namdMyNode->simParams->restartFilename);
1873  if ( namdMyNode->simParams->restartSave ) {
1874  sprintf(timestepstr,".%d",timestep);
1875  strcat(restart_name, timestepstr);
1876  }
1877  strcat(restart_name, ".vel");
1878 #else
1879  char *restart_name = NULL;
1880  if ( namdMyNode->simParams->restartSave )
1881  restart_name = buildFileName(velType);
1882  else
1883  restart_name = buildFileName(velType,timestep);
1884 #endif
1885 
1886  NAMD_backup_file(restart_name,".old");
1887 
1888  //Always output a binary file
1889  write_binary_file_master(restart_name, n);
1890 
1891  delete [] restart_name;
1892 }
1893 
1894 void ParOutput::output_restart_velocities_slave(int timestep, int fID, int tID, Vector *vecs, int64 offset){
1895 #if OUTPUT_SINGLE_FILE
1896  char timestepstr[20];
1897 
1898  int baselen = strlen(namdMyNode->simParams->restartFilename);
1899  char *restart_name = new char[baselen+26];
1900 
1901  strcpy(restart_name, namdMyNode->simParams->restartFilename);
1902  if ( namdMyNode->simParams->restartSave ) {
1903  sprintf(timestepstr,".%d",timestep);
1904  strcat(restart_name, timestepstr);
1905  }
1906  strcat(restart_name, ".vel");
1907 #else
1908  char *restart_name = NULL;
1909  if ( namdMyNode->simParams->restartSave )
1910  restart_name = buildFileName(velType);
1911  else
1912  restart_name = buildFileName(velType,timestep);
1913 
1914  NAMD_backup_file(restart_name,".old");
1915 #endif
1916 
1917  //Always output a binary file
1918  write_binary_file_slave(restart_name, fID, tID, vecs, offset);
1919 
1920  delete [] restart_name;
1921 }
1922 
1923 void ParOutput::output_final_velocities_master(int n){
1924 #if OUTPUT_SINGLE_FILE
1925  char *output_name = new char[strlen(namdMyNode->simParams->outputFilename)+8];
1926  // Build the output filename
1927  strcpy(output_name, namdMyNode->simParams->outputFilename);
1928  strcat(output_name, ".vel");
1929 #else
1930  char *output_name = buildFileName(velType);
1931 #endif
1932 
1933  NAMD_backup_file(output_name);
1934 
1935  //Write the velocities to a binary file
1936  write_binary_file_master(output_name, n);
1937 }
1938 
1939 void ParOutput::output_final_velocities_slave(int fID, int tID, Vector *vecs, int64 offset){
1940 #if OUTPUT_SINGLE_FILE
1941  char *output_name = new char[strlen(namdMyNode->simParams->outputFilename)+8];
1942  // Build the output filename
1943  strcpy(output_name, namdMyNode->simParams->outputFilename);
1944  strcat(output_name, ".vel");
1945 #else
1946  char *output_name = buildFileName(velType);
1947  NAMD_backup_file(output_name);
1948 #endif
1949 
1950  //Write the velocities to a binary file
1951  write_binary_file_slave(output_name, fID, tID, vecs, offset);
1952 
1953  delete [] output_name;
1954 }
1955 
1956 void ParOutput::write_binary_file_master(char *fname, int n){
1957  char errmsg[256];
1958  int fd; // File descriptor
1959  int32 n32 = n;
1960 
1961  fd = NAMD_open(fname);
1962 
1963  sprintf(errmsg, "Error on write to binary file %s", fname);
1964 
1965  #if !OUTPUT_SINGLE_FILE
1966  // Write out extra fields as MAGIC number etc.
1967  int32 tmpInt = OUTPUT_MAGIC_NUMBER;
1968  float tmpFlt = OUTPUT_FILE_VERSION;
1969  NAMD_write(fd, (char *) &tmpInt, sizeof(int32), errmsg);
1970  NAMD_write(fd, (char *) &tmpFlt, sizeof(float), errmsg);
1971  tmpInt = namdMyNode->simParams->numoutputprocs;
1972  NAMD_write(fd, (char *) &tmpInt, sizeof(int32), errmsg);
1973  #endif
1974 
1975  // Write out the number of atoms and the vectors
1976  NAMD_write(fd, (char *) &n32, sizeof(int32), errmsg);
1977 
1978  NAMD_close(fd, fname);
1979 }
1980 
1981 void ParOutput::write_binary_file_slave(char *fname, int fID, int tID, Vector *vecs, int64 offset){
1982  char errmsg[256];
1983 
1984 #if OUTPUT_SINGLE_FILE
1985  //the mode has to be "r+" because the file already exists
1986  FILE *ofp = fopen(fname, "rb+");
1987  if ( ! ofp ) {
1988  sprintf(errmsg, "Error on opening binary file %s", fname);
1989  NAMD_err(errmsg);
1990  }
1991 
1992  //if the output is a single file, then the file position needs to be set correctly
1993 #ifdef WIN32
1994  if ( _fseeki64(ofp, offset, SEEK_SET) )
1995 #else
1996  if ( fseeko(ofp, offset, SEEK_SET) )
1997 #endif
1998  {
1999  sprintf(errmsg, "Error on seeking binary file %s", fname);
2000  NAMD_err(errmsg);
2001  }
2002 #else
2003  //the mode has to be "w+" because the file doesn't exist yet
2004  FILE *ofp = fopen(fname, "wb+");
2005  if ( ! ofp ) {
2006  sprintf(errmsg, "Error on opening binary file %s", fname);
2007  NAMD_err(errmsg);
2008  }
2009 
2010  // Write out extra fields as MAGIC number etc.
2011  int32 tmpInt = OUTPUT_MAGIC_NUMBER;
2012  float tmpFlt = OUTPUT_FILE_VERSION;
2013  fwrite(&tmpInt, sizeof(int32), 1, ofp);
2014  fwrite(&tmpFlt, sizeof(float), 1, ofp);
2015  fwrite(&outputID, sizeof(int), 1, ofp);
2016  fwrite(&fID, sizeof(int), 1, ofp);
2017  fwrite(&tID, sizeof(int), 1, ofp);
2018 #endif
2019 
2020  int parN = tID-fID+1;
2021  if ( fwrite(vecs, sizeof(Vector), parN, ofp) != parN ) {
2022  sprintf(errmsg, "Error on writing to binary file %s", fname);
2023  NAMD_err(errmsg);
2024  }
2025 
2026  if ( fclose(ofp) ) {
2027  sprintf(errmsg, "Error on closing binary file %s", fname);
2028  NAMD_err(errmsg);
2029  }
2030 }
2032 
2033 
2035 void ParOutput::forceMaster(int timestep, int n){
2036  SimParameters *simParams = Node::Object()->simParameters;
2037 
2038  if ( timestep >= 0 ) {
2039 
2040  // Output force DCD trajectory
2041  if ( simParams->forceDcdFrequency &&
2042  ((timestep % simParams->forceDcdFrequency) == 0) )
2043  {
2044  output_forcedcdfile_master(timestep, n);
2045  }
2046 
2047  }
2048 
2049  // Output forces
2050  if (timestep == FORCE_OUTPUT)
2051  {
2052  int realstep = simParams->firstTimestep;
2053  iout << "WRITING FORCES TO OUTPUT FILE AT STEP "
2054  << realstep << "\n" << endi;
2055  fflush(stdout);
2056  output_forces_master(n);
2057  }
2058 
2059  // Close trajectory files in velocityMaster above
2060 }
2061 
2062 //output atoms' forces from id fID to tID.
2063 void ParOutput::forceSlave(int timestep, int fID, int tID, Vector *vecs){
2064  SimParameters *simParams = Node::Object()->simParameters;
2065 
2066  if ( timestep >= 0 ) {
2067 
2068  // Output force DCD trajectory
2069  if ( simParams->forceDcdFrequency &&
2070  ((timestep % simParams->forceDcdFrequency) == 0) )
2071  {
2072  output_forcedcdfile_slave(timestep, fID, tID, vecs);
2073  }
2074 
2075  }
2076 
2077  // Output forces
2078  if (timestep == FORCE_OUTPUT)
2079  {
2080  int64 offset = sizeof(int)+sizeof(Vector)*((int64)fID);
2081  output_forces_slave(fID, tID, vecs, offset);
2082  }
2083 
2084  // Close trajectory files in velocitySlave above
2085 }
2086 
2087 void ParOutput::output_forcedcdfile_master(int timestep, int n){
2088  int ret_code; // Return code from DCD calls
2089  SimParameters *simParams = Node::Object()->simParameters;
2090 
2091  // If this is the last time we will be writing coordinates,
2092  // close the file before exiting
2093  if ( timestep == END_OF_RUN ) {
2094  if ( ! forcedcdFirst ) {
2095  iout << "CLOSING FORCE DCD FILE\n" << endi;
2096  close_dcd_write(forcedcdFileID);
2097  } else {
2098  iout << "FORCE DCD FILE WAS NOT CREATED\n" << endi;
2099  }
2100  return;
2101  }
2102 
2103  if (forcedcdFirst)
2104  {
2105  // Open the DCD file
2106  iout << "OPENING FORCE DCD FILE\n" << endi;
2107 
2108  #if OUTPUT_SINGLE_FILE
2109  char *forcedcdFilename = simParams->forceDcdFilename;
2110  #else
2111  char *forcedcdFilename = buildFileName(forcedcdType);
2112  #endif
2113 
2114  forcedcdFileID=open_dcd_write(forcedcdFilename);
2115 
2116  if (forcedcdFileID == DCD_FILEEXISTS)
2117  {
2118  char err_msg[257];
2119  sprintf(err_msg, "Force DCD file %s already exists!",forcedcdFilename);
2120  NAMD_err(err_msg);
2121  }
2122  else if (forcedcdFileID < 0)
2123  {
2124  char err_msg[257];
2125  sprintf(err_msg, "Couldn't open force DCD file %s",forcedcdFilename);
2126  NAMD_err(err_msg);
2127  }
2128 
2129  #if !OUTPUT_SINGLE_FILE
2130  // Write out extra fields as MAGIC number etc.
2131  int32 tmpInt = OUTPUT_MAGIC_NUMBER;
2132  float tmpFlt = OUTPUT_FILE_VERSION;
2133  NAMD_write(forcedcdFileID, (char *) &tmpInt, sizeof(int32));
2134  NAMD_write(forcedcdFileID, (char *) &tmpFlt, sizeof(float));
2135  tmpInt = simParams->numoutputprocs;
2136  NAMD_write(forcedcdFileID, (char *) &tmpInt, sizeof(int32));
2137  #endif
2138 
2139  int NSAVC, NFILE, NPRIV, NSTEP;
2140  NSAVC = simParams->forceDcdFrequency;
2141  NPRIV = timestep;
2142  NSTEP = NPRIV - NSAVC;
2143  NFILE = 0;
2144 
2145  // Write out the header
2146  const int with_unitcell = 0;
2147  ret_code = write_dcdheader(forcedcdFileID,
2148  forcedcdFilename,
2149  n, NFILE, NPRIV, NSAVC, NSTEP,
2150  simParams->dt/TIMEFACTOR, with_unitcell);
2151 
2152 
2153  if (ret_code<0)
2154  {
2155  NAMD_err("Writing of force DCD header failed!!");
2156  }
2157 
2158  #if !OUTPUT_SINGLE_FILE
2159  //In this case, the file name is dynamically allocated
2160  delete [] forcedcdFilename;
2161  #endif
2162 
2163  forcedcdFirst = FALSE;
2164  }
2165 
2166  // Write out the values for this timestep
2167  iout << "WRITING FORCES TO DCD FILE AT STEP "
2168  << timestep << "\n" << endi;
2169  fflush(stdout);
2170 
2171  //In the case of writing to multiple files, only the header
2172  //of the dcd file needs to be updated. Note that the format of
2173  //the new dcd file has also changed! -Chao Mei
2174 
2175 #if OUTPUT_SINGLE_FILE
2176  //write X,Y,Z headers
2177  int totalAtoms = namdMyNode->molecule->numAtoms;
2178  write_dcdstep_par_XYZUnits(forcedcdFileID, totalAtoms);
2179 #endif
2180 
2181  //update the header
2182  update_dcdstep_par_header(forcedcdFileID);
2183 
2184 }
2185 
2186 void ParOutput::output_forcedcdfile_slave(int timestep, int fID, int tID, Vector *vecs){
2187  int ret_code; // Return code from DCD calls
2188  SimParameters *simParams = Node::Object()->simParameters;
2189 
2190  // If this is the last time we will be writing coordinates,
2191  // close the file before exiting
2192  if ( timestep == END_OF_RUN ) {
2193  if ( ! forcedcdFirst ) {
2194  close_dcd_write(forcedcdFileID);
2195  }
2196 #if OUTPUT_SINGLE_FILE
2197  delete [] forcedcdX;
2198  delete [] forcedcdY;
2199  delete [] forcedcdZ;
2200 #endif
2201  return;
2202  }
2203 
2204  int parN = tID-fID+1;
2205 
2206  if (forcedcdFirst)
2207  {
2208 
2209  #if OUTPUT_SINGLE_FILE
2210  char *forcedcdFilename = namdMyNode->simParams->forceDcdFilename;
2211  #else
2212  char *forcedcdFilename = buildFileName(forcedcdType);
2213  #endif
2214  forcedcdFileID=open_dcd_write_par_slave(forcedcdFilename);
2215  if(forcedcdFileID < 0)
2216  {
2217  char err_msg[257];
2218  sprintf(err_msg, "Couldn't open force DCD file %s",forcedcdFilename);
2219  NAMD_err(err_msg);
2220  }
2221  #if OUTPUT_SINGLE_FILE
2222  //If outputting to a single file, dcd files conforms to the old format
2223  //as data are organized as 3 seperate arrays of X,Y,Z, while in the new
2224  //format used in outputing multiple files, the data are organized as an
2225  //array of XYZs.
2226  forcedcdX = new float[parN];
2227  forcedcdY = new float[parN];
2228  forcedcdZ = new float[parN];
2229  //seek to beginning of X,Y,Z sections which means skipping header.
2230  //Cell data is not needed because this is force trajectory
2231  int skipbytes = get_dcdheader_size();
2232  seek_dcdfile(forcedcdFileID, skipbytes, SEEK_SET);
2233  #endif
2234 
2235  #if !OUTPUT_SINGLE_FILE
2236  delete [] forcedcdFilename;
2237  // Write out extra fields as MAGIC number etc.
2238  int32 tmpInt = OUTPUT_MAGIC_NUMBER;
2239  float tmpFlt = OUTPUT_FILE_VERSION;
2240  NAMD_write(forcedcdFileID, (char *) &tmpInt, sizeof(int32));
2241  NAMD_write(forcedcdFileID, (char *) &tmpFlt, sizeof(float));
2242  NAMD_write(forcedcdFileID, (char *) &outputID, sizeof(int));
2243  NAMD_write(forcedcdFileID, (char *) &fID, sizeof(int));
2244  NAMD_write(forcedcdFileID, (char *) &tID, sizeof(int));
2245  #endif
2246 
2247  forcedcdFirst = FALSE;
2248  }
2249 
2250 #if OUTPUT_SINGLE_FILE
2251  //The following seek will set the stream position to the
2252  //beginning of the place where a new timestep output should
2253  //be performed.
2254  CmiAssert(sizeof(off_t)==8);
2255  int totalAtoms = namdMyNode->molecule->numAtoms;
2256 
2257  for(int i=0; i<parN; i++){
2258  forcedcdX[i] = vecs[i].x;
2259  forcedcdY[i] = vecs[i].y;
2260  forcedcdZ[i] = vecs[i].z;
2261  }
2262 
2263  write_dcdstep_par_slave(forcedcdFileID, fID, tID, totalAtoms, forcedcdX, forcedcdY, forcedcdZ);
2264  //same with the slave output for coordiantes trajectory file
2265  //but cell data is not needed because this is force trajectory
2266  int atomsRemains = (totalAtoms-1)-(tID+1)+1;
2267  off_t offset = ((off_t)atomsRemains)*sizeof(float)+1*sizeof(int);
2268  seek_dcdfile(forcedcdFileID, offset, SEEK_CUR);
2269 #else
2270  //write the timestep
2271  NAMD_write(forcedcdFileID, (char *)&timestep, sizeof(int));
2272  //write the values for this timestep
2273  NAMD_write(forcedcdFileID, (char *)vecs, sizeof(Vector)*parN);
2274 #endif
2275 }
2276 
2277 void ParOutput::output_forces_master(int n){
2278 #if OUTPUT_SINGLE_FILE
2279  char *output_name = new char[strlen(namdMyNode->simParams->outputFilename)+8];
2280  // Build the output filename
2281  strcpy(output_name, namdMyNode->simParams->outputFilename);
2282  strcat(output_name, ".force");
2283 #else
2284  char *output_name = buildFileName(forceType);
2285 #endif
2286 
2287  NAMD_backup_file(output_name);
2288 
2289  //Write the force to a binary file
2290  write_binary_file_master(output_name, n);
2291 }
2292 
2293 void ParOutput::output_forces_slave(int fID, int tID, Vector *vecs, int64 offset){
2294 #if OUTPUT_SINGLE_FILE
2295  char *output_name = new char[strlen(namdMyNode->simParams->outputFilename)+8];
2296  // Build the output filename
2297  strcpy(output_name, namdMyNode->simParams->outputFilename);
2298  strcat(output_name, ".force");
2299 #else
2300  char *output_name = buildFileName(forceType);
2301  NAMD_backup_file(output_name);
2302 #endif
2303 
2304  //Write the forces to a binary file
2305  write_binary_file_slave(output_name, fID, tID, vecs, offset);
2306 
2307  delete [] output_name;
2308 }
2310 
2311 
2313 void ParOutput::coordinateMaster(int timestep, int n, Lattice &lat){
2314  SimParameters *simParams = Node::Object()->simParameters;
2315 
2316  if ( timestep >= 0 ) {
2317  // Output a DCD trajectory
2318  if ( simParams->dcdFrequency &&
2319  ((timestep % simParams->dcdFrequency) == 0) )
2320  {
2321  output_dcdfile_master(timestep, n,
2322  simParams->dcdUnitCell ? &lat : NULL);
2323  }
2324 
2325  // Output a restart file
2326  if ( simParams->restartFrequency &&
2327  ((timestep % simParams->restartFrequency) == 0) )
2328  {
2329  iout << "WRITING COORDINATES TO RESTART FILE AT STEP "
2330  << timestep << "\n" << endi;
2331  output_restart_coordinates_master(timestep, n);
2332  iout << "FINISHED WRITING RESTART COORDINATES\n" <<endi;
2333  fflush(stdout);
2334  }
2335 
2336 /* Interactive MD is not supported in Parallel IO
2337  // Interactive MD
2338  if ( simParams->IMDon &&
2339  ( ((timestep % simParams->IMDfreq) == 0) ||
2340  (timestep == simParams->firstTimestep) ) )
2341  {
2342  IMDOutput *imd = Node::Object()->imd;
2343  wrap_coor(fcoor,lattice,&fcoor_wrapped);
2344  if (imd != NULL) imd->gather_coordinates(timestep, n, fcoor);
2345  }
2346 */
2347  }
2348 
2349 /* EVAL_MEASURE of a timestep is not supported in Parallel IO
2350  if (timestep == EVAL_MEASURE)
2351  {
2352  #ifdef NAMD_TCL
2353  wrap_coor(coor,lattice,&coor_wrapped);
2354  Node::Object()->getScript()->measure(coor);
2355  #endif
2356  }
2357 */
2358  // Output final coordinates
2359  if (timestep == FILE_OUTPUT || timestep == END_OF_RUN)
2360  {
2361  int realstep = ( timestep == FILE_OUTPUT ?
2362  simParams->firstTimestep : simParams->N );
2363  iout << "WRITING COORDINATES TO OUTPUT FILE AT STEP "
2364  << realstep << "\n" << endi;
2365  fflush(stdout);
2366  output_final_coordinates_master(n);
2367  }
2368 
2369  // Close trajectory files
2370  if (timestep == END_OF_RUN)
2371  {
2372  if (simParams->dcdFrequency) output_dcdfile_master(END_OF_RUN,0,NULL);
2373  }
2374 }
2375 
2376 void ParOutput::coordinateSlave(int timestep, int fID, int tID, Vector *vecs, FloatVector *fvecs){
2377  SimParameters *simParams = Node::Object()->simParameters;
2378 
2379  if ( timestep >= 0 ) {
2380  // Output a DCD trajectory
2381  if ( simParams->dcdFrequency &&
2382  ((timestep % simParams->dcdFrequency) == 0) )
2383  {
2384  output_dcdfile_slave(timestep, fID, tID, fvecs);
2385  }
2386 
2387  // Output a restart file
2388  if ( simParams->restartFrequency &&
2389  ((timestep % simParams->restartFrequency) == 0) )
2390  {
2391  int64 offset = sizeof(int)+sizeof(Vector)*((int64)fID);
2392  output_restart_coordinates_slave(timestep, fID, tID, vecs, offset);
2393  }
2394 
2395 /* Interactive MD is not supported in Parallel IO
2396  // Interactive MD
2397  if ( simParams->IMDon &&
2398  ( ((timestep % simParams->IMDfreq) == 0) ||
2399  (timestep == simParams->firstTimestep) ) )
2400  {
2401  IMDOutput *imd = Node::Object()->imd;
2402  wrap_coor(fcoor,lattice,&fcoor_wrapped);
2403  if (imd != NULL) imd->gather_coordinates(timestep, n, fcoor);
2404  }
2405 */
2406  }
2407 
2408 /* EVAL_MEASURE of a timestep is not supported in Parallel IO
2409  if (timestep == EVAL_MEASURE)
2410  {
2411  #ifdef NAMD_TCL
2412  wrap_coor(coor,lattice,&coor_wrapped);
2413  Node::Object()->getScript()->measure(coor);
2414  #endif
2415  }
2416 */
2417  // Output final coordinates
2418  if (timestep == FILE_OUTPUT || timestep == END_OF_RUN)
2419  {
2420  int64 offset = sizeof(int)+sizeof(Vector)*((int64)fID);
2421  output_final_coordinates_slave(fID, tID, vecs, offset);
2422  }
2423 
2424  // Close trajectory files
2425  if (timestep == END_OF_RUN)
2426  {
2427  if (simParams->dcdFrequency) output_dcdfile_slave(END_OF_RUN,0,0,NULL);
2428  }
2429 }
2430 
2431 void ParOutput::output_dcdfile_master(int timestep, int n, const Lattice *lattice){
2432  int ret_code; // Return code from DCD calls
2433  SimParameters *simParams = namdMyNode->simParams;
2434 
2435  // If this is the last time we will be writing coordinates,
2436  // close the file before exiting
2437  if ( timestep == END_OF_RUN ) {
2438  if ( ! dcdFirst ) {
2439  iout << "CLOSING COORDINATE DCD FILE\n" << endi;
2440  close_dcd_write(dcdFileID);
2441  } else {
2442  iout << "COORDINATE DCD FILE WAS NOT CREATED\n" << endi;
2443  }
2444  return;
2445  }
2446 
2447  if (dcdFirst)
2448  {
2449  // Open the DCD file
2450  iout << "OPENING COORDINATE DCD FILE\n" << endi;
2451 
2452  #if OUTPUT_SINGLE_FILE
2453  char *dcdFilename = simParams->dcdFilename;
2454  #else
2455  char *dcdFilename = buildFileName(dcdType);
2456  #endif
2457 
2458 
2459  dcdFileID=open_dcd_write(dcdFilename);
2460 
2461  if (dcdFileID == DCD_FILEEXISTS)
2462  {
2463  char err_msg[257];
2464  sprintf(err_msg, "DCD file %s already exists!!",dcdFilename);
2465  NAMD_err(err_msg);
2466  }
2467  else if (dcdFileID < 0)
2468  {
2469  char err_msg[257];
2470  sprintf(err_msg, "Couldn't open DCD file %s",dcdFilename);
2471  NAMD_err(err_msg);
2472  }
2473 
2474  #if !OUTPUT_SINGLE_FILE
2475  // Write out extra fields as MAGIC number etc.
2476  int32 tmpInt = OUTPUT_MAGIC_NUMBER;
2477  float tmpFlt = OUTPUT_FILE_VERSION;
2478  NAMD_write(dcdFileID, (char *) &tmpInt, sizeof(int32));
2479  NAMD_write(dcdFileID, (char *) &tmpFlt, sizeof(float));
2480  tmpInt = simParams->numoutputprocs;
2481  NAMD_write(dcdFileID, (char *) &tmpInt, sizeof(int32));
2482  #endif
2483 
2484  int NSAVC, NFILE, NPRIV, NSTEP;
2485  NSAVC = simParams->dcdFrequency;
2486  NPRIV = timestep;
2487  NSTEP = NPRIV - NSAVC;
2488  NFILE = 0;
2489 
2490  // Write out the header
2491  ret_code = write_dcdheader(dcdFileID,
2492  dcdFilename,
2493  n, NFILE, NPRIV, NSAVC, NSTEP,
2494  simParams->dt/TIMEFACTOR, lattice != NULL);
2495 
2496 
2497  if (ret_code<0)
2498  {
2499  NAMD_err("Writing of DCD header failed!!");
2500  }
2501 
2502  #if !OUTPUT_SINGLE_FILE
2503  //dcdFilename needs to be freed as it is dynamically allocated
2504  delete [] dcdFilename;
2505  #endif
2506 
2507  dcdFirst = FALSE;
2508  }
2509 
2510  // Write out the values for this timestep
2511  iout << "WRITING COORDINATES TO DCD FILE AT STEP "
2512  << timestep << "\n" << endi;
2513  fflush(stdout);
2514 
2515  //In the case of writing to multiple files, the header of the
2516  //dcd file needs to be updated. In addition, the lattice data
2517  //needs to be written if necessary. Note that the format of
2518  //the new dcd file has also changed! -Chao Mei
2519 
2520  // Write out the Cell data
2521  if (lattice) {
2522  double unitcell[6];
2523  lattice_to_unitcell(lattice,unitcell);
2524  write_dcdstep_par_cell(dcdFileID, unitcell);
2525  }
2526 
2527 #if OUTPUT_SINGLE_FILE
2528  //write X,Y,Z headers
2529  int totalAtoms = namdMyNode->molecule->numAtoms;
2530  write_dcdstep_par_XYZUnits(dcdFileID, totalAtoms);
2531 #endif
2532 
2533  //update the header
2534  update_dcdstep_par_header(dcdFileID);
2535 }
2536 void ParOutput::output_dcdfile_slave(int timestep, int fID, int tID, FloatVector *fvecs){
2537  int ret_code; // Return code from DCD calls
2538  SimParameters *simParams = Node::Object()->simParameters;
2539 
2540  // If this is the last time we will be writing coordinates,
2541  // close the file before exiting
2542  if ( timestep == END_OF_RUN ) {
2543  if ( ! dcdFirst ) {
2544  close_dcd_write(dcdFileID);
2545  }
2546 #if OUTPUT_SINGLE_FILE
2547  delete [] dcdX;
2548  delete [] dcdY;
2549  delete [] dcdZ;
2550 #endif
2551  return;
2552  }
2553 
2554  int parN = tID-fID+1;
2555 
2556  if (dcdFirst)
2557  {
2558 
2559  #if OUTPUT_SINGLE_FILE
2560  char *dcdFilename = simParams->dcdFilename;
2561  #else
2562  char *dcdFilename = buildFileName(dcdType);
2563  #endif
2564  dcdFileID=open_dcd_write_par_slave(dcdFilename);
2565  if(dcdFileID < 0)
2566  {
2567  char err_msg[257];
2568  sprintf(err_msg, "Couldn't open DCD file %s", dcdFilename);
2569  NAMD_err(err_msg);
2570  }
2571 
2572  #if OUTPUT_SINGLE_FILE
2573  dcdX = new float[parN];
2574  dcdY = new float[parN];
2575  dcdZ = new float[parN];
2576  //seek to beginning of X,Y,Z sections which means skipping header
2577  //skip the cell data if necessary
2578  int skipbytes = get_dcdheader_size();
2579  if(simParams->dcdUnitCell) {
2580  skipbytes += sizeof(int)*2 + 6*sizeof(double);
2581  }
2582  seek_dcdfile(dcdFileID, skipbytes, SEEK_SET);
2583  #endif
2584 
2585  #if !OUTPUT_SINGLE_FILE
2586  delete [] dcdFilename;
2587 
2588  // Write out extra fields as MAGIC number etc.
2589  int32 tmpInt = OUTPUT_MAGIC_NUMBER;
2590  float tmpFlt = OUTPUT_FILE_VERSION;
2591  NAMD_write(dcdFileID, (char *) &tmpInt, sizeof(int32));
2592  NAMD_write(dcdFileID, (char *) &tmpFlt, sizeof(float));
2593  NAMD_write(dcdFileID, (char *) &outputID, sizeof(int));
2594  NAMD_write(dcdFileID, (char *) &fID, sizeof(int));
2595  NAMD_write(dcdFileID, (char *) &tID, sizeof(int));
2596  #endif
2597  dcdFirst = FALSE;
2598  }
2599 
2600 #if OUTPUT_SINGLE_FILE
2601  //The following seek will set the stream position to the
2602  //beginning of the place where a new timestep output should
2603  //be performed.
2604  CmiAssert(sizeof(off_t)==8);
2605  int totalAtoms = namdMyNode->molecule->numAtoms;
2606 
2607  for(int i=0; i<parN; i++){
2608  dcdX[i] = fvecs[i].x;
2609  dcdY[i] = fvecs[i].y;
2610  dcdZ[i] = fvecs[i].z;
2611  }
2612 
2613  write_dcdstep_par_slave(dcdFileID, fID, tID, totalAtoms, dcdX, dcdY, dcdZ);
2614 
2615  //1. Always foward to the beginning position of X,Y,Z sections in the
2616  //next timeframe.
2617  //2. SHOULD AVOID USING SEEK_END: although the file size is updated on the
2618  //master proc, slave procs should rely on this update by seeking
2619  //from SEEK_END because such info may not be updated in a timely manner
2620  //when this slave proc needs it.
2621 
2622  //We know the current position of file after slave writing is at the
2623  //end of Z output for atom "tID" (tID is the atom id indexed from 0)
2624  //(totalAtoms-1) is the last atom id, (tID+1) is the next atom id
2625  int atomsRemains = (totalAtoms-1)-(tID+1)+1;
2626  off_t offset = ((off_t)atomsRemains)*sizeof(float)+1*sizeof(int);
2627  //then skip the cell data if necessary
2628  if(simParams->dcdUnitCell) {
2629  offset += sizeof(int)*2 + 6*sizeof(double);
2630  }
2631  seek_dcdfile(dcdFileID, offset, SEEK_CUR);
2632 #else
2633  //write the timestep
2634  NAMD_write(dcdFileID, (char *)&timestep, sizeof(int));
2635  //write the values for this timestep
2636  NAMD_write(dcdFileID, (char *)fvecs, sizeof(FloatVector)*parN);
2637 #endif
2638 }
2639 
2640 void ParOutput::output_restart_coordinates_master(int timestep, int n){
2641 #if OUTPUT_SINGLE_FILE
2642  char timestepstr[20];
2643 
2644  int baselen = strlen(namdMyNode->simParams->restartFilename);
2645  char *restart_name = new char[baselen+26];
2646 
2647  strcpy(restart_name, namdMyNode->simParams->restartFilename);
2648  if ( namdMyNode->simParams->restartSave ) {
2649  sprintf(timestepstr,".%d",timestep);
2650  strcat(restart_name, timestepstr);
2651  }
2652  strcat(restart_name, ".coor");
2653 #else
2654  char *restart_name = NULL;
2655  if ( namdMyNode->simParams->restartSave )
2656  restart_name = buildFileName(coorType);
2657  else
2658  restart_name = buildFileName(coorType,timestep);
2659 #endif
2660 
2661  NAMD_backup_file(restart_name,".old");
2662 
2663  // Generate a binary restart file
2664  write_binary_file_master(restart_name, n);
2665 
2666  delete [] restart_name;
2667 }
2668 void ParOutput::output_restart_coordinates_slave(int timestep, int fID, int tID, Vector *vecs, int64 offset){
2669 #if OUTPUT_SINGLE_FILE
2670  char timestepstr[20];
2671 
2672  int baselen = strlen(namdMyNode->simParams->restartFilename);
2673  char *restart_name = new char[baselen+26];
2674 
2675  strcpy(restart_name, namdMyNode->simParams->restartFilename);
2676  if ( namdMyNode->simParams->restartSave ) {
2677  sprintf(timestepstr,".%d",timestep);
2678  strcat(restart_name, timestepstr);
2679  }
2680  strcat(restart_name, ".coor");
2681 #else
2682  char *restart_name = NULL;
2683  if ( namdMyNode->simParams->restartSave )
2684  restart_name = buildFileName(coorType);
2685  else
2686  restart_name = buildFileName(coorType,timestep);
2687 
2688  NAMD_backup_file(restart_name,".old");
2689 #endif
2690 
2691  // Generate a binary restart file
2692  write_binary_file_slave(restart_name, fID, tID, vecs, offset);
2693 
2694  delete [] restart_name;
2695 }
2696 
2697 void ParOutput::output_final_coordinates_master(int n){
2698 #if OUTPUT_SINGLE_FILE
2699  char *output_name = new char[strlen(namdMyNode->simParams->outputFilename)+8];
2700 
2701  // Built the output filename
2702  strcpy(output_name, namdMyNode->simParams->outputFilename);
2703  strcat(output_name, ".coor");
2704 #else
2705  char *output_name = buildFileName(coorType);
2706 #endif
2707 
2708  NAMD_backup_file(output_name);
2709 
2710  // Write the coordinates to a binary file
2711  write_binary_file_master(output_name, n);
2712 
2713  delete [] output_name;
2714 }
2715 void ParOutput::output_final_coordinates_slave(int fID, int tID, Vector *vecs, int64 offset){
2716 #if OUTPUT_SINGLE_FILE
2717  char *output_name = new char[strlen(namdMyNode->simParams->outputFilename)+8];
2718 
2719  // Built the output filename
2720  strcpy(output_name, namdMyNode->simParams->outputFilename);
2721  strcat(output_name, ".coor");
2722 #else
2723  char *output_name = buildFileName(coorType);
2724  NAMD_backup_file(output_name);
2725 #endif
2726 
2727  // Write the coordinates to a binary file
2728  write_binary_file_slave(output_name, fID, tID, vecs, offset);
2729 
2730  delete [] output_name;
2731 }
2733 
2735 #if !OUTPUT_SINGLE_FILE
2736 char *ParOutput::buildFileName(OUTPUTFILETYPE type, int timestep){
2737  char *filename = NULL;
2738  const char *typeName = NULL;
2739  switch(type) {
2740  case dcdType:
2741  typeName = "dcd";
2742  break;
2743  case forcedcdType:
2744  typeName = "forcedcd";
2745  break;
2746  case veldcdType:
2747  typeName = "veldcd";
2748  break;
2749  case coorType:
2750  typeName = "coor";
2751  break;
2752  case forceType:
2753  typeName = "force";
2754  break;
2755  case velType:
2756  typeName = "vel";
2757  break;
2758  default:
2759  typeName = "invalid";
2760  break;
2761  }
2762  int baselen = strlen(namdMyNode->simParams->outputFilename);
2763  filename = new char[baselen+32];
2764  memset(filename, 0, baselen+32);
2765 
2766 #if 0
2767  strcpy(filename, namdMyNode->simParams->outputFilename);
2768 
2769  //check if the directory exists or not
2770  if(access(filename, F_OK)!=0) {
2771  int ret = MKDIR(filename);
2772  if(ret!=0) {
2773  char errmsg[512];
2774  sprintf(errmsg, "Error in creating top-level directory %s!", filename);
2775  NAMD_err(errmsg);
2776  }
2777  }
2778 
2779  strcat(filename, PATHSEPSTR);
2780  strcat(filename, typeName);
2781 
2782  //check if the directory exists or not
2783  if(access(filename, F_OK)!=0) {
2784  int ret = MKDIR(filename);
2785  if(ret!=0) {
2786  char errmsg[512];
2787  sprintf(errmsg, "Error in creating middle-level directory %s!", filename);
2788  NAMD_err(errmsg);
2789  }
2790  }
2791 #else
2792  sprintf(filename, "%s%s%s", namdMyNode->simParams->outputFilename, PATHSEPSTR, typeName);
2793 #endif
2794 
2795  char tmpstr[20];
2796  if(outputID == -1) {
2797  //indicating the output from master
2798  if(timestep!=-9999) {
2799  //not the default value
2800  sprintf(tmpstr, "%smeta.%d", PATHSEPSTR, timestep);
2801  }else{
2802  sprintf(tmpstr, "%smeta", PATHSEPSTR);
2803  }
2804  }else{
2805  //indicating the output from slave
2806  sprintf(tmpstr, "%s%d", PATHSEPSTR, outputID);
2807  strcat(filename, tmpstr);
2808  #if 0
2809  if(access(filename, F_OK)!=0) {
2810  int ret = MKDIR(filename);
2811  if(ret!=0) {
2812  char errmsg[512];
2813  sprintf(errmsg, "Error in creating last-level directory %s!", filename);
2814  NAMD_err(errmsg);
2815  }
2816  }
2817  #endif
2818  if(timestep!=-9999) {
2819  //not the default value
2820  sprintf(tmpstr, "%s%s_%d.%d", PATHSEPSTR,typeName,outputID,timestep);
2821  }else{
2822  sprintf(tmpstr,"%s%s_%d", PATHSEPSTR,typeName,outputID);
2823  }
2824  }
2825 
2826  strcat(filename, tmpstr);
2827  return filename;
2828 }
2829 #endif
2830 #endif
Node::Object
static Node * Object()
Definition: Node.h:86
Molecule::dcdSelectionParams
DCDParams dcdSelectionParams[16]
Definition: Molecule.h:481
ReplicaDcdInitMsg::dcdSelectIndex
uint16_t dcdSelectIndex
Definition: DataExchanger.h:39
open_dcd_write_par_slave
int open_dcd_write_par_slave(char *dcdname)
Definition: dcdlib.C:689
dcd_params::frequency
int frequency
Definition: common.h:262
Output::velocityNeeded
static int velocityNeeded(int)
Definition: Output.C:500
Lattice::c
NAMD_HOST_DEVICE Vector c() const
Definition: Lattice.h:270
NAMD_err
void NAMD_err(const char *err_msg)
Definition: common.C:170
Molecule::get_atom_index_from_dcd_selection
int get_atom_index_from_dcd_selection(const int index, const int atomIndex)
Definition: Molecule.h:877
SimParameters::close_veldcdfile
void close_veldcdfile()
Definition: Output.C:829
Lattice::c_p
NAMD_HOST_DEVICE int c_p() const
Definition: Lattice.h:291
sendReplicaDcdInit
void sendReplicaDcdInit(int dstPart, ReplicaDcdInitMsg *msg, int msgsize)
Definition: DataExchanger.C:45
Output::setReplicaDcdIndex
void setReplicaDcdIndex(int index)
Definition: Output.C:838
update_dcdstep_par_header
int update_dcdstep_par_header(int fd)
Definition: dcdlib.C:839
FILE_OUTPUT
#define FILE_OUTPUT
Definition: Output.h:25
Node::imd
IMDOutput * imd
Definition: Node.h:186
dcd_params::tag
uint16 tag
Definition: common.h:259
EVAL_MEASURE
#define EVAL_MEASURE
Definition: Output.h:27
Vector
Definition: Vector.h:72
Node::output
Output * output
Definition: Node.h:185
Node::simParameters
SimParameters * simParameters
Definition: Node.h:181
Output::~Output
~Output()
Definition: Output.C:179
Real
float Real
Definition: common.h:118
int32
int32_t int32
Definition: common.h:38
endi
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
Vector::z
BigReal z
Definition: Vector.h:74
FloatVector::y
float y
Definition: Vector.h:26
PATHSEPSTR
#define PATHSEPSTR
Definition: Output.C:47
FALSE
#define FALSE
Definition: common.h:127
MKDIR
#define MKDIR(X)
Definition: Output.C:48
iWARN
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:82
Molecule::get_clusterSize
int get_clusterSize(int anum) const
Definition: Molecule.h:1093
Output::coordinate
void coordinate(int, int, Vector *, FloatVector *, Lattice &)
Definition: Output.C:345
Output::replicaDcdSelectInit
void replicaDcdSelectInit(int index, const char *tag, const char *filename)
Definition: Output.C:855
IMDOutput::gather_coordinates
void gather_coordinates(int timestep, int N, FloatVector *coords)
Definition: IMDOutput.C:29
OUTPUT_FILE_VERSION
#define OUTPUT_FILE_VERSION
Definition: Output.h:32
NAMD_close
#define NAMD_close
Definition: Output.C:53
iout
#define iout
Definition: InfoStream.h:51
seek_dcdfile
#define seek_dcdfile
Definition: Output.C:125
Output::forceNeeded
static int forceNeeded(int)
Definition: Output.C:595
Molecule
Molecule stores the structural information for the system.
Definition: Molecule.h:175
Lattice::b_p
NAMD_HOST_DEVICE int b_p() const
Definition: Lattice.h:290
SimParameters::close_dcdfile
void close_dcdfile()
Definition: Output.C:813
open_dcd_write
int open_dcd_write(const char *dcdname)
Definition: dcdlib.C:662
wrap_coor_int_dcd_selection
void wrap_coor_int_dcd_selection(xVector *coor, Lattice &lattice, xDone *done, int index)
Definition: Output.C:288
NAMD_open
#define NAMD_open
Definition: Output.C:51
uint16
uint16_t uint16
Definition: common.h:41
DCD_FILEEXISTS
#define DCD_FILEEXISTS
Definition: dcdlib.h:52
Molecule::find_or_create_dcd_selection_index
uint16_t find_or_create_dcd_selection_index(const char *keystr)
Molecule::get_dcd_selection_size
const int get_dcd_selection_size(const int index)
Definition: Molecule.h:886
wrap_coor
void wrap_coor(Vector *coor, Lattice &lattice, double *done)
Definition: Output.C:337
Output::recvReplicaDcdInit
void recvReplicaDcdInit(ReplicaDcdInitMsg *msg)
Definition: Output.C:869
Output
Definition: Output.h:35
dcd_params::outFilename
char outFilename[NAMD_FILENAME_BUFFER_SIZE]
Definition: common.h:261
write_dcdstep_par_XYZUnits
int write_dcdstep_par_XYZUnits(int fd, int N)
Definition: dcdlib.C:810
get_dcdheader_size
int get_dcdheader_size()
Definition: dcdlib.C:1005
FloatVector::x
float x
Definition: Vector.h:26
write_dcdstep_par_cell
int write_dcdstep_par_cell(int fd, double *cell)
Definition: dcdlib.C:785
OUTPUT_MAGIC_NUMBER
#define OUTPUT_MAGIC_NUMBER
Definition: Output.h:31
FloatVector::z
float z
Definition: Vector.h:26
sendReplicaDcdData
void sendReplicaDcdData(int dstPart, ReplicaDcdDataMsg *msg, int msgsize)
Definition: DataExchanger.C:56
Vector::length
NAMD_HOST_DEVICE BigReal length(void) const
Definition: Vector.h:202
Output::replicaDcdInit
void replicaDcdInit(int index, const char *filename)
Definition: Output.C:843
NAMD_bug
void NAMD_bug(const char *err_msg)
Definition: common.C:195
namdMyNode
#define namdMyNode
Definition: Output.C:128
Bool
int Bool
Definition: common.h:142
Molecule::get_dcd_selection_index_from_atom_id
int get_dcd_selection_index_from_atom_id(const int index, const int atomIndex)
Definition: Molecule.h:881
Output::velocity
void velocity(int, int, Vector *)
Definition: Output.C:531
Output::recvReplicaDcdData
void recvReplicaDcdData(ReplicaDcdDataMsg *msg)
Definition: Output.C:880
write_dcdstep
int write_dcdstep(int fd, int N, float *X, float *Y, float *Z, double *cell)
Definition: dcdlib.C:736
wrap_coor_dcd_selection
void wrap_coor_dcd_selection(Vector *coor, Lattice &lattice, double *done, int index)
Definition: Output.C:329
Vector::x
BigReal x
Definition: Vector.h:74
Lattice::a_p
NAMD_HOST_DEVICE int a_p() const
Definition: Lattice.h:289
Molecule::numAtoms
int numAtoms
Definition: Molecule.h:585
END_OF_RUN
#define END_OF_RUN
Definition: Output.h:26
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:147
NAMD_write
#define NAMD_write
Definition: Output.C:52
Lattice::b
NAMD_HOST_DEVICE Vector b() const
Definition: Lattice.h:269
write_dcdheader
int write_dcdheader(int fd, const char *filename, int N, int NFILE, int NPRIV, int NSAVC, int NSTEP, double DELTA, int with_unitcell)
Definition: dcdlib.C:915
Output::coordinateNeeded
static std::pair< int, int > coordinateNeeded(int)
Definition: Output.C:199
NAMD_backup_file
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:235
OUTPUTFILETYPE
OUTPUTFILETYPE
Definition: common.h:247
simParams
#define simParams
Definition: Output.C:129
O_LARGEFILE
#define O_LARGEFILE
Definition: Output.C:56
dcd_params::size
int size
Definition: common.h:266
Output::Output
Output()
Definition: Output.C:169
sendReplicaDcdAck
void sendReplicaDcdAck(int dstPart, ReplicaDcdAckMsg *msg)
Definition: DataExchanger.C:67
ReplicaDcdDataMsg::dcdSelectIndex
uint16_t dcdSelectIndex
Definition: DataExchanger.h:48
Vector::y
BigReal y
Definition: Vector.h:74
maxDCD
static const short maxDCD
Definition: Output.C:42
Lattice::wrap_nearest_delta
NAMD_HOST_DEVICE Vector wrap_nearest_delta(Position pos1) const
Definition: Lattice.h:233
Node::getScript
ScriptTcl * getScript()
Definition: Node.C:1599
write_dcdstep_par_slave
int write_dcdstep_par_slave(int fd, int parL, int parU, int N, float *X, float *Y, float *Z)
Definition: dcdlib.C:864
PDBVELFACTOR
#define PDBVELFACTOR
Definition: common.h:57
TIMEFACTOR
#define TIMEFACTOR
Definition: common.h:55
ReplicaDcdDataMsg::unitcell
double unitcell[6]
Definition: DataExchanger.h:51
Molecule::get_cluster
int get_cluster(int anum) const
Definition: Molecule.h:1092
Lattice::a
NAMD_HOST_DEVICE Vector a() const
Definition: Lattice.h:268
close_dcd_write
void close_dcd_write(int fd)
Definition: dcdlib.C:1063
Molecule::is_water
Bool is_water(int)
int64
int64_t int64
Definition: common.h:39
ScriptTcl::measure
void measure(Vector *)
Definition: ScriptTcl.C:1400
Lattice::wrap_delta
NAMD_HOST_DEVICE Vector wrap_delta(const Position &pos1) const
Definition: Lattice.h:222
Node::molecule
Molecule * molecule
Definition: Node.h:179
TRUE
#define TRUE
Definition: common.h:128
wrap_coor_int
void wrap_coor_int(xVector *coor, Lattice &lattice, xDone *done)
Definition: Output.C:254
Output::force
void force(int, int, Vector *)
Definition: Output.C:621
FORCE_OUTPUT
#define FORCE_OUTPUT
Definition: Output.h:28
lattice_to_unitcell
static void lattice_to_unitcell(const Lattice *lattice, double *unitcell)
Definition: Output.C:132
PDBVELINVFACTOR
#define PDBVELINVFACTOR
Definition: common.h:58
for
for(int i=0;i< n1;++i)
Definition: ComputeFullDirectBase.h:34
Generated by   doxygen 1.8.14