← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/colvars/colvaratoms.cpp

  • 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/colvaratoms.cpp
 
1
/// -*- c++ -*-
 
2
 
1
3
#include "colvarmodule.h"
2
4
#include "colvarparse.h"
3
5
#include "colvaratoms.h"
11
13
 
12
14
// Note: "conf" is the configuration of the cvc who is using this atom group;
13
15
// "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...)
14
 
cvm::atom_group::atom_group (std::string const &conf,
 
16
cvm::atom_group::atom_group(std::string const &conf,
15
17
                             char const        *key)
16
 
  : b_center (false), b_rotate (false), b_user_defined_fit (false),
17
 
    b_fit_gradients (false),
18
 
    ref_pos_group (NULL),
19
 
    noforce (false)
 
18
  : b_center(false), b_rotate(false), b_user_defined_fit(false),
 
19
    b_fit_gradients(false),
 
20
    ref_pos_group(NULL),
 
21
    noforce(false)
20
22
{
21
 
  cvm::log ("Defining atom group \""+
22
 
            std::string (key)+"\".\n");
 
23
  cvm::log("Defining atom group \""+
 
24
            std::string(key)+"\".\n");
23
25
  // real work is done by parse
24
 
  parse (conf, key);
 
26
  parse(conf, key);
25
27
}
26
28
 
27
29
 
28
 
cvm::atom_group::atom_group (std::vector<cvm::atom> const &atoms)
29
 
  : b_dummy (false), b_center (false), b_rotate (false),
30
 
    b_fit_gradients (false), ref_pos_group (NULL),
31
 
    noforce (false)
 
30
cvm::atom_group::atom_group(std::vector<cvm::atom> const &atoms)
 
31
  : b_dummy(false), b_center(false), b_rotate(false),
 
32
    b_fit_gradients(false), ref_pos_group(NULL),
 
33
    noforce(false)
32
34
{
33
 
  this->reserve (atoms.size());
 
35
  this->reserve(atoms.size());
34
36
  for (size_t i = 0; i < atoms.size(); i++) {
35
 
    this->push_back (atoms[i]);
 
37
    this->push_back(atoms[i]);
36
38
  }
37
39
  total_mass = 0.0;
38
40
  for (cvm::atom_iter ai = this->begin();
43
45
 
44
46
 
45
47
cvm::atom_group::atom_group()
46
 
  : b_dummy (false), b_center (false), b_rotate (false),
47
 
    b_fit_gradients (false), ref_pos_group (NULL),
48
 
    noforce (false)
 
48
  : b_dummy(false), b_center(false), b_rotate(false),
 
49
    b_user_defined_fit(false), b_fit_gradients(false),
 
50
    ref_pos_group(NULL), noforce(false)
49
51
{
50
52
  total_mass = 0.0;
51
53
}
60
62
}
61
63
 
62
64
 
63
 
void cvm::atom_group::add_atom (cvm::atom const &a)
 
65
void cvm::atom_group::add_atom(cvm::atom const &a)
64
66
{
65
67
  if (b_dummy) {
66
 
    cvm::fatal_error ("Error: cannot add atoms to a dummy group.\n");
 
68
    cvm::error("Error: cannot add atoms to a dummy group.\n", INPUT_ERROR);
67
69
  } else {
68
 
    this->push_back (a);
 
70
    this->push_back(a);
69
71
    total_mass += a.mass;
70
72
  }
71
73
}
78
80
       ai != this->end(); ai++) {
79
81
    total_mass += ai->mass;
80
82
  }
81
 
  cvm::log ("Re-initialized atom group "+name+":"+cvm::to_str (i)+"/"+
82
 
            cvm::to_str (j)+". "+ cvm::to_str (this->size())+
83
 
            " atoms: total mass = "+cvm::to_str (this->total_mass)+".\n");
 
83
  cvm::log("Re-initialized atom group "+name+":"+cvm::to_str(i)+"/"+
 
84
            cvm::to_str(j)+". "+ cvm::to_str(this->size())+
 
85
            " atoms: total mass = "+cvm::to_str(this->total_mass)+".\n");
84
86
}
85
87
 
86
 
void cvm::atom_group::parse (std::string const &conf,
 
88
int cvm::atom_group::parse(std::string const &conf,
87
89
                             char const        *key)
88
90
{
89
91
  std::string group_conf;
92
94
  // not the config string of this group, but of its parent object
93
95
  // (which has already taken care of the delimiters)
94
96
  save_delimiters = false;
95
 
  key_lookup (conf, key, group_conf, dummy_pos);
 
97
  key_lookup(conf, key, group_conf, dummy_pos);
96
98
  // restoring the normal value, because we do want keywords checked
97
99
  // inside "group_conf"
98
100
  save_delimiters = true;
99
101
 
100
102
  if (group_conf.size() == 0) {
101
 
    cvm::fatal_error ("Error: atom group \""+
102
 
                      std::string (key)+"\" is set, but "
103
 
                      "has no definition.\n");
 
103
    cvm::error("Error: atom group \""+std::string(key)+
 
104
                "\" is set, but has no definition.\n",
 
105
                INPUT_ERROR);
 
106
    return COLVARS_ERROR;
104
107
  }
105
108
 
106
109
  cvm::increase_depth();
107
110
 
108
 
  cvm::log ("Initializing atom group \""+std::string (key)+"\".\n");
 
111
  cvm::log("Initializing atom group \""+std::string(key)+"\".\n");
109
112
 
110
113
  // whether or not to include messages in the log
111
114
  // colvarparse::Parse_Mode mode = parse_silent;
121
124
    std::string numbers_conf = "";
122
125
    size_t pos = 0;
123
126
    std::vector<int> atom_indexes;
124
 
    while (key_lookup (group_conf, "atomNumbers", numbers_conf, pos)) {
 
127
    while (key_lookup(group_conf, "atomNumbers", numbers_conf, pos)) {
125
128
      if (numbers_conf.size()) {
126
 
        std::istringstream is (numbers_conf);
 
129
        std::istringstream is(numbers_conf);
127
130
        int ia;
128
131
        while (is >> ia) {
129
 
          atom_indexes.push_back (ia);
 
132
          atom_indexes.push_back(ia);
130
133
        }
131
134
      }
132
135
 
133
136
      if (atom_indexes.size()) {
134
 
        this->reserve (this->size()+atom_indexes.size());
 
137
        this->reserve(this->size()+atom_indexes.size());
135
138
        for (size_t i = 0; i < atom_indexes.size(); i++) {
136
 
          this->push_back (cvm::atom (atom_indexes[i]));
 
139
          this->push_back(cvm::atom(atom_indexes[i]));
137
140
        }
138
 
      } else
139
 
        cvm::fatal_error ("Error: no numbers provided for \""
140
 
                          "atomNumbers\".\n");
 
141
        if (cvm::get_error()) return COLVARS_ERROR;
 
142
      } else {
 
143
        cvm::error("Error: no numbers provided for \""
 
144
                    "atomNumbers\".\n", INPUT_ERROR);
 
145
        return COLVARS_ERROR;
 
146
      }
141
147
 
142
148
      atom_indexes.clear();
143
149
    }
144
150
 
145
151
    std::string index_group_name;
146
 
    if (get_keyval (group_conf, "indexGroup", index_group_name)) {
 
152
    if (get_keyval(group_conf, "indexGroup", index_group_name)) {
147
153
      // use an index group from the index file read globally
148
154
      std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
149
155
      std::list<std::vector<int> >::iterator index_groups_i = cvm::index_groups.begin();
150
 
      for ( ; names_i != cvm::index_group_names.end() ; names_i++, index_groups_i++) {
 
156
      for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) {
151
157
        if (*names_i == index_group_name)
152
158
          break;
153
159
      }
154
160
      if (names_i == cvm::index_group_names.end()) {
155
 
        cvm::fatal_error ("Error: could not find index group "+
156
 
                          index_group_name+" among those provided by the index file.\n");
 
161
        cvm::error("Error: could not find index group "+
 
162
                    index_group_name+" among those provided by the index file.\n",
 
163
                    INPUT_ERROR);
 
164
        return COLVARS_ERROR;
157
165
      }
158
 
      this->reserve (index_groups_i->size());
 
166
      this->reserve(index_groups_i->size());
159
167
      for (size_t i = 0; i < index_groups_i->size(); i++) {
160
 
        this->push_back (cvm::atom ((*index_groups_i)[i]));
 
168
        this->push_back(cvm::atom((*index_groups_i)[i]));
161
169
      }
 
170
      if (cvm::get_error()) return COLVARS_ERROR;
162
171
    }
163
172
  }
164
173
 
165
174
  {
166
175
    std::string range_conf = "";
167
176
    size_t pos = 0;
168
 
    while (key_lookup (group_conf, "atomNumbersRange",
 
177
    while (key_lookup(group_conf, "atomNumbersRange",
169
178
                       range_conf, pos)) {
170
179
 
171
180
      if (range_conf.size()) {
172
 
        std::istringstream is (range_conf);
 
181
        std::istringstream is(range_conf);
173
182
        int initial, final;
174
183
        char dash;
175
184
        if ( (is >> initial) && (initial > 0) &&
176
185
             (is >> dash) && (dash == '-') &&
177
186
             (is >> final) && (final > 0) ) {
178
187
          for (int anum = initial; anum <= final; anum++) {
179
 
            this->push_back (cvm::atom (anum));
 
188
            this->push_back(cvm::atom(anum));
180
189
          }
 
190
          if (cvm::get_error()) return COLVARS_ERROR;
181
191
          range_conf = "";
182
192
          continue;
183
193
        }
184
194
      }
185
195
 
186
 
      cvm::fatal_error ("Error: no valid definition for \""
187
 
                        "atomNumbersRange\", \""+
188
 
                        range_conf+"\".\n");
 
196
      cvm::error("Error: no valid definition for \"atomNumbersRange\", \""+
 
197
                    range_conf+"\".\n", INPUT_ERROR);
189
198
    }
190
199
  }
191
200
 
192
201
  std::vector<std::string> psf_segids;
193
 
  get_keyval (group_conf, "psfSegID", psf_segids, std::vector<std::string> (), mode);
 
202
  get_keyval(group_conf, "psfSegID", psf_segids, std::vector<std::string> (), mode);
194
203
  for (std::vector<std::string>::iterator psii = psf_segids.begin();
195
 
       psii < psf_segids.end(); psii++) {
 
204
       psii < psf_segids.end(); ++psii) {
196
205
 
197
206
    if ( (psii->size() == 0) || (psii->size() > 4) ) {
198
 
      cvm::fatal_error ("Error: invalid segmend identifier provided, \""+
199
 
                        (*psii)+"\".\n");
 
207
      cvm::error("Error: invalid segmend identifier provided, \""+
 
208
                  (*psii)+"\".\n", INPUT_ERROR);
200
209
    }
201
210
  }
202
211
 
205
214
    size_t pos = 0;
206
215
    size_t range_count = 0;
207
216
    std::vector<std::string>::iterator psii = psf_segids.begin();
208
 
    while (key_lookup (group_conf, "atomNameResidueRange",
 
217
    while (key_lookup(group_conf, "atomNameResidueRange",
209
218
                       range_conf, pos)) {
210
219
      range_count++;
211
220
 
212
221
      if (range_count > psf_segids.size()) {
213
 
        cvm::fatal_error ("Error: more instances of \"atomNameResidueRange\" than "
214
 
                          "values of \"psfSegID\".\n");
 
222
        cvm::error("Error: more instances of \"atomNameResidueRange\" than "
 
223
                    "values of \"psfSegID\".\n", INPUT_ERROR);
215
224
      }
216
225
 
217
 
      std::string const &psf_segid = psf_segids.size() ? *psii : std::string ("");
 
226
      std::string const &psf_segid = psf_segids.size() ? *psii : std::string("");
218
227
 
219
228
      if (range_conf.size()) {
220
229
 
221
 
        std::istringstream is (range_conf);
 
230
        std::istringstream is(range_conf);
222
231
        std::string atom_name;
223
232
        int initial, final;
224
233
        char dash;
227
236
             (is >> dash) && (dash == '-') &&
228
237
             (is >> final) && (final > 0) ) {
229
238
          for (int resid = initial; resid <= final; resid++) {
230
 
            this->push_back (cvm::atom (resid, atom_name, psf_segid));
 
239
            this->push_back(cvm::atom(resid, atom_name, psf_segid));
231
240
          }
 
241
          if (cvm::get_error()) return COLVARS_ERROR;
232
242
          range_conf = "";
233
243
        } else {
234
 
          cvm::fatal_error ("Error: cannot parse definition for \""
235
 
                            "atomNameResidueRange\", \""+
236
 
                            range_conf+"\".\n");
 
244
          cvm::error("Error: cannot parse definition for \""
 
245
                      "atomNameResidueRange\", \""+
 
246
                      range_conf+"\".\n");
237
247
        }
238
248
 
239
249
      } else {
240
 
        cvm::fatal_error ("Error: atomNameResidueRange with empty definition.\n");
 
250
        cvm::error("Error: atomNameResidueRange with empty definition.\n");
241
251
      }
242
252
 
243
253
      if (psf_segid.size())
244
 
        psii++;
 
254
        ++psii;
245
255
    }
246
256
  }
247
257
 
248
258
  {
249
259
    // read the atoms from a file
250
260
    std::string atoms_file_name;
251
 
    if (get_keyval (group_conf, "atomsFile", atoms_file_name, std::string (""), mode)) {
 
261
    if (get_keyval(group_conf, "atomsFile", atoms_file_name, std::string(""), mode)) {
252
262
 
253
263
      std::string atoms_col;
254
 
      if (!get_keyval (group_conf, "atomsCol", atoms_col, std::string (""), mode)) {
255
 
        cvm::fatal_error ("Error: parameter atomsCol is required if atomsFile is set.\n");
 
264
      if (!get_keyval(group_conf, "atomsCol", atoms_col, std::string(""), mode)) {
 
265
        cvm::error("Error: parameter atomsCol is required if atomsFile is set.\n",
 
266
                    INPUT_ERROR);
256
267
      }
257
268
 
258
269
      double atoms_col_value;
259
 
      bool const atoms_col_value_defined = get_keyval (group_conf, "atomsColValue", atoms_col_value, 0.0, mode);
260
 
      if (atoms_col_value_defined && (!atoms_col_value))
261
 
        cvm::fatal_error ("Error: atomsColValue, "
262
 
                          "if provided, must be non-zero.\n");
 
270
      bool const atoms_col_value_defined = get_keyval(group_conf, "atomsColValue", atoms_col_value, 0.0, mode);
 
271
      if (atoms_col_value_defined && (!atoms_col_value)) {
 
272
        cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR);
 
273
      }
263
274
 
264
 
      cvm::load_atoms (atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
 
275
      cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
265
276
    }
266
277
  }
267
278
 
 
279
  // Catch any errors from all the initialization steps above
 
280
  if (cvm::get_error()) return COLVARS_ERROR;
 
281
 
268
282
  for (std::vector<cvm::atom>::iterator a1 = this->begin();
269
 
       a1 != this->end(); a1++) {
 
283
       a1 != this->end(); ++a1) {
270
284
    std::vector<cvm::atom>::iterator a2 = a1;
271
285
    ++a2;
272
 
    for ( ; a2 != this->end(); a2++) {
 
286
    for ( ; a2 != this->end(); ++a2) {
273
287
      if (a1->id == a2->id) {
274
288
        if (cvm::debug())
275
 
          cvm::log ("Discarding doubly counted atom with number "+
276
 
                    cvm::to_str (a1->id+1)+".\n");
277
 
        a2 = this->erase (a2);
 
289
          cvm::log("Discarding doubly counted atom with number "+
 
290
                    cvm::to_str(a1->id+1)+".\n");
 
291
        a2 = this->erase(a2);
278
292
        if (a2 == this->end())
279
293
          break;
280
294
      }
281
295
    }
282
296
  }
283
297
 
284
 
  if (get_keyval (group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) {
 
298
  if (get_keyval(group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) {
285
299
    b_dummy = true;
286
300
    this->total_mass = 1.0;
287
301
  } else
288
302
    b_dummy = false;
289
303
 
290
 
  if (b_dummy && (this->size()))
291
 
    cvm::fatal_error ("Error: cannot set up group \""+
292
 
                      std::string (key)+"\" as a dummy atom "
293
 
                      "and provide it with atom definitions.\n");
 
304
  if (b_dummy && (this->size())) {
 
305
    cvm::error("Error: cannot set up group \""+
 
306
                std::string(key)+"\" as a dummy atom "
 
307
                "and provide it with atom definitions.\n", INPUT_ERROR);
 
308
  }
294
309
 
295
 
#if (! defined (COLVARS_STANDALONE))
 
310
#if(! defined(COLVARS_STANDALONE))
296
311
  if ( (!b_dummy) && (!cvm::b_analysis) && (!(this->size())) ) {
297
 
    cvm::fatal_error ("Error: no atoms defined for atom group \""+
298
 
                      std::string (key)+"\".\n");
 
312
    cvm::error("Error: no atoms defined for atom group \""+
 
313
                      std::string(key)+"\".\n");
299
314
  }
300
315
#endif
301
316
 
302
317
  if (!b_dummy) {
 
318
 
 
319
    // calculate total mass (TODO: this is the step that most needs deferred re-initialization)
303
320
    this->total_mass = 0.0;
304
321
    for (cvm::atom_iter ai = this->begin();
305
322
         ai != this->end(); ai++) {
306
323
      this->total_mass += ai->mass;
307
324
    }
308
 
  }
309
325
 
310
 
  if (!b_dummy) {
 
326
    // whether these atoms will ever receive forces or not
311
327
    bool enable_forces = true;
312
328
    // disableForces is deprecated
313
 
    if (get_keyval (group_conf, "enableForces", enable_forces, true, mode)) {
 
329
    if (get_keyval(group_conf, "enableForces", enable_forces, true, mode)) {
314
330
      noforce = !enable_forces;
315
331
    } else {
316
 
      get_keyval (group_conf, "disableForces", noforce, false, mode);
 
332
      get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent);
317
333
    }
 
334
 
 
335
    get_keyval(group_conf, "weights", weights, weights, colvarparse::parse_silent);
318
336
  }
319
337
 
320
338
  // FITTING OPTIONS
321
339
 
322
 
  bool b_defined_center = get_keyval (group_conf, "centerReference", b_center, false, mode);
323
 
  bool b_defined_rotate = get_keyval (group_conf, "rotateReference", b_rotate, false, mode);
 
340
  bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false, mode);
 
341
  bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false, mode);
324
342
  // is the user setting explicit options?
325
343
  b_user_defined_fit = b_defined_center || b_defined_rotate;
326
344
 
327
 
  get_keyval (group_conf, "enableFitGradients", b_fit_gradients, true, mode);
 
345
  get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true, mode);
328
346
 
329
347
  if (b_center || b_rotate) {
330
348
 
331
349
    if (b_dummy)
332
 
      cvm::fatal_error ("Error: centerReference or rotateReference "
333
 
                        "cannot be defined for a dummy atom.\n");
 
350
      cvm::error("Error: centerReference or rotateReference "
 
351
                  "cannot be defined for a dummy atom.\n");
334
352
 
335
 
    if (key_lookup (group_conf, "refPositionsGroup")) {
 
353
    if (key_lookup(group_conf, "refPositionsGroup")) {
336
354
      // instead of this group, define another group to compute the fit
337
355
      if (ref_pos_group) {
338
 
        cvm::fatal_error ("Error: the atom group \""+
339
 
                          std::string (key)+"\" has already a reference group "
340
 
                          "for the rototranslational fit, which was communicated by the "
341
 
                          "colvar component.  You should not use refPositionsGroup "
342
 
                          "in this case.\n");
 
356
        cvm::error("Error: the atom group \""+
 
357
                    std::string(key)+"\" has already a reference group "
 
358
                    "for the rototranslational fit, which was communicated by the "
 
359
                    "colvar component.  You should not use refPositionsGroup "
 
360
                    "in this case.\n");
343
361
      }
344
 
      cvm::log ("Within atom group \""+std::string (key)+"\":\n");
345
 
      ref_pos_group = new atom_group (group_conf, "refPositionsGroup");
 
362
      cvm::log("Within atom group \""+std::string(key)+"\":\n");
 
363
      ref_pos_group = new atom_group(group_conf, "refPositionsGroup");
 
364
 
 
365
      // regardless of the configuration, fit gradients must be calculated by refPositionsGroup
 
366
      ref_pos_group->b_fit_gradients = this->b_fit_gradients;
 
367
      this->b_fit_gradients = false;
346
368
    }
347
369
 
348
370
    atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
349
371
 
350
 
    get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode);
 
372
    get_keyval(group_conf, "refPositions", ref_pos, ref_pos, mode);
351
373
 
352
374
    std::string ref_pos_file;
353
 
    if (get_keyval (group_conf, "refPositionsFile", ref_pos_file, std::string (""), mode)) {
 
375
    if (get_keyval(group_conf, "refPositionsFile", ref_pos_file, std::string(""), mode)) {
354
376
 
355
377
      if (ref_pos.size()) {
356
 
        cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and "
357
 
                          "\"refPositions\" at the same time.\n");
 
378
        cvm::error("Error: cannot specify \"refPositionsFile\" and "
 
379
                    "\"refPositions\" at the same time.\n");
358
380
      }
359
381
 
360
382
      std::string ref_pos_col;
361
383
      double ref_pos_col_value;
362
384
 
363
 
      if (get_keyval (group_conf, "refPositionsCol", ref_pos_col, std::string (""), mode)) {
 
385
      if (get_keyval(group_conf, "refPositionsCol", ref_pos_col, std::string(""), mode)) {
364
386
        // if provided, use PDB column to select coordinates
365
 
        bool found = get_keyval (group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode);
 
387
        bool found = get_keyval(group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode);
366
388
        if (found && !ref_pos_col_value)
367
 
          cvm::fatal_error ("Error: refPositionsColValue, "
368
 
                            "if provided, must be non-zero.\n");
 
389
          cvm::error("Error: refPositionsColValue, "
 
390
                      "if provided, must be non-zero.\n");
369
391
      } else {
370
392
        // if not, rely on existing atom indices for the group
371
393
        group_for_fit->create_sorted_ids();
372
 
        ref_pos.resize (group_for_fit->size());
 
394
        ref_pos.resize(group_for_fit->size());
373
395
      }
374
396
 
375
 
      cvm::load_coords (ref_pos_file.c_str(), ref_pos, group_for_fit->sorted_ids,
 
397
      cvm::load_coords(ref_pos_file.c_str(), ref_pos, group_for_fit->sorted_ids,
376
398
                        ref_pos_col, ref_pos_col_value);
377
399
    }
378
400
 
380
402
 
381
403
      if (b_rotate) {
382
404
        if (ref_pos.size() != group_for_fit->size())
383
 
          cvm::fatal_error ("Error: the number of reference positions provided ("+
384
 
                            cvm::to_str (ref_pos.size())+
385
 
                            ") does not match the number of atoms within \""+
386
 
                            std::string (key)+
387
 
                            "\" ("+cvm::to_str (group_for_fit->size())+
388
 
                            "): to perform a rotational fit, "+
389
 
                            "these numbers should be equal.\n");
 
405
          cvm::error("Error: the number of reference positions provided("+
 
406
                      cvm::to_str(ref_pos.size())+
 
407
                      ") does not match the number of atoms within \""+
 
408
                      std::string(key)+
 
409
                      "\" ("+cvm::to_str(group_for_fit->size())+
 
410
                      "): to perform a rotational fit, "+
 
411
                      "these numbers should be equal.\n", INPUT_ERROR);
390
412
      }
391
413
 
392
414
      // save the center of geometry of ref_pos and subtract it
393
415
      center_ref_pos();
394
416
 
395
417
    } else {
396
 
#if (! defined (COLVARS_STANDALONE))
397
 
      cvm::fatal_error ("Error: no reference positions provided.\n");
 
418
#if(! defined(COLVARS_STANDALONE))
 
419
      cvm::error("Error: no reference positions provided.\n");
398
420
#endif
399
421
    }
400
422
 
401
423
    if (b_fit_gradients) {
402
 
      group_for_fit->fit_gradients.assign (group_for_fit->size(), cvm::atom_pos (0.0, 0.0, 0.0));
403
 
      rot.request_group1_gradients (group_for_fit->size());
 
424
      group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::atom_pos(0.0, 0.0, 0.0));
 
425
      rot.request_group1_gradients(group_for_fit->size());
404
426
    }
405
427
 
406
428
    if (b_rotate && !noforce) {
407
 
      cvm::log ("Warning: atom group \""+std::string (key)+
 
429
      cvm::log("Warning: atom group \""+std::string(key)+
408
430
                "\" will be aligned to a fixed orientation given by the reference positions provided.  "
409
431
                "If the internal structure of the group changes too much (i.e. its RMSD is comparable "
410
432
                "to its radius of gyration), the optimal rotation and its gradients may become discontinuous.  "
411
433
                "If that happens, use refPositionsGroup (or a different definition for it if already defined) "
412
434
                "to align the coordinates.\n");
413
435
      // initialize rot member data
414
 
      rot.request_group1_gradients (this->size());
 
436
      rot.request_group1_gradients(this->size());
415
437
    }
416
438
 
417
439
  }
418
440
 
419
441
  if (cvm::debug())
420
 
    cvm::log ("Done initializing atom group with name \""+
421
 
              std::string (key)+"\".\n");
422
 
 
423
 
  this->check_keywords (group_conf, key);
424
 
 
425
 
  cvm::log ("Atom group \""+std::string (key)+"\" defined, "+
426
 
            cvm::to_str (this->size())+" atoms initialized: total mass = "+
427
 
            cvm::to_str (this->total_mass)+".\n");
 
442
    cvm::log("Done initializing atom group with name \""+
 
443
              std::string(key)+"\".\n");
 
444
 
 
445
  this->check_keywords(group_conf, key);
 
446
  if (cvm::get_error()) {
 
447
    cvm::error("Error setting up atom group \""+std::string(key)+"\".");
 
448
    return COLVARS_ERROR;
 
449
  }
 
450
 
 
451
  cvm::log("Atom group \""+std::string(key)+"\" defined, "+
 
452
            cvm::to_str(this->size())+" atoms initialized: total mass = "+
 
453
            cvm::to_str(this->total_mass)+".\n");
428
454
 
429
455
  cvm::decrease_depth();
 
456
 
 
457
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
430
458
}
431
459
 
432
460
 
433
 
void cvm::atom_group::create_sorted_ids (void)
 
461
int cvm::atom_group::create_sorted_ids(void)
434
462
{
435
463
  // Only do the work if the vector is not yet populated
436
464
  if (sorted_ids.size())
437
 
    return;
 
465
    return COLVARS_OK;
438
466
 
439
467
  std::list<int> temp_id_list;
440
 
  for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
441
 
    temp_id_list.push_back (ai->id);
 
468
  cvm::atom_iter ai;
 
469
  for (ai = this->begin(); ai != this->end(); ai++) {
 
470
    temp_id_list.push_back(ai->id);
442
471
  }
443
472
  temp_id_list.sort();
444
473
  temp_id_list.unique();
445
474
  if (temp_id_list.size() != this->size()) {
446
 
    cvm::fatal_error ("Error: duplicate atom IDs in atom group? (found " +
447
 
                      cvm::to_str(temp_id_list.size()) +
448
 
                      " unique atom IDs instead of" +
449
 
                      cvm::to_str(this->size()) + ").\n");
450
 
  }
451
 
  sorted_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
452
 
  return;
 
475
    cvm::error("Error: duplicate atom IDs in atom group? (found " +
 
476
                cvm::to_str(temp_id_list.size()) +
 
477
                " unique atom IDs instead of" +
 
478
                cvm::to_str(this->size()) + ").\n");
 
479
    return COLVARS_ERROR;
 
480
  }
 
481
  sorted_ids = std::vector<int> (temp_id_list.size());
 
482
  unsigned int id_i = 0;
 
483
  std::list<int>::iterator li;
 
484
  for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
 
485
    sorted_ids[id_i] = *li;
 
486
    id_i++;
 
487
  }
 
488
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
453
489
}
454
490
 
455
491
void cvm::atom_group::center_ref_pos()
459
495
  // the rotational fit
460
496
  // This is called either by atom_group::parse or by CVCs that set
461
497
  // reference positions (eg. RMSD, eigenvector)
462
 
  ref_pos_cog = cvm::atom_pos (0.0, 0.0, 0.0);
463
 
  std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
464
 
  for ( ; pi != ref_pos.end(); pi++) {
 
498
  ref_pos_cog = cvm::atom_pos(0.0, 0.0, 0.0);
 
499
  std::vector<cvm::atom_pos>::iterator pi;
 
500
  for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {
465
501
    ref_pos_cog += *pi;
466
502
  }
467
503
  ref_pos_cog /= (cvm::real) ref_pos.size();
468
 
  for (std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
469
 
       pi != ref_pos.end(); pi++) {
 
504
  for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {
470
505
    (*pi) -= ref_pos_cog;
471
506
  }
472
507
}
500
535
  if (b_rotate) {
501
536
    // rotate the group (around the center of geometry if b_center is
502
537
    // true, around the origin otherwise)
503
 
    rot.calc_optimal_rotation (fit_group->positions(), ref_pos);
 
538
    rot.calc_optimal_rotation(fit_group->positions(), ref_pos);
504
539
 
505
540
    for (cvm::atom_iter ai = this->begin();
506
541
         ai != this->end(); ai++) {
507
 
      ai->pos = rot.rotate (ai->pos);
 
542
      ai->pos = rot.rotate(ai->pos);
508
543
    }
509
544
  }
510
545
 
517
552
  }
518
553
}
519
554
 
520
 
void cvm::atom_group::apply_translation (cvm::rvector const &t)
 
555
void cvm::atom_group::apply_translation(cvm::rvector const &t)
521
556
{
522
557
  if (b_dummy) return;
523
558
 
527
562
  }
528
563
}
529
564
 
530
 
void cvm::atom_group::apply_rotation (cvm::rotation const &rot)
 
565
void cvm::atom_group::apply_rotation(cvm::rotation const &rot)
531
566
{
532
567
  if (b_dummy) return;
533
568
 
534
569
  for (cvm::atom_iter ai = this->begin();
535
570
       ai != this->end(); ai++) {
536
 
    ai->pos = rot.rotate (ai->pos);
 
571
    ai->pos = rot.rotate(ai->pos);
537
572
  }
538
573
}
539
574
 
547
582
    for (cvm::atom_iter ai = this->begin();
548
583
         ai != this->end(); ai++) {
549
584
      ai->read_velocity();
550
 
      ai->vel = rot.rotate (ai->vel);
 
585
      ai->vel = rot.rotate(ai->vel);
551
586
    }
552
587
 
553
588
  } else {
568
603
    for (cvm::atom_iter ai = this->begin();
569
604
         ai != this->end(); ai++) {
570
605
      ai->read_system_force();
571
 
      ai->system_force = rot.rotate (ai->system_force);
 
606
      ai->system_force = rot.rotate(ai->system_force);
572
607
    }
573
608
 
574
609
  } else {
585
620
  if (b_dummy)
586
621
    return dummy_atom_pos;
587
622
 
588
 
  cvm::atom_pos cog (0.0, 0.0, 0.0);
 
623
  cvm::atom_pos cog(0.0, 0.0, 0.0);
589
624
  for (cvm::atom_const_iter ai = this->begin();
590
625
       ai != this->end(); ai++) {
591
626
    cog += ai->pos;
599
634
  if (b_dummy)
600
635
    return dummy_atom_pos;
601
636
 
602
 
  cvm::atom_pos com (0.0, 0.0, 0.0);
 
637
  cvm::atom_pos com(0.0, 0.0, 0.0);
603
638
  for (cvm::atom_const_iter ai = this->begin();
604
639
       ai != this->end(); ai++) {
605
640
    com += ai->mass * ai->pos;
609
644
}
610
645
 
611
646
 
612
 
void cvm::atom_group::set_weighted_gradient (cvm::rvector const &grad)
 
647
void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad)
613
648
{
614
649
  if (b_dummy) return;
615
650
 
627
662
  if ((!b_center) && (!b_rotate)) return; // no fit
628
663
 
629
664
  if (cvm::debug())
630
 
    cvm::log ("Calculating fit gradients.\n");
 
665
    cvm::log("Calculating fit gradients.\n");
631
666
 
632
667
  atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
633
 
  group_for_fit->fit_gradients.assign (group_for_fit->size(), cvm::rvector (0.0, 0.0, 0.0));
 
668
  group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::rvector(0.0, 0.0, 0.0));
634
669
 
635
670
  if (b_center) {
636
671
    // add the center of geometry contribution to the gradients
637
672
    for (size_t i = 0; i < this->size(); i++) {
638
673
      // need to bring the gradients in original frame first
639
674
      cvm::rvector const atom_grad = b_rotate ?
640
 
        (rot.inverse()).rotate ((*this)[i].grad) :
 
675
        (rot.inverse()).rotate((*this)[i].grad) :
641
676
        (*this)[i].grad;
642
677
      for (size_t j = 0; j < group_for_fit->size(); j++) {
643
678
        group_for_fit->fit_gradients[j] +=
644
 
          (-1.0)/(cvm::real (group_for_fit->size())) *
 
679
          (-1.0)/(cvm::real(group_for_fit->size())) *
645
680
          atom_grad;
646
681
      }
647
682
    }
655
690
 
656
691
    for (size_t i = 0; i < this->size(); i++) {
657
692
 
658
 
      cvm::atom_pos const pos_orig = rot_inv.rotate ((b_center ? ((*this)[i].pos - cog) : ((*this)[i].pos)));
 
693
      cvm::atom_pos const pos_orig = rot_inv.rotate((b_center ? ((*this)[i].pos - cog) : ((*this)[i].pos)));
659
694
 
660
695
      for (size_t j = 0; j < group_for_fit->size(); j++) {
661
696
        // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
662
697
        cvm::quaternion const dxdq =
663
 
          rot.q.position_derivative_inner (pos_orig, (*this)[i].grad);
 
698
          rot.q.position_derivative_inner(pos_orig, (*this)[i].grad);
664
699
        // multiply by \cdot {\partial q}/\partial\vec{x}_j and add it to the fit gradients
665
700
        for (size_t iq = 0; iq < 4; iq++) {
666
701
          group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq];
669
704
    }
670
705
  }
671
706
  if (cvm::debug())
672
 
    cvm::log ("Done calculating fit gradients.\n");
 
707
    cvm::log("Done calculating fit gradients.\n");
673
708
}
674
709
 
675
710
 
676
711
std::vector<cvm::atom_pos> cvm::atom_group::positions() const
677
712
{
678
 
  if (b_dummy)
679
 
    cvm::fatal_error ("Error: positions are not available "
 
713
  if (b_dummy) {
 
714
    cvm::error("Error: positions are not available "
680
715
                      "from a dummy atom group.\n");
 
716
  }
681
717
 
682
 
  std::vector<cvm::atom_pos> x (this->size(), 0.0);
 
718
  std::vector<cvm::atom_pos> x(this->size(), 0.0);
683
719
  cvm::atom_const_iter ai = this->begin();
684
720
  std::vector<cvm::atom_pos>::iterator xi = x.begin();
685
 
  for ( ; ai != this->end(); xi++, ai++) {
 
721
  for ( ; ai != this->end(); ++xi, ++ai) {
686
722
    *xi = ai->pos;
687
723
  }
688
724
  return x;
689
725
}
690
726
 
691
 
std::vector<cvm::atom_pos> cvm::atom_group::positions_shifted (cvm::rvector const &shift) const
 
727
std::vector<cvm::atom_pos> cvm::atom_group::positions_shifted(cvm::rvector const &shift) const
692
728
{
693
 
  if (b_dummy)
694
 
    cvm::fatal_error ("Error: positions are not available "
695
 
                      "from a dummy atom group.\n");
 
729
  if (b_dummy) {
 
730
    cvm::error("Error: positions are not available "
 
731
                 "from a dummy atom group.\n");
 
732
  }
696
733
 
697
 
  std::vector<cvm::atom_pos> x (this->size(), 0.0);
 
734
  std::vector<cvm::atom_pos> x(this->size(), 0.0);
698
735
  cvm::atom_const_iter ai = this->begin();
699
736
  std::vector<cvm::atom_pos>::iterator xi = x.begin();
700
 
  for ( ; ai != this->end(); xi++, ai++) {
 
737
  for ( ; ai != this->end(); ++xi, ++ai) {
701
738
    *xi = (ai->pos + shift);
702
739
  }
703
740
  return x;
705
742
 
706
743
std::vector<cvm::rvector> cvm::atom_group::velocities() const
707
744
{
708
 
  if (b_dummy)
709
 
    cvm::fatal_error ("Error: velocities are not available "
710
 
                      "from a dummy atom group.\n");
 
745
  if (b_dummy) {
 
746
    cvm::error("Error: velocities are not available "
 
747
                 "from a dummy atom group.\n");
 
748
  }
711
749
 
712
 
  std::vector<cvm::rvector> v (this->size(), 0.0);
 
750
  std::vector<cvm::rvector> v(this->size(), 0.0);
713
751
  cvm::atom_const_iter ai = this->begin();
714
752
  std::vector<cvm::atom_pos>::iterator vi = v.begin();
715
753
  for ( ; ai != this->end(); vi++, ai++) {
720
758
 
721
759
std::vector<cvm::rvector> cvm::atom_group::system_forces() const
722
760
{
723
 
  if (b_dummy)
724
 
    cvm::fatal_error ("Error: system forces are not available "
725
 
                      "from a dummy atom group.\n");
 
761
  if (b_dummy) {
 
762
    cvm::error("Error: system forces are not available "
 
763
                 "from a dummy atom group.\n");
 
764
  }
726
765
 
727
 
  std::vector<cvm::rvector> f (this->size(), 0.0);
 
766
  std::vector<cvm::rvector> f(this->size(), 0.0);
728
767
  cvm::atom_const_iter ai = this->begin();
729
768
  std::vector<cvm::atom_pos>::iterator fi = f.begin();
730
 
  for ( ; ai != this->end(); fi++, ai++) {
 
769
  for ( ; ai != this->end(); ++fi, ++ai) {
731
770
    *fi = ai->system_force;
732
771
  }
733
772
  return f;
735
774
 
736
775
cvm::rvector cvm::atom_group::system_force() const
737
776
{
738
 
  if (b_dummy)
739
 
    cvm::fatal_error ("Error: system forces are not available "
740
 
                      "from a dummy atom group.\n");
 
777
  if (b_dummy) {
 
778
    cvm::error("Error: system forces are not available "
 
779
                "from a dummy atom group.\n");
 
780
  }
741
781
 
742
 
  cvm::rvector f (0.0);
 
782
  cvm::rvector f(0.0);
743
783
  for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
744
784
    f += ai->system_force;
745
785
  }
747
787
}
748
788
 
749
789
 
750
 
void cvm::atom_group::apply_colvar_force (cvm::real const &force)
 
790
void cvm::atom_group::apply_colvar_force(cvm::real const &force)
751
791
{
752
792
  if (b_dummy)
753
793
    return;
754
794
 
755
 
  if (noforce)
756
 
    cvm::fatal_error ("Error: sending a force to a group that has "
757
 
                      "\"enableForces\" set to off.\n");
 
795
  if (noforce) {
 
796
    cvm::error("Error: sending a force to a group that has "
 
797
                "\"enableForces\" set to off.\n");
 
798
    return;
 
799
  }
758
800
 
759
801
  if (b_rotate) {
760
802
 
762
804
    cvm::rotation const rot_inv = rot.inverse();
763
805
    for (cvm::atom_iter ai = this->begin();
764
806
         ai != this->end(); ai++) {
765
 
      ai->apply_force (rot_inv.rotate (force * ai->grad));
 
807
      ai->apply_force(rot_inv.rotate(force * ai->grad));
766
808
    }
767
809
 
768
810
  } else {
769
811
 
770
812
    for (cvm::atom_iter ai = this->begin();
771
813
         ai != this->end(); ai++) {
772
 
      ai->apply_force (force * ai->grad);
 
814
      ai->apply_force(force * ai->grad);
773
815
    }
774
816
  }
775
817
 
782
824
      // rotate forces back to the original frame
783
825
      cvm::rotation const rot_inv = rot.inverse();
784
826
      for (size_t j = 0; j < group_for_fit->size(); j++) {
785
 
        (*group_for_fit)[j].apply_force (rot_inv.rotate (force * group_for_fit->fit_gradients[j]));
 
827
        (*group_for_fit)[j].apply_force(rot_inv.rotate(force * group_for_fit->fit_gradients[j]));
786
828
      }
787
829
    } else {
788
830
      for (size_t j = 0; j < group_for_fit->size(); j++) {
789
 
        (*group_for_fit)[j].apply_force (force * group_for_fit->fit_gradients[j]);
 
831
        (*group_for_fit)[j].apply_force(force * group_for_fit->fit_gradients[j]);
790
832
      }
791
833
    }
792
834
  }
794
836
}
795
837
 
796
838
 
797
 
void cvm::atom_group::apply_force (cvm::rvector const &force)
 
839
void cvm::atom_group::apply_force(cvm::rvector const &force)
798
840
{
799
841
  if (b_dummy)
800
842
    return;
801
843
 
802
 
  if (noforce)
803
 
    cvm::fatal_error ("Error: sending a force to a group that has "
804
 
                      "\"disableForces\" defined.\n");
 
844
  if (noforce) {
 
845
    cvm::error("Error: sending a force to a group that has "
 
846
                "\"disableForces\" defined.\n");
 
847
    return;
 
848
  }
805
849
 
806
850
  if (b_rotate) {
807
851
 
808
852
    cvm::rotation const rot_inv = rot.inverse();
809
853
    for (cvm::atom_iter ai = this->begin();
810
854
         ai != this->end(); ai++) {
811
 
      ai->apply_force (rot_inv.rotate ((ai->mass/this->total_mass) * force));
 
855
      ai->apply_force(rot_inv.rotate((ai->mass/this->total_mass) * force));
812
856
    }
813
857
 
814
858
  } else {
815
859
 
816
860
    for (cvm::atom_iter ai = this->begin();
817
861
         ai != this->end(); ai++) {
818
 
      ai->apply_force ((ai->mass/this->total_mass) * force);
 
862
      ai->apply_force((ai->mass/this->total_mass) * force);
819
863
    }
820
864
  }
821
865
}
822
866
 
823
867
 
824
 
void cvm::atom_group::apply_forces (std::vector<cvm::rvector> const &forces)
 
868
void cvm::atom_group::apply_forces(std::vector<cvm::rvector> const &forces)
825
869
{
826
870
  if (b_dummy)
827
871
    return;
828
872
 
829
873
  if (noforce)
830
 
    cvm::fatal_error ("Error: sending a force to a group that has "
831
 
                      "\"disableForces\" defined.\n");
 
874
    cvm::error("Error: sending a force to a group that has "
 
875
                "\"disableForces\" defined.\n");
832
876
 
833
877
  if (forces.size() != this->size()) {
834
 
    cvm::fatal_error ("Error: trying to apply an array of forces to an atom "
835
 
                      "group which does not have the same length.\n");
 
878
    cvm::error("Error: trying to apply an array of forces to an atom "
 
879
                "group which does not have the same length.\n");
836
880
  }
837
881
 
838
882
  if (b_rotate) {
840
884
    cvm::rotation const rot_inv = rot.inverse();
841
885
    cvm::atom_iter ai = this->begin();
842
886
    std::vector<cvm::rvector>::const_iterator fi = forces.begin();
843
 
    for ( ; ai != this->end(); fi++, ai++) {
844
 
      ai->apply_force (rot_inv.rotate (*fi));
 
887
    for ( ; ai != this->end(); ++fi, ++ai) {
 
888
      ai->apply_force(rot_inv.rotate(*fi));
845
889
    }
846
890
 
847
891
  } else {
848
892
 
849
893
    cvm::atom_iter ai = this->begin();
850
894
    std::vector<cvm::rvector>::const_iterator fi = forces.begin();
851
 
    for ( ; ai != this->end(); fi++, ai++) {
852
 
      ai->apply_force (*fi);
 
895
    for ( ; ai != this->end(); ++fi, ++ai) {
 
896
      ai->apply_force(*fi);
853
897
    }
854
898
  }
855
899
}