ConstraintSet.C 20.2 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17 18 19
// ConstraintSet.C
//
////////////////////////////////////////////////////////////////////////////////

#include "ConstraintSet.h"
20
#include "PositionConstraint.h"
21 22 23 24 25 26 27 28
#include "DistanceConstraint.h"
#include "AngleConstraint.h"
#include "TorsionConstraint.h"
#include "Atom.h"
#include "AtomSet.h"
#include "Context.h"
#include <iostream>
#include <iomanip>
29
#include <cstdlib> // atof
30
#include <cstring> // strcmp
31 32
using namespace std;

33
const int constraints_maxiter = 50;
34

35 36 37 38 39 40 41
////////////////////////////////////////////////////////////////////////////////
ConstraintSet::~ConstraintSet(void)
{
  for ( int ic = 0; ic < constraint_list.size(); ic++ )
    delete constraint_list[ic];
}

42
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
43
bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv)
44
{
45
  enum constraint_type { unknown, position_type, distance_type,
46 47
                         angle_type, torsion_type }
    type = unknown;
48
  const double position_tolerance = 1.0e-7;
49 50 51
  const double distance_tolerance = 1.0e-7;
  const double angle_tolerance = 1.0e-4;
  const bool onpe0 = ctxt_.onpe0();
52

Francois Gygi committed
53 54
  // argv[0] == "constraint"
  // argv[1] == "define"
55
  // argv[2] == {"position", "distance", "angle", "torsion"}
Francois Gygi committed
56 57 58 59
  // argv[3] == constraint name
  // argv[4-(5,6,7)] == atom names
  // argv[{5,6,7}] == {distance,angle,angle}
  // argv[{6,7,8}] == velocity
60 61 62 63 64

  if ( argc < 2 )
  {
    if ( onpe0 )
    {
65
      cout << " Use: constraint define position constraint_name atom_name"
66
           << endl;
67 68
      cout << " Use: constraint define distance constraint_name "
           << "atom_name1 atom_name2 distance_value [velocity]"
69
           << endl;
70 71
      cout << "      constraint define angle constraint_name "
           << "name1 name2 name3 angle_value [velocity]"
72
           << endl;
73 74
      cout << "      constraint define torsion constraint_name "
           << "name1 name2 name3 name4 angle_value"
75
           << " [velocity] "
76 77
           << endl;
    }
Francois Gygi committed
78
    return false;
79
  }
Francois Gygi committed
80
  const string constraint_type = argv[2];
81 82 83 84 85
  if ( constraint_type == "position" )
  {
    type = position_type;
  }
  else if ( constraint_type == "distance" )
86 87 88
  {
    type = distance_type;
  }
Francois Gygi committed
89
  else if ( constraint_type == "angle" )
90 91 92
  {
    type = angle_type;
  }
Francois Gygi committed
93
  else if ( constraint_type == "torsion" )
94 95 96 97 98
  {
    type = torsion_type;
  }
  else
  {
99
    if ( onpe0 )
100
      cout << " Incorrect constraint type " << constraint_type << endl;
101 102
    return false;
  }
103

104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
  if ( type == position_type )
  {
    // define position name A

    if ( argc != 5 )
    {
      if ( onpe0 )
        cout << " Incorrect number of arguments for position constraint"
             << endl;
      return false;
    }
    string name = argv[3];
    string name1 = argv[4];

    Atom *a1 = atoms.findAtom(name1);

    if ( a1 == 0 )
    {
      if ( onpe0 )
      {
        cout << " ConstraintSet: could not find atom " << name1 << endl;
        cout << " ConstraintSet: could not define constraint" << endl;
      }
      return false;
    }

    // check if constraint is already defined
    bool found = false;
    Constraint *pc = 0;
    for ( int i = 0; i < constraint_list.size(); i++ )
    {
      pc = constraint_list[i];
      assert(pc != 0);
      // check if a constraint with same name or with same atom is defined
      if ( pc->type() == "position" )
        found = ( pc->name() == name ) || ( pc->names(0) == name1 );
    }

    if ( found )
    {
      if ( onpe0 )
        cout << " ConstraintSet: constraint is already defined:\n"
             << " cannot define constraint" << endl;
      return false;
    }
    else
    {
      PositionConstraint *c =
        new PositionConstraint(name,name1,position_tolerance);
      constraint_list.push_back(c);
    }
  }
  else if ( type == distance_type )
157
  {
158 159
    // define distance name A B value
    // define distance name A B value velocity
160

Francois Gygi committed
161
    if ( argc < 7 || argc > 8 )
162
    {
163
      if ( onpe0 )
164 165
        cout << " Incorrect number of arguments for distance constraint"
             << endl;
166 167
      return false;
    }
168
    double distance, velocity=0.0;
Francois Gygi committed
169 170 171
    string name = argv[3];
    string name1 = argv[4];
    string name2 = argv[5];
172

173 174
    Atom *a1 = atoms.findAtom(name1);
    Atom *a2 = atoms.findAtom(name2);
175

176 177
    if ( a1 == 0 || a2 == 0 )
    {
Francois Gygi committed
178 179 180
      if ( onpe0 )
      {
        if ( a1 == 0 )
181
          cout << " ConstraintSet: could not find atom " << name1 << endl;
Francois Gygi committed
182
        if ( a2 == 0 )
183 184
          cout << " ConstraintSet: could not find atom " << name2 << endl;
        cout << " ConstraintSet: could not define constraint" << endl;
Francois Gygi committed
185 186
      }
      return false;
187 188 189
    }
    if ( name1 == name2 )
    {
190
      if ( onpe0 )
191 192
        cout << " ConstraintSet: cannot define distance constraint between "
             << name1 << " and " << name2 << endl;
Francois Gygi committed
193
      return false;
194
    }
195

196 197 198 199 200 201 202 203 204 205 206 207
    // define distance value
    if ( !strcmp(argv[6],"*") )
    {
      // use current distance
      distance = length(a1->position()-a2->position());
      if ( onpe0 )
        cout << "ConstraintSet::define_constraint: using current distance "
             << distance << endl;
    }
    else
      distance = atof(argv[6]);

Francois Gygi committed
208
    if ( argc == 8 )
209
    {
Francois Gygi committed
210
      velocity = atof(argv[7]);
211
    }
212

213 214 215
    if ( distance <= 0.0 )
    {
      if ( onpe0 )
216 217
        cout << " ConstraintSet: distance must be positive" << endl
             << " ConstraintSet: could not define constraint" << endl;
Francois Gygi committed
218
      return false;
219
    }
220

Francois Gygi committed
221
    // check if constraint is already defined
222 223 224 225 226 227 228
    bool found = false;
    Constraint *pc = 0;
    for ( int i = 0; i < constraint_list.size(); i++ )
    {
      pc = constraint_list[i];
      assert(pc != 0);
      // check if a constraint (name1,name2) or (name2,name1) is defined
229 230 231 232
      if ( pc->type() == "distance" )
        found = ( pc->names(0) == name1 && pc->names(1) == name2 ) ||
                ( pc->names(1) == name1 && pc->names(0) == name2 ) ||
                ( pc->name() == name );
233
    }
234

235 236 237
    if ( found )
    {
      if ( onpe0 )
238 239 240 241
        cout << " ConstraintSet: a distance constraint named " << name << endl
             << " or involving atoms " << name1 << " " << name2 << endl
             << " is already defined" << endl
             << " ConstraintSet: cannot define constraint" << endl;
Francois Gygi committed
242
      return false;
243 244 245 246
    }
    else
    {
      DistanceConstraint *c =
Francois Gygi committed
247
        new DistanceConstraint(name,name1,name2,distance,
248
                               velocity,distance_tolerance);
249

250 251 252 253 254
      constraint_list.push_back(c);
    }
  }
  else if ( type == angle_type )
  {
255 256
    // constraint define angle name A B C value
    // constraint define angle name A B C value velocity
257

Francois Gygi committed
258
    if ( argc < 8  || argc > 9 )
259 260
    {
      if ( onpe0 )
261
        cout << " Incorrect number of arguments for angle constraint"
262 263 264
             << endl;
      return false;
    }
Francois Gygi committed
265 266 267 268
    string name = argv[3];
    string name1 = argv[4];
    string name2 = argv[5];
    string name3 = argv[6];
269

270 271 272
    Atom *a1 = atoms.findAtom(name1);
    Atom *a2 = atoms.findAtom(name2);
    Atom *a3 = atoms.findAtom(name3);
273

274 275
    if ( a1 == 0 || a2 == 0 || a3 == 0 )
    {
276
      if ( onpe0 )
Francois Gygi committed
277 278
      {
        if ( a1 == 0 )
279
          cout << " ConstraintSet: could not find atom " << name1 << endl;
Francois Gygi committed
280
        if ( a2 == 0 )
281
          cout << " ConstraintSet: could not find atom " << name2 << endl;
Francois Gygi committed
282
        if ( a3 == 0 )
283 284
          cout << " ConstraintSet: could not find atom " << name3 << endl;
        cout << " ConstraintSet: could not define constraint" << endl;
Francois Gygi committed
285 286
      }
      return false;
287
    }
288

289 290
    if ( name1 == name2 || name1 == name3 || name2 == name3)
    {
291
      if ( onpe0 )
292 293
        cout << " ConstraintSet: cannot define angle constraint between "
             << name1 << " " << name2 << " and " << name3 << endl;
Francois Gygi committed
294
      return false;
295
    }
296

297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323
    // define angle value
    double angle = 0.0;
    if ( !strcmp(argv[7],"*") )
    {
      // use current angle
      D3vector r12(a1->position()-a2->position());
      D3vector r32(a3->position()-a2->position());
      if ( norm2(r12) == 0.0 || norm2(r32) == 0.0 )
      {
        if ( onpe0 )
        {
          cout << " ConstraintSet: cannot define angle constraint." << endl;
          cout << " ConstraintSet: atoms are too close" << endl;
          return false;
        }
      }
      const double sp = normalized(r12) * normalized(r32);
      const double c = max(-1.0,min(1.0,sp));
      angle = (180.0/M_PI)*acos(c);

      if ( onpe0 )
        cout << "ConstraintSet::define_constraint: using current angle "
             << angle << endl;
    }
    else
      angle = atof(argv[7]);

Francois Gygi committed
324 325
    double velocity = 0.0;
    if ( argc == 9 )
326
    {
Francois Gygi committed
327
      velocity = atof(argv[8]);
328
    }
329

330 331 332
    if ( angle < 0.0 || angle > 180.0 )
    {
      if ( onpe0 )
333 334
        cout << " ConstraintSet: angle must be in [0,180]" << endl
             << " ConstraintSet: could not define constraint" << endl;
Francois Gygi committed
335
      return false;
336
    }
337

338 339 340 341 342 343 344
    // check if equivalent constraint is already defined
    bool found = false;
    Constraint *pc = 0;
    for ( int i = 0; i < constraint_list.size(); i++ )
    {
      pc = constraint_list[i];
      assert(pc != 0);
345
      // check if a constraint (name1,name2,name3) or
346
      // (name3,name2,name1) is defined
347 348 349 350 351 352 353 354
      if ( pc->type() == "angle" )
        found = ( pc->names(0) == name1 &&
                  pc->names(1) == name2 &&
                  pc->names(2) == name3 ) ||
                ( pc->names(0) == name3 &&
                  pc->names(1) == name2 &&
                  pc->names(2) == name1 ) ||
                ( pc->name() == name );
355
    }
356

357 358 359
    if ( found )
    {
      if ( onpe0 )
360 361 362 363
        cout << " ConstraintSet: an angle constraint named " << name << endl
             << " or involving atoms "
             << name1 << " " << name2 << " " << name3 << endl
             << " is already defined" << endl
364
             << " ConstraintSet: cannot define constraint" << endl;
Francois Gygi committed
365
      return false;
366 367 368 369
    }
    else
    {
      AngleConstraint *c =
Francois Gygi committed
370 371
      new AngleConstraint(name, name1,name2,name3,angle,
        velocity,angle_tolerance);
372 373 374 375 376
      constraint_list.push_back(c);
    }
  }
  else if ( type == torsion_type )
  {
377 378
    // constraint define torsion name A B C D angle
    // constraint define torsion name A B C D angle velocity
379

Francois Gygi committed
380
    if ( argc < 9  || argc > 10 )
381 382
    {
      if ( onpe0 )
383
        cout << " Incorrect number of arguments for torsion constraint"
384 385 386
             << endl;
      return false;
    }
Francois Gygi committed
387 388 389 390 391
    string name = argv[3];
    string name1 = argv[4];
    string name2 = argv[5];
    string name3 = argv[6];
    string name4 = argv[7];
392

393 394 395 396
    Atom *a1 = atoms.findAtom(name1);
    Atom *a2 = atoms.findAtom(name2);
    Atom *a3 = atoms.findAtom(name3);
    Atom *a4 = atoms.findAtom(name4);
397

398 399
    if ( a1 == 0 || a2 == 0 || a3 == 0 || a4 == 0 )
    {
Francois Gygi committed
400
      if ( onpe0 )
401
      {
Francois Gygi committed
402
        if ( a1 == 0 )
403
          cout << " ConstraintSet: could not find atom " << name1 << endl;
Francois Gygi committed
404
        if ( a2 == 0 )
405
          cout << " ConstraintSet: could not find atom " << name2 << endl;
Francois Gygi committed
406
        if ( a3 == 0 )
407
          cout << " ConstraintSet: could not find atom " << name3 << endl;
Francois Gygi committed
408
        if ( a4 == 0 )
409 410
          cout << " ConstraintSet: could not find atom " << name4 << endl;
        cout << " ConstraintSet: could not define constraint" << endl;
Francois Gygi committed
411 412
      }
      return false;
413 414 415 416
    }
    if ( name1 == name2 || name1 == name3 || name1 == name4 ||
         name2 == name3 || name2 == name4 || name3 == name4 )
    {
417
      if ( onpe0 )
418
        cout << " ConstraintSet: cannot define torsion constraint using "
419
             << name1 << " " << name2 << " " << name3 << " " << name4
420
             << endl;
Francois Gygi committed
421
      return false;
422
    }
423

424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465
    // define angle value
    double angle = 0.0;
    if ( !strcmp(argv[8],"*") )
    {
      // use current angle
      D3vector r12(a1->position()-a2->position());
      D3vector r32(a3->position()-a2->position());
      D3vector r43(a4->position()-a3->position());
      if ( norm2(r12) == 0.0 || norm2(r32) == 0.0 || norm2(r43) == 0.0 )
      {
        if ( onpe0 )
        {
          cout << " ConstraintSet: cannot define torsion constraint." << endl;
          cout << " ConstraintSet: atoms are too close" << endl;
          return false;
        }
      }

      D3vector e12(normalized(r12));
      D3vector e32(normalized(r32));
      D3vector e23(-e32);
      D3vector e43(normalized(r43));

      const double sin123 = length(e12^e32);
      const double sin234 = length(e23^e43);

      if ( sin123 != 0.0 && sin234 != 0.0 )
      {
        D3vector e123 = normalized(e12^e32);
        D3vector e234 = normalized(e23^e43);
        double cc = max(min(e123*e234,1.0),-1.0);
        double ss = max(min((e123^e234)*e32,1.0),-1.0);
        angle = (180.0/M_PI) * atan2(ss,cc);
      }

      if ( onpe0 )
        cout << "ConstraintSet::define_constraint: using current angle "
             << angle << endl;
    }
    else
      angle = atof(argv[8]);

466 467 468 469 470
    if ( angle > 180.0 )
      while ( angle > 180.0 ) angle -= 360.0;
    else if ( angle < -180.0 )
      while ( angle < -180.0 ) angle += 360.0;

Francois Gygi committed
471 472
    double velocity = 0.0;
    if ( argc == 10 )
473
    {
Francois Gygi committed
474
      velocity = atof(argv[9]);
475
    }
476

477 478 479 480 481 482 483
    // check if equivalent constraint is already defined
    bool found = false;
    Constraint *pc = 0;
    for ( int i = 0; i < constraint_list.size(); i++ )
    {
      pc = constraint_list[i];
      assert(pc != 0);
484
      // check if an equivalent constraint (name1,name2,name3,name4) or
485
      // (name4,name3,name2,name1) is defined
486
      if ( pc->type() == "torsion" )
487 488 489 490 491 492 493 494 495
        found = ( pc->names(0) == name1 &&
                  pc->names(1) == name2 &&
                  pc->names(2) == name3 &&
                  pc->names(3) == name4 ) ||
                ( pc->names(0) == name4 &&
                  pc->names(1) == name3 &&
                  pc->names(2) == name2 &&
                  pc->names(3) == name1 ) ||
                ( pc->name() == name );
496
    }
497

498 499 500
    if ( found )
    {
      if ( onpe0 )
501 502 503
        cout << " ConstraintSet: a torsion constraint named " << name << endl
             << " or involving atoms "
             << name1 << " " << name2 << " " << name3 << " " << name4 << endl
504
             << " is already defined" << endl
505
             << " ConstraintSet: cannot define constraint" << endl;
Francois Gygi committed
506
      return false;
507 508 509 510
    }
    else
    {
      TorsionConstraint *c =
Francois Gygi committed
511
      new TorsionConstraint(name,name1,name2,name3,name4,
512 513 514 515 516 517 518
                            angle,velocity,angle_tolerance);
      constraint_list.push_back(c);
    }
  }
  else
  {
    if ( onpe0 )
519
      cout << " ConstraintSet::set_constraint: internal error" << endl;
520 521
    return false;
  }
522

523 524 525
  // update total number of blocked degrees of freedom
  ndofs_ += constraint_list.back()->dofs();

526 527 528 529
  return true;
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
530
bool ConstraintSet::set_constraint(int argc, char **argv)
531 532
{
  const bool onpe0 = ctxt_.onpe0();
Francois Gygi committed
533 534 535 536 537 538 539 540 541 542 543
  assert(argc==4||argc==5);
  // argv[0] == "constraint"
  // argv[1] == "set"
  // argv[2] == constraint_name
  // argv[3] == value
  // argv[4] (optional) == velocity
  string name = argv[2];
  const double value = atof(argv[3]);
  double velocity = 0.0;
  const bool set_velocity = ( argc == 5 );
  if ( set_velocity ) velocity = atof(argv[4]);
544

Francois Gygi committed
545 546
    // check if constraint is already defined
  bool found = false;
547 548 549 550 551 552 553
  vector<Constraint*>::iterator i = constraint_list.begin();
  while ( !found && i != constraint_list.end() )
  {
    Constraint *pc = *i;
    assert(pc != 0);

    if ( pc->name() == name )
Francois Gygi committed
554
    {
555
      found = true;
Francois Gygi committed
556 557 558 559
      pc->set_value(value);
      if ( set_velocity ) pc->set_velocity(velocity);
    }
    i++;
560
  }
561 562 563 564

  if ( !found )
  {
    if ( onpe0 )
565
      cout << " ConstraintSet: no such constraint" << endl;
566 567
    return false;
  }
Francois Gygi committed
568 569 570 571 572 573 574 575 576 577 578 579 580
  return true;
}

////////////////////////////////////////////////////////////////////////////////
bool ConstraintSet::delete_constraint(int argc, char **argv)
{
  assert(argc==3);
  // argv[0] == "constraint"
  // argv[1] == "delete"
  // argv[2] == constraint_name
  string name = argv[2];
  const bool onpe0 = ctxt_.onpe0();

581 582 583 584 585 586 587 588 589 590 591 592 593
  bool found = false;
  // note next loop in reverse: avoid use of invalidated iterators
  // after erase operation

  vector<Constraint*>::iterator i = constraint_list.begin();
  while ( !found && i != constraint_list.end() )
  {
    Constraint *pc = *i;
    assert(pc != 0);

    // note structure of if else test to avoid incrementing
    // invalidated iterator after erase (see Meyers STL, p.45)
    if ( pc->name() == name )
Francois Gygi committed
594
    {
595
      found = true;
596 597 598 599

      // update total number of blocked degrees of freedom
      ndofs_ -= pc->dofs();

600 601 602 603 604 605 606 607 608
      delete pc;

      // remove constraint pointer from the list
      // note: iterator is incremented before erasing, remains valid
      constraint_list.erase(i++);
    }
    else
    {
      i++;
Francois Gygi committed
609
    }
610
  }
611

Francois Gygi committed
612
  if ( !found )
613
  {
614
    if ( onpe0 ) cout << " No such constraint" << endl;
615 616 617 618 619 620 621 622
    return false;
  }
  return true;
}

////////////////////////////////////////////////////////////////////////////////
void ConstraintSet::list_constraints(ostream &os)
{
Francois Gygi committed
623
  if ( !constraint_list.empty() )
624
  {
Francois Gygi committed
625 626 627 628 629 630 631
    os << " <constraint_set>" << endl;
    for ( int i = 0; i < constraint_list.size(); i++ )
    {
      Constraint *c = constraint_list[i];
      os << *c << endl;
    }
    os << " </constraint_set>" << endl;
632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654
  }
}

////////////////////////////////////////////////////////////////////////////////
void ConstraintSet::enforce(AtomSet& atoms)
{
  vector<vector<double> > r0,rp,v0;
  setup(atoms);
  atoms.get_positions(r0);
  rp=r0;
  atoms.get_velocities(v0);
  enforce_r(r0,rp);
  atoms.set_positions(rp);
  enforce_v(r0,v0);
  atoms.set_velocities(v0);
}
////////////////////////////////////////////////////////////////////////////////
void ConstraintSet::enforce_r(const vector<vector<double> > &r0,
                              vector<vector<double> > &rp)
{
  const bool onpe0 = ctxt_.onpe0();
  int iter = 0;
  bool done = false;
655
  while ( !done && (iter < constraints_maxiter) )
656 657 658 659 660 661 662 663 664 665
  {
    done = true;
    for ( int i = 0; i < constraint_list.size(); i++ )
    {
      Constraint *c = constraint_list[i];
      bool b = c->enforce_r(r0,rp);
      done &= b;
    }
    iter++;
  }
666

667 668 669
  if ( !done )
  {
    if ( onpe0 )
670 671
      cout << " ConstraintSet: could not enforce position constraints in "
           << constraints_maxiter << " iterations" << endl;
672 673 674 675 676 677 678 679 680 681
  }
}

////////////////////////////////////////////////////////////////////////////////
void ConstraintSet::enforce_v(const vector<vector<double> > &r0,
                              vector<vector<double> > &v0)
{
  const bool onpe0 = ctxt_.onpe0();
  int iter = 0;
  bool done = false;
682
  while ( !done && (iter < constraints_maxiter) )
683 684 685 686 687 688 689 690 691
  {
    done = true;
    for ( int i = 0; i < constraint_list.size(); i++ )
    {
      bool b = constraint_list[i]->enforce_v(r0,v0);
      done &= b;
    }
    iter++;
  }
692

693 694 695
  if ( !done )
  {
    if ( onpe0 )
696 697
      cout << " ConstraintSet: could not enforce velocity constraints in "
           << constraints_maxiter << " iterations" << endl;
698 699 700 701
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
702 703
void ConstraintSet::compute_forces(const vector<vector<double> > &r0,
 const vector<vector<double> > &f)
704 705 706
{
  for ( int i = 0; i < constraint_list.size(); i++ )
  {
Francois Gygi committed
707
    constraint_list[i]->compute_force(r0,f);
708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727
  }
}

////////////////////////////////////////////////////////////////////////////////
void ConstraintSet::update_constraints(double dt)
{
  for ( int i = 0; i < constraint_list.size(); i++ )
  {
    constraint_list[i]->update(dt);
  }
}

////////////////////////////////////////////////////////////////////////////////
void ConstraintSet::setup(AtomSet& atoms)
{
  for ( int i = 0; i < constraint_list.size(); i++ )
  {
    constraint_list[i]->setup(atoms);
  }
}
Francois Gygi committed
728 729 730 731 732 733 734 735 736

////////////////////////////////////////////////////////////////////////////////
void ConstraintSet::reset(void)
{
  for ( int i = 0; i < constraint_list.size(); i++ )
    delete constraint_list[i];
  ndofs_ = 0;
  constraint_list.resize(0);
}