← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/colvars/colvarcomp_coordnums.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/colvarcomp_coordnums.cpp
 
1
/// -*- c++ -*-
 
2
 
1
3
#include <cmath>
2
4
 
3
5
#include "colvarmodule.h"
10
12
 
11
13
 
12
14
template<bool calculate_gradients>
13
 
cvm::real colvar::coordnum::switching_function (cvm::real const &r0,
 
15
cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
14
16
                                                int const &en,
15
17
                                                int const &ed,
16
18
                                                cvm::atom &A1,
17
19
                                                cvm::atom &A2)
18
20
{
19
 
  cvm::rvector const diff = cvm::position_distance (A1.pos, A2.pos);
 
21
  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
20
22
  cvm::real const l2 = diff.norm2()/(r0*r0);
21
23
 
22
24
  // Assume en and ed are even integers, and avoid sqrt in the following
23
25
  int const en2 = en/2;
24
26
  int const ed2 = ed/2;
25
27
 
26
 
  cvm::real const xn = std::pow (l2, en2);
27
 
  cvm::real const xd = std::pow (l2, ed2);
 
28
  cvm::real const xn = std::pow(l2, en2);
 
29
  cvm::real const xd = std::pow(l2, ed2);
28
30
  cvm::real const func = (1.0-xn)/(1.0-xd);
29
31
 
30
32
  if (calculate_gradients) {
39
41
 
40
42
 
41
43
template<bool calculate_gradients>
42
 
cvm::real colvar::coordnum::switching_function (cvm::rvector const &r0_vec,
 
44
cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
43
45
                                                int const &en,
44
46
                                                int const &ed,
45
47
                                                cvm::atom &A1,
46
48
                                                cvm::atom &A2)
47
49
{
48
 
  cvm::rvector const diff = cvm::position_distance (A1.pos, A2.pos);
49
 
  cvm::rvector const scal_diff (diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
 
50
  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
 
51
  cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
50
52
  cvm::real const l2 = scal_diff.norm2();
51
53
 
52
54
  // Assume en and ed are even integers, and avoid sqrt in the following
53
55
  int const en2 = en/2;
54
56
  int const ed2 = ed/2;
55
57
 
56
 
  cvm::real const xn = std::pow (l2, en2);
57
 
  cvm::real const xd = std::pow (l2, ed2);
 
58
  cvm::real const xn = std::pow(l2, en2);
 
59
  cvm::real const xd = std::pow(l2, ed2);
58
60
  cvm::real const func = (1.0-xn)/(1.0-xd);
59
61
 
60
62
  if (calculate_gradients) {
61
63
    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
62
 
    cvm::rvector const dl2dx ((2.0/(r0_vec.x*r0_vec.x))*diff.x,
63
 
                              (2.0/(r0_vec.y*r0_vec.y))*diff.y,
64
 
                              (2.0/(r0_vec.z*r0_vec.z))*diff.z);
 
64
    cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
 
65
                             (2.0/(r0_vec.y*r0_vec.y))*diff.y,
 
66
                             (2.0/(r0_vec.z*r0_vec.z))*diff.z);
65
67
    A1.grad += (-1.0)*dFdl2*dl2dx;
66
68
    A2.grad +=        dFdl2*dl2dx;
67
69
  }
70
72
 
71
73
 
72
74
 
73
 
colvar::coordnum::coordnum (std::string const &conf)
74
 
  : distance (conf), b_anisotropic (false), b_group2_center_only (false)
 
75
colvar::coordnum::coordnum(std::string const &conf)
 
76
  : distance(conf), b_anisotropic(false), b_group2_center_only(false)
75
77
{
76
78
  function_type = "coordnum";
77
 
  x.type (colvarvalue::type_scalar);
 
79
  x.type(colvarvalue::type_scalar);
78
80
 
79
81
  // group1 and group2 are already initialized by distance()
80
82
  if (group1.b_dummy)
81
 
    cvm::fatal_error ("Error: only group2 is allowed to be a dummy atom\n");
 
83
    cvm::fatal_error("Error: only group2 is allowed to be a dummy atom\n");
82
84
 
83
85
 
84
86
  // need to specify this explicitly because the distance() constructor
85
87
  // has set it to true
86
88
  b_inverse_gradients = false;
87
89
 
88
 
  bool const b_scale = get_keyval (conf, "cutoff", r0,
89
 
                                   cvm::real (4.0 * cvm::unit_angstrom()));
 
90
  bool const b_scale = get_keyval(conf, "cutoff", r0,
 
91
                                   cvm::real(4.0 * cvm::unit_angstrom()));
90
92
 
91
 
  if (get_keyval (conf, "cutoff3", r0_vec,
92
 
                  cvm::rvector (4.0, 4.0, 4.0), parse_silent)) {
 
93
  if (get_keyval(conf, "cutoff3", r0_vec,
 
94
                  cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {
93
95
 
94
96
    if (b_scale)
95
 
      cvm::fatal_error ("Error: cannot specify \"scale\" and "
 
97
      cvm::fatal_error("Error: cannot specify \"scale\" and "
96
98
                        "\"scale3\" at the same time.\n");
97
99
    b_anisotropic = true;
98
100
    // remove meaningless negative signs
101
103
    if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
102
104
  }
103
105
 
104
 
  get_keyval (conf, "expNumer", en, int (6) );
105
 
  get_keyval (conf, "expDenom", ed, int (12));
 
106
  get_keyval(conf, "expNumer", en, int(6) );
 
107
  get_keyval(conf, "expDenom", ed, int(12));
106
108
 
107
109
  if ( (en%2) || (ed%2) ) {
108
 
    cvm::fatal_error ("Error: odd exponents provided, can only use even ones.\n");
 
110
    cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
109
111
  }
110
112
 
111
 
  get_keyval (conf, "group2CenterOnly", b_group2_center_only, group2.b_dummy);
 
113
  get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2.b_dummy);
112
114
}
113
115
 
114
116
 
115
117
colvar::coordnum::coordnum()
116
 
  : b_anisotropic (false), b_group2_center_only (false)
 
118
  : b_anisotropic(false), b_group2_center_only(false)
117
119
{
118
120
  function_type = "coordnum";
119
 
  x.type (colvarvalue::type_scalar);
 
121
  x.type(colvarvalue::type_scalar);
120
122
}
121
123
 
122
124
 
132
134
 
133
135
    if (b_anisotropic) {
134
136
      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
135
 
        x.real_value += switching_function<false> (r0_vec, en, ed, *ai1, group2_com_atom);
 
137
        x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, group2_com_atom);
136
138
    } else {
137
139
      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
138
 
        x.real_value += switching_function<false> (r0, en, ed, *ai1, group2_com_atom);
 
140
        x.real_value += switching_function<false>(r0, en, ed, *ai1, group2_com_atom);
139
141
    }
140
142
 
141
143
  } else {
143
145
    if (b_anisotropic) {
144
146
      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
145
147
        for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
146
 
          x.real_value += switching_function<false> (r0_vec, en, ed, *ai1, *ai2);
 
148
          x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, *ai2);
147
149
        }
148
150
    } else {
149
151
      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
150
152
        for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
151
 
          x.real_value += switching_function<false> (r0, en, ed, *ai1, *ai2);
 
153
          x.real_value += switching_function<false>(r0, en, ed, *ai1, *ai2);
152
154
        }
153
155
    }
154
156
  }
166
168
 
167
169
    if (b_anisotropic) {
168
170
      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
169
 
        switching_function<true> (r0_vec, en, ed, *ai1, group2_com_atom);
 
171
        switching_function<true>(r0_vec, en, ed, *ai1, group2_com_atom);
170
172
    } else {
171
173
      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
172
 
        switching_function<true> (r0, en, ed, *ai1, group2_com_atom);
 
174
        switching_function<true>(r0, en, ed, *ai1, group2_com_atom);
173
175
    }
174
176
 
175
 
    group2.set_weighted_gradient (group2_com_atom.grad);
 
177
    group2.set_weighted_gradient(group2_com_atom.grad);
176
178
 
177
179
  } else {
178
180
 
179
181
    if (b_anisotropic) {
180
182
      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
181
183
        for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
182
 
          switching_function<true> (r0_vec, en, ed, *ai1, *ai2);
 
184
          switching_function<true>(r0_vec, en, ed, *ai1, *ai2);
183
185
        }
184
186
    } else {
185
187
      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
186
188
        for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
187
 
          switching_function<true> (r0, en, ed, *ai1, *ai2);
 
189
          switching_function<true>(r0, en, ed, *ai1, *ai2);
188
190
        }
189
191
    }
190
192
  }
191
 
 
192
 
  //   if (cvm::debug()) {
193
 
  //     for (size_t i = 0; i < group1.size(); i++) {
194
 
  //       cvm::log ("atom["+cvm::to_str (group1[i].id+1)+"] gradient: "+
195
 
  //                 cvm::to_str (group1[i].grad)+"\n");
196
 
  //     }
197
 
 
198
 
  //     for (size_t i = 0; i < group2.size(); i++) {
199
 
  //       cvm::log ("atom["+cvm::to_str (group2[i].id+1)+"] gradient: "+
200
 
  //                 cvm::to_str (group2[i].grad)+"\n");
201
 
  //     }
202
 
  //   }
203
193
}
204
194
 
205
 
void colvar::coordnum::apply_force (colvarvalue const &force)
 
195
void colvar::coordnum::apply_force(colvarvalue const &force)
206
196
{
207
197
  if (!group1.noforce)
208
 
    group1.apply_colvar_force (force.real_value);
 
198
    group1.apply_colvar_force(force.real_value);
209
199
 
210
200
  if (!group2.noforce)
211
 
    group2.apply_colvar_force (force.real_value);
 
201
    group2.apply_colvar_force(force.real_value);
212
202
}
213
203
 
214
204
 
215
205
 
216
206
// h_bond member functions
217
207
 
218
 
colvar::h_bond::h_bond (std::string const &conf)
219
 
  : cvc (conf)
 
208
colvar::h_bond::h_bond(std::string const &conf)
 
209
  : cvc(conf)
220
210
{
221
211
  if (cvm::debug())
222
 
    cvm::log ("Initializing h_bond object.\n");
 
212
    cvm::log("Initializing h_bond object.\n");
223
213
 
224
214
  function_type = "h_bond";
225
 
  x.type (colvarvalue::type_scalar);
 
215
  x.type(colvarvalue::type_scalar);
226
216
 
227
217
  int a_num, d_num;
228
 
  get_keyval (conf, "acceptor", a_num, -1);
229
 
  get_keyval (conf, "donor",    d_num, -1);
 
218
  get_keyval(conf, "acceptor", a_num, -1);
 
219
  get_keyval(conf, "donor",    d_num, -1);
230
220
 
231
221
  if ( (a_num == -1) || (d_num == -1) ) {
232
 
    cvm::fatal_error ("Error: either acceptor or donor undefined.\n");
 
222
    cvm::fatal_error("Error: either acceptor or donor undefined.\n");
233
223
  }
234
224
 
235
 
  acceptor = cvm::atom (a_num);
236
 
  donor    = cvm::atom (d_num);
237
 
  atom_groups.push_back (new cvm::atom_group);
238
 
  atom_groups[0]->add_atom (acceptor);
239
 
  atom_groups[0]->add_atom (donor);
 
225
  cvm::atom acceptor = cvm::atom(a_num);
 
226
  cvm::atom donor    = cvm::atom(d_num);
 
227
  atom_groups.push_back(new cvm::atom_group);
 
228
  atom_groups[0]->add_atom(acceptor);
 
229
  atom_groups[0]->add_atom(donor);
240
230
 
241
 
  get_keyval (conf, "cutoff",   r0, (3.3 * cvm::unit_angstrom()));
242
 
  get_keyval (conf, "expNumer", en, 6);
243
 
  get_keyval (conf, "expDenom", ed, 8);
 
231
  get_keyval(conf, "cutoff",   r0, (3.3 * cvm::unit_angstrom()));
 
232
  get_keyval(conf, "expNumer", en, 6);
 
233
  get_keyval(conf, "expDenom", ed, 8);
244
234
 
245
235
  if ( (en%2) || (ed%2) ) {
246
 
    cvm::fatal_error ("Error: odd exponents provided, can only use even ones.\n");
 
236
    cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
247
237
  }
248
238
 
249
239
  if (cvm::debug())
250
 
    cvm::log ("Done initializing h_bond object.\n");
 
240
    cvm::log("Done initializing h_bond object.\n");
251
241
}
252
242
 
253
243
 
254
 
colvar::h_bond::h_bond (cvm::atom const &acceptor_i,
255
 
                        cvm::atom const &donor_i,
256
 
                        cvm::real r0_i, int en_i, int ed_i)
257
 
  : acceptor (acceptor_i),
258
 
    donor (donor_i),
259
 
    r0 (r0_i), en (en_i), ed (ed_i)
 
244
colvar::h_bond::h_bond(cvm::atom const &acceptor,
 
245
                       cvm::atom const &donor,
 
246
                       cvm::real r0_i, int en_i, int ed_i)
 
247
  : r0(r0_i), en(en_i), ed(ed_i)
260
248
{
261
249
  function_type = "h_bond";
262
 
  x.type (colvarvalue::type_scalar);
 
250
  x.type(colvarvalue::type_scalar);
263
251
 
264
 
  atom_groups.push_back (new cvm::atom_group);
265
 
  atom_groups[0]->add_atom (acceptor);
266
 
  atom_groups[0]->add_atom (donor);
 
252
  atom_groups.push_back(new cvm::atom_group);
 
253
  atom_groups[0]->add_atom(acceptor);
 
254
  atom_groups[0]->add_atom(donor);
267
255
}
268
256
 
269
257
colvar::h_bond::h_bond()
270
 
  : cvc ()
 
258
  : cvc()
271
259
{
272
260
  function_type = "h_bond";
273
 
  x.type (colvarvalue::type_scalar);
 
261
  x.type(colvarvalue::type_scalar);
274
262
}
275
263
 
276
264
 
277
265
colvar::h_bond::~h_bond()
278
266
{
279
 
  for (int i=0; i<atom_groups.size(); i++) {
280
 
    delete atom_groups[i];
281
 
  }
 
267
  delete atom_groups[0];
282
268
}
283
269
 
284
270
 
285
271
void colvar::h_bond::calc_value()
286
272
{
287
 
  x.real_value = colvar::coordnum::switching_function<false> (r0, en, ed, acceptor, donor);
 
273
  x.real_value = colvar::coordnum::switching_function<false>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
288
274
}
289
275
 
290
276
 
291
277
void colvar::h_bond::calc_gradients()
292
278
{
293
 
  colvar::coordnum::switching_function<true> (r0, en, ed, acceptor, donor);
294
 
  (*atom_groups[0])[0].grad = acceptor.grad;
295
 
  (*atom_groups[0])[1].grad = donor.grad;
 
279
  colvar::coordnum::switching_function<true>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
296
280
}
297
281
 
298
282
 
299
 
void colvar::h_bond::apply_force (colvarvalue const &force)
 
283
void colvar::h_bond::apply_force(colvarvalue const &force)
300
284
{
301
 
  acceptor.apply_force (force.real_value * acceptor.grad);
302
 
  donor.apply_force    (force.real_value * donor.grad);
 
285
  (atom_groups[0])->apply_colvar_force(force);
303
286
}
304
287
 
305
288
 
306
289
 
307
290
 
308
 
colvar::selfcoordnum::selfcoordnum (std::string const &conf)
309
 
 : distance (conf, false)
 
291
colvar::selfcoordnum::selfcoordnum(std::string const &conf)
 
292
 : distance(conf, false)
310
293
{
311
294
  function_type = "selfcoordnum";
312
 
  x.type (colvarvalue::type_scalar);
 
295
  x.type(colvarvalue::type_scalar);
313
296
 
314
297
  // group1 is already initialized by distance()
315
298
 
317
300
  // has set it to true
318
301
  b_inverse_gradients = false;
319
302
 
320
 
  get_keyval (conf, "cutoff", r0, cvm::real (4.0 * cvm::unit_angstrom()));
321
 
  get_keyval (conf, "expNumer", en, int (6) );
322
 
  get_keyval (conf, "expDenom", ed, int (12));
 
303
  get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom()));
 
304
  get_keyval(conf, "expNumer", en, int(6) );
 
305
  get_keyval(conf, "expDenom", ed, int(12));
323
306
 
324
307
  if ( (en%2) || (ed%2) ) {
325
 
    cvm::fatal_error ("Error: odd exponents provided, can only use even ones.\n");
 
308
    cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
326
309
  }
327
310
}
328
311
 
330
313
colvar::selfcoordnum::selfcoordnum()
331
314
{
332
315
  function_type = "selfcoordnum";
333
 
  x.type (colvarvalue::type_scalar);
 
316
  x.type(colvarvalue::type_scalar);
334
317
}
335
318
 
336
319
 
337
320
void colvar::selfcoordnum::calc_value()
338
321
{
339
322
  x.real_value = 0.0;
340
 
 
341
 
  for (size_t i = 0; i < group1.size() - 1; i++)
342
 
    for (size_t j = i + 1; j < group1.size(); j++)
343
 
      x.real_value += colvar::coordnum::switching_function<false> (r0, en, ed, group1[i], group1[j]);
 
323
  for (size_t i = 0; i < group1.size() - 1; i++) {
 
324
    for (size_t j = i + 1; j < group1.size(); j++) {
 
325
      x.real_value += colvar::coordnum::switching_function<false>(r0, en, ed, group1[i], group1[j]);
 
326
    }
 
327
  }
344
328
}
345
329
 
346
330
 
347
331
void colvar::selfcoordnum::calc_gradients()
348
332
{
349
 
  for (size_t i = 0; i < group1.size() - 1; i++)
350
 
    for (size_t j = i + 1; j < group1.size(); j++)
351
 
      colvar::coordnum::switching_function<true> (r0, en, ed, group1[i], group1[j]);
 
333
  for (size_t i = 0; i < group1.size() - 1; i++) {
 
334
    for (size_t j = i + 1; j < group1.size(); j++) {
 
335
      colvar::coordnum::switching_function<true>(r0, en, ed, group1[i], group1[j]);
 
336
    }
 
337
  }
352
338
}
353
339
 
354
 
void colvar::selfcoordnum::apply_force (colvarvalue const &force)
 
340
void colvar::selfcoordnum::apply_force(colvarvalue const &force)
355
341
{
356
 
  if (!group1.noforce)
357
 
    group1.apply_colvar_force (force.real_value);
 
342
  if (!group1.noforce) {
 
343
    group1.apply_colvar_force(force.real_value);
 
344
  }
358
345
}
359
346