← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/colvars/colvarmodule.h

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/colvars/colvarmodule.h
 
1
/// -*- c++ -*-
 
2
 
1
3
#ifndef COLVARMODULE_H
2
4
#define COLVARMODULE_H
3
5
 
4
6
#ifndef COLVARS_VERSION
5
 
#define COLVARS_VERSION "2014-05-08"
 
7
#define COLVARS_VERSION "2015-02-27"
6
8
#endif
7
9
 
8
10
#ifndef COLVARS_DEBUG
18
20
/// shared between all object instances) to be accessed from other
19
21
/// objects.
20
22
 
 
23
// Internal method return codes
 
24
#define COLVARS_NOT_IMPLEMENTED -2
 
25
#define COLVARS_ERROR -1
 
26
#define COLVARS_OK 0
 
27
 
 
28
// On error, values of the colvars module error register
 
29
#define GENERAL_ERROR   1
 
30
#define FILE_ERROR      (1<<1)
 
31
#define MEMORY_ERROR    (1<<2)
 
32
#define BUG_ERROR       (1<<3) // Inconsistent state indicating bug
 
33
#define INPUT_ERROR     (1<<4) // out of bounds or inconsistent input
 
34
#define DELETE_COLVARS  (1<<5) // Instruct the caller to delete cvm
 
35
#define FATAL_ERROR     (1<<6) // Should be set, or not, together with other bits
 
36
 
21
37
 
22
38
#include <iostream>
23
39
#include <iomanip>
29
45
#include <vector>
30
46
#include <list>
31
47
 
 
48
#ifdef NAMD_VERSION
 
49
// use Lustre-friendly wrapper to POSIX write()
 
50
#include "fstream_namd.h"
 
51
#endif
32
52
 
33
53
class colvarparse;
34
54
class colvar;
35
55
class colvarbias;
36
56
class colvarproxy;
 
57
class colvarscript;
37
58
 
38
59
 
39
60
/// \brief Collective variables module (main class)
55
76
public:
56
77
 
57
78
  friend class colvarproxy;
 
79
  // TODO colvarscript should be unaware of colvarmodule's internals
 
80
  friend class colvarscript;
58
81
 
59
82
  /// Defining an abstract real number allows to switch precision
60
83
  typedef  double    real;
62
85
  typedef  int       residue_id;
63
86
 
64
87
  class rvector;
65
 
  template <class T,
66
 
            size_t const length> class vector1d;
67
 
  template <class T,
68
 
            size_t const outer_length,
69
 
            size_t const inner_length> class matrix2d;
 
88
  template <class T> class vector1d;
 
89
  template <class T> class matrix2d;
70
90
  class quaternion;
71
91
  class rotation;
72
92
 
85
105
  typedef std::vector<atom>::iterator       atom_iter;
86
106
  typedef std::vector<atom>::const_iterator atom_const_iter;
87
107
 
 
108
  /// Module-wide error state
 
109
  /// see constants at the top of this file
 
110
  static int errorCode;
 
111
  static inline void set_error_bits(int code)
 
112
  {
 
113
    errorCode |= code;
 
114
  }
 
115
  static inline int get_error()
 
116
  {
 
117
    return errorCode;
 
118
  }
 
119
  static inline void clear_error()
 
120
  {
 
121
    errorCode = 0;
 
122
  }
 
123
 
88
124
  /// Current step number
89
125
  static size_t it;
90
126
  /// Starting step number for this run
115
151
  /// Prefix for all output files for this run
116
152
  static std::string output_prefix;
117
153
 
118
 
  /// Prefix for files from a previous run (including restart/output)
119
 
  static std::string input_prefix;
120
 
 
121
154
  /// input restart file name (determined from input_prefix)
122
155
  static std::string restart_in_name;
123
156
 
129
162
  /// Array of named (reusable) collective variable components
130
163
  static std::vector<cvc *>     cvcs;
131
164
  /// Named cvcs register themselves at initialization time
132
 
  inline void register_cvc (cvc *p) {
 
165
  inline void register_cvc(cvc *p) {
133
166
    cvcs.push_back(p);
134
167
  }
135
168
  */
142
175
  /// \brief Number of metadynamics biases initialized (in normal
143
176
  /// conditions should be 1)
144
177
  static size_t n_meta_biases;
145
 
  /// \brief Number of harmonic biases initialized (no limit on the
 
178
  /// \brief Number of restraint biases initialized (no limit on the
146
179
  /// number)
147
 
  static size_t n_harm_biases;
 
180
  static size_t n_rest_biases;
148
181
  /// \brief Number of histograms initialized (no limit on the
149
182
  /// number)
150
183
  static size_t n_histo_biases;
158
191
 
159
192
  /// \brief Constructor \param config_name Configuration file name
160
193
  /// \param restart_name (optional) Restart file name
161
 
  colvarmodule (char const *config_name,
162
 
                colvarproxy *proxy_in);
 
194
  colvarmodule(colvarproxy *proxy);
163
195
 
164
196
  /// Destructor
165
197
  ~colvarmodule();
166
198
 
167
 
  /// Initialize collective variables
168
 
  void init_colvars (std::string const &conf);
169
 
 
170
 
  /// Initialize collective variable biases
171
 
  void init_biases (std::string const &conf);
172
 
 
173
 
  /// Re-initialize data at the beginning of a run. For use with
174
 
  /// MD codes that can change system parameters like atom masses
175
 
  /// between run commands.
176
 
  void setup();
 
199
  /// Actual function called by the destructor
 
200
  int reset();
 
201
 
 
202
  /// Open a config file, load its contents, and pass it to config_string()
 
203
  int read_config_file(char const *config_file_name);
 
204
 
 
205
  /// \brief Parse a config string assuming it is a complete configuration
 
206
  /// (i.e. calling all parse functions)
 
207
  int read_config_string(std::string const &conf);
 
208
 
 
209
  /// \brief Parse a "clean" config string (no comments)
 
210
  int parse_config(std::string &conf);
 
211
 
 
212
 
 
213
  // Parse functions (setup internal data based on a string)
 
214
 
 
215
  /// Parse the few module's global parameters
 
216
  int parse_global_params(std::string const &conf);
 
217
 
 
218
  /// Parse and initialize collective variables
 
219
  int parse_colvars(std::string const &conf);
 
220
 
 
221
  /// Parse and initialize collective variable biases
 
222
  int parse_biases(std::string const &conf);
 
223
 
 
224
  /// Parse and initialize collective variable biases of a specific type
 
225
  template <class bias_type>
 
226
  int parse_biases_type(std::string const &conf, char const *keyword, size_t &bias_count);
 
227
 
 
228
  /// Test error condition and keyword parsing
 
229
  /// on error, delete new bias
 
230
  bool check_new_bias(std::string &conf, char const *key);
 
231
 
 
232
  // "Setup" functions (change internal data based on related data
 
233
  // from the proxy that may change during program execution)
 
234
  // No additional parsing is done within these functions
 
235
 
 
236
  /// (Re)initialize internal data (currently used by LAMMPS)
 
237
  /// Also calls setup() member functions of colvars and biases
 
238
  int setup();
 
239
 
 
240
  /// (Re)initialize and (re)read the input state file calling read_restart()
 
241
  int setup_input();
 
242
 
 
243
  /// (Re)initialize the output trajectory and state file (does not write it yet)
 
244
  int setup_output();
 
245
 
 
246
#ifdef NAMD_VERSION
 
247
  typedef ofstream_namd ofstream;
 
248
#else
 
249
  typedef std::ofstream ofstream;
 
250
#endif
 
251
 
 
252
  /// Read the input restart file
 
253
  std::istream & read_restart(std::istream &is);
 
254
  /// Write the output restart file
 
255
  std::ostream & write_restart(std::ostream &os);
 
256
 
 
257
  /// Open a trajectory file if requested (and leave it open)
 
258
  int open_traj_file(std::string const &file_name);
 
259
  /// Close it
 
260
  int close_traj_file();
 
261
  /// Write in the trajectory file
 
262
  std::ostream & write_traj(std::ostream &os);
 
263
  /// Write explanatory labels in the trajectory file
 
264
  std::ostream & write_traj_label(std::ostream &os);
 
265
 
 
266
  /// Write all FINAL output files
 
267
  int write_output_files();
 
268
  /// Backup a file before writing it
 
269
  static int backup_file(char const *filename);
 
270
 
 
271
  /// Look up a bias by name; returns NULL if not found
 
272
  static colvarbias * bias_by_name(std::string const &name);
 
273
 
 
274
  /// Look up a colvar by name; returns NULL if not found
 
275
  static colvar * colvar_by_name(std::string const &name);
177
276
 
178
277
  /// Load new configuration for the given bias -
179
278
  /// currently works for harmonic (force constant and/or centers)
180
 
  void change_configuration (std::string const &bias_name, std::string const &conf);
 
279
  int change_configuration(std::string const &bias_name, std::string const &conf);
181
280
 
182
281
  /// Read a colvar value
183
282
  std::string read_colvar(std::string const &name);
184
283
 
185
284
  /// Calculate change in energy from using alt. config. for the given bias -
186
285
  /// currently works for harmonic (force constant and/or centers)
187
 
  real energy_difference (std::string const &bias_name, std::string const &conf);
 
286
  real energy_difference(std::string const &bias_name, std::string const &conf);
 
287
 
 
288
  /// Give the total number of bins for a given bias.
 
289
  int bias_bin_num(std::string const &bias_name);
 
290
  /// Calculate the bin index for a given bias.
 
291
  int bias_current_bin(std::string const &bias_name);
 
292
  //// Give the count at a given bin index.
 
293
  int bias_bin_count(std::string const &bias_name, size_t bin_index);
 
294
  //// Share among replicas.
 
295
  int bias_share(std::string const &bias_name);
188
296
 
189
297
  /// Calculate collective variables and biases
190
 
  void calc();
191
 
  /// Read the input restart file
192
 
  std::istream & read_restart (std::istream &is);
193
 
  /// Write the output restart file
194
 
  std::ostream & write_restart (std::ostream &os);
195
 
  /// Write all output files (called by the proxy)
196
 
  void write_output_files();
197
 
  /// \brief Call colvarproxy::backup_file()
198
 
  static void backup_file (char const *filename);
 
298
  int calc();
199
299
 
200
300
  /// Perform analysis
201
 
  void analyze();
 
301
  int analyze();
202
302
  /// \brief Read a collective variable trajectory (post-processing
203
303
  /// only, not called at runtime)
204
 
  bool read_traj (char const *traj_filename,
 
304
  int read_traj(char const *traj_filename,
205
305
                  size_t      traj_read_begin,
206
306
                  size_t      traj_read_end);
207
307
 
208
 
  /// Get the pointer of a colvar from its name (returns NULL if not found)
209
 
  static colvar * colvar_p (std::string const &name);
210
 
 
211
308
  /// Quick conversion of an object to a string
212
 
  template<typename T> static std::string to_str (T const &x,
 
309
  template<typename T> static std::string to_str(T const &x,
213
310
                                                  size_t const &width = 0,
214
311
                                                  size_t const &prec = 0);
215
312
  /// Quick conversion of a vector of objects to a string
216
 
  template<typename T> static std::string to_str (std::vector<T> const &x,
 
313
  template<typename T> static std::string to_str(std::vector<T> const &x,
217
314
                                                  size_t const &width = 0,
218
315
                                                  size_t const &prec = 0);
219
316
 
220
317
  /// Reduce the number of characters in a string
221
 
  static inline std::string wrap_string (std::string const &s,
 
318
  static inline std::string wrap_string(std::string const &s,
222
319
                                         size_t const &nchars)
223
320
  {
224
321
    if (!s.size())
225
 
      return std::string (nchars, ' ');
 
322
      return std::string(nchars, ' ');
226
323
    else
227
 
      return ( (s.size() <= size_t (nchars)) ?
228
 
               (s+std::string (nchars-s.size(), ' ')) :
229
 
               (std::string (s, 0, nchars)) );
 
324
      return ( (s.size() <= size_t(nchars)) ?
 
325
               (s+std::string(nchars-s.size(), ' ')) :
 
326
               (std::string(s, 0, nchars)) );
230
327
  }
231
328
 
232
329
  /// Number of characters to represent a time step
262
359
  static void request_system_force();
263
360
 
264
361
  /// Print a message to the main log
265
 
  static void log (std::string const &message);
 
362
  static void log(std::string const &message);
266
363
 
267
364
  /// Print a message to the main log and exit with error code
268
 
  static void fatal_error (std::string const &message);
 
365
  static void fatal_error(std::string const &message);
 
366
 
 
367
  /// Print a message to the main log and set global error code
 
368
  static void error(std::string const &message, int code = GENERAL_ERROR);
269
369
 
270
370
  /// Print a message to the main log and exit normally
271
 
  static void exit (std::string const &message);
 
371
  static void exit(std::string const &message);
272
372
 
 
373
  // Replica exchange commands.
 
374
  static bool replica_enabled();
 
375
  static int replica_index();
 
376
  static int replica_num();
 
377
  static void replica_comm_barrier();
 
378
  static int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
 
379
  static int replica_comm_send(char* msg_data, int msg_len, int dest_rep);
273
380
 
274
381
  /// \brief Get the distance between two atomic positions with pbcs handled
275
382
  /// correctly
276
 
  static rvector position_distance (atom_pos const &pos1,
 
383
  static rvector position_distance(atom_pos const &pos1,
277
384
                                    atom_pos const &pos2);
278
385
 
279
386
 
283
390
  /// Note: in the case of periodic boundary conditions, this provides
284
391
  /// an analytical square distance (while taking the square of
285
392
  /// position_distance() would produce leads to a cusp)
286
 
  static real position_dist2 (atom_pos const &pos1,
 
393
  static real position_dist2(atom_pos const &pos1,
287
394
                              atom_pos const &pos2);
288
395
 
289
396
  /// \brief Get the closest periodic image to a reference position
290
397
  /// \param pos The position to look for the closest periodic image
291
398
  /// \param ref_pos (optional) The reference position
292
 
  static void select_closest_image (atom_pos &pos,
 
399
  static void select_closest_image(atom_pos &pos,
293
400
                                    atom_pos const &ref_pos);
294
401
 
295
402
  /// \brief Perform select_closest_image() on a set of atomic positions
296
403
  ///
297
404
  /// After that, distance vectors can then be calculated directly,
298
405
  /// without using position_distance()
299
 
  static void select_closest_images (std::vector<atom_pos> &pos,
 
406
  static void select_closest_images(std::vector<atom_pos> &pos,
300
407
                                     atom_pos const &ref_pos);
301
408
 
302
409
 
307
414
  static std::list<std::vector<int> > index_groups;
308
415
 
309
416
  /// \brief Read a Gromacs .ndx file
310
 
  static void read_index_file (char const *filename);
 
417
  static int read_index_file(char const *filename);
311
418
 
312
419
 
313
420
  /// \brief Create atoms from a file \param filename name of the file
314
421
  /// (usually a PDB) \param atoms array of the atoms to be allocated
315
422
  /// \param pdb_field (optiona) if "filename" is a PDB file, use this
316
423
  /// field to determine which are the atoms to be set
317
 
  static void load_atoms (char const *filename,
 
424
  static int load_atoms(char const *filename,
318
425
                          std::vector<atom> &atoms,
319
426
                          std::string const &pdb_field,
320
427
                          double const pdb_field_value = 0.0);
321
428
 
322
429
  /// \brief Load the coordinates for a group of atoms from a file
323
430
  /// (PDB or XYZ)
324
 
  static void load_coords (char const *filename,
 
431
  static int load_coords(char const *filename,
325
432
                           std::vector<atom_pos> &pos,
326
433
                           const std::vector<int> &indices,
327
434
                           std::string const &pdb_field,
329
436
 
330
437
  /// \brief Load the coordinates for a group of atoms from an
331
438
  /// XYZ file
332
 
  static void load_coords_xyz (char const *filename,
 
439
  static int load_coords_xyz(char const *filename,
333
440
                              std::vector<atom_pos> &pos,
334
441
                              const std::vector<int> &indices);
335
442
 
345
452
  std::string   restart_out_name;
346
453
 
347
454
  /// Pseudo-random number with Gaussian distribution
348
 
  static real rand_gaussian (void);
 
455
  static real rand_gaussian(void);
 
456
 
349
457
protected:
350
458
 
351
459
  /// Configuration file
355
463
  colvarparse *parse;
356
464
 
357
465
  /// Name of the trajectory file
358
 
  std::string   cv_traj_name;
 
466
  std::string cv_traj_name;
359
467
 
360
468
  /// Collective variables output trajectory file
361
 
  std::ofstream cv_traj_os;
 
469
  colvarmodule::ofstream cv_traj_os;
362
470
 
363
471
  /// Appending to the existing trajectory file?
364
 
  bool          cv_traj_append;
 
472
  bool cv_traj_append;
365
473
 
366
474
  /// Output restart file
367
 
  std::ofstream restart_out_os;
 
475
  colvarmodule::ofstream restart_out_os;
 
476
 
 
477
  /// \brief Counter for the current depth in the object hierarchy (useg e.g. in output
 
478
  static size_t depth;
 
479
 
 
480
  /// Use scripted colvars forces?
 
481
  static bool use_scripted_forces;
 
482
 
 
483
public:
368
484
 
369
485
  /// \brief Pointer to the proxy object, used to retrieve atomic data
370
486
  /// from the hosting program; it is static in order to be accessible
371
487
  /// from static functions in the colvarmodule class
372
488
  static colvarproxy *proxy;
373
489
 
374
 
  /// \brief Counter for the current depth in the object hierarchy (useg e.g. in outpu
375
 
  static size_t depth;
376
 
 
377
 
public:
378
 
 
379
490
  /// Increase the depth (number of indentations in the output)
380
491
  static void increase_depth();
381
492
 
382
493
  /// Decrease the depth (number of indentations in the output)
383
494
  static void decrease_depth();
 
495
 
 
496
  static inline bool scripted_forces() { return use_scripted_forces; }
384
497
};
385
498
 
386
499
 
395
508
std::istream & operator >> (std::istream &is, cvm::rvector &v);
396
509
 
397
510
 
398
 
template<typename T> std::string cvm::to_str (T const &x,
 
511
template<typename T> std::string cvm::to_str(T const &x,
399
512
                                              size_t const &width,
400
513
                                              size_t const &prec) {
401
514
  std::ostringstream os;
402
 
  if (width) os.width (width);
 
515
  if (width) os.width(width);
403
516
  if (prec) {
404
 
    os.setf (std::ios::scientific, std::ios::floatfield);
405
 
    os.precision (prec);
 
517
    os.setf(std::ios::scientific, std::ios::floatfield);
 
518
    os.precision(prec);
406
519
  }
407
520
  os << x;
408
521
  return os.str();
409
522
}
410
523
 
411
 
template<typename T> std::string cvm::to_str (std::vector<T> const &x,
 
524
template<typename T> std::string cvm::to_str(std::vector<T> const &x,
412
525
                                              size_t const &width,
413
526
                                              size_t const &prec) {
414
 
  if (!x.size()) return std::string ("");
 
527
  if (!x.size()) return std::string("");
415
528
  std::ostringstream os;
416
529
  if (prec) {
417
 
    os.setf (std::ios::scientific, std::ios::floatfield);
 
530
    os.setf(std::ios::scientific, std::ios::floatfield);
418
531
  }
419
532
  os << "{ ";
420
 
  if (width) os.width (width);
421
 
  if (prec) os.precision (prec);
 
533
  if (width) os.width(width);
 
534
  if (prec) os.precision(prec);
422
535
  os << x[0];
423
536
  for (size_t i = 1; i < x.size(); i++) {
424
537
    os << ", ";
425
 
    if (width) os.width (width);
426
 
    if (prec) os.precision (prec);
 
538
    if (width) os.width(width);
 
539
    if (prec) os.precision(prec);
427
540
    os << x[i];
428
541
  }
429
542
  os << " }";
454
567
  return proxy->dt();
455
568
}
456
569
 
 
570
// Replica exchange commands
 
571
inline bool cvm::replica_enabled() {
 
572
  return proxy->replica_enabled();
 
573
}
 
574
inline int cvm::replica_index() {
 
575
  return proxy->replica_index();
 
576
}
 
577
inline int cvm::replica_num() {
 
578
  return proxy->replica_num();
 
579
}
 
580
inline void cvm::replica_comm_barrier() {
 
581
  return proxy->replica_comm_barrier();
 
582
}
 
583
inline int cvm::replica_comm_recv(char* msg_data, int buf_len, int src_rep) {
 
584
  return proxy->replica_comm_recv(msg_data,buf_len,src_rep);
 
585
}
 
586
inline int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep) {
 
587
  return proxy->replica_comm_send(msg_data,msg_len,dest_rep);
 
588
}
 
589
 
 
590
 
457
591
inline void cvm::request_system_force()
458
592
{
459
 
  proxy->request_system_force (true);
 
593
  proxy->request_system_force(true);
460
594
}
461
595
 
462
 
inline void cvm::select_closest_image (atom_pos &pos,
 
596
inline void cvm::select_closest_image(atom_pos &pos,
463
597
                                       atom_pos const &ref_pos)
464
598
{
465
 
  proxy->select_closest_image (pos, ref_pos);
 
599
  proxy->select_closest_image(pos, ref_pos);
466
600
}
467
601
 
468
 
inline void cvm::select_closest_images (std::vector<atom_pos> &pos,
 
602
inline void cvm::select_closest_images(std::vector<atom_pos> &pos,
469
603
                                        atom_pos const &ref_pos)
470
604
{
471
 
  proxy->select_closest_images (pos, ref_pos);
 
605
  proxy->select_closest_images(pos, ref_pos);
472
606
}
473
607
 
474
 
inline cvm::rvector cvm::position_distance (atom_pos const &pos1,
 
608
inline cvm::rvector cvm::position_distance(atom_pos const &pos1,
475
609
                                            atom_pos const &pos2)
476
610
{
477
 
  return proxy->position_distance (pos1, pos2);
 
611
  return proxy->position_distance(pos1, pos2);
478
612
}
479
613
 
480
 
inline cvm::real cvm::position_dist2 (cvm::atom_pos const &pos1,
 
614
inline cvm::real cvm::position_dist2(cvm::atom_pos const &pos1,
481
615
                                      cvm::atom_pos const &pos2)
482
616
{
483
 
  return proxy->position_dist2 (pos1, pos2);
 
617
  return proxy->position_dist2(pos1, pos2);
484
618
}
485
619
 
486
 
inline cvm::real cvm::rand_gaussian (void)
 
620
inline cvm::real cvm::rand_gaussian(void)
487
621
{
488
622
  return proxy->rand_gaussian();
489
623
}
490
624
 
491
625
#endif
492
 
 
493
 
 
494
 
// Emacs
495
 
// Local Variables:
496
 
// mode: C++
497
 
// End: