ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimSetup.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 780 by mmeineke, Mon Sep 22 21:23:25 2003 UTC vs.
Revision 1234 by tim, Fri Jun 4 03:15:31 2004 UTC

# Line 1 | Line 1
1   #include <algorithm>
2 < #include <cstdlib>
2 > #include <stdlib.h>
3   #include <iostream>
4 < #include <cmath>
4 > #include <math.h>
5   #include <string>
6   #include <sprng.h>
7
7   #include "SimSetup.hpp"
8   #include "ReadWrite.hpp"
9   #include "parse_me.h"
10   #include "Integrator.hpp"
11   #include "simError.h"
12 + #include "RigidBody.hpp"
13 + #include "OOPSEMinimizer.hpp"
14 + #include "ConstraintElement.hpp"
15 + #include "ConstraintPair.hpp"
16 + #include "ConstraintManager.hpp"
17  
18   #ifdef IS_MPI
19   #include "mpiBASS.h"
# Line 22 | Line 26
26   #define NVT_ENS        1
27   #define NPTi_ENS       2
28   #define NPTf_ENS       3
29 + #define NPTxyz_ENS     4
30  
26 #define FF_DUFF 0
27 #define FF_LJ   1
28 #define FF_EAM  2
31  
32 + #define FF_DUFF  0
33 + #define FF_LJ    1
34 + #define FF_EAM   2
35 + #define FF_H2O   3
36 +
37   using namespace std;
38  
39 + /**
40 + * Check whether dividend is divisble by divisor or not
41 + */
42 + bool isDivisible(double dividend, double divisor){
43 +  double tolerance = 0.000001;
44 +  double quotient;
45 +  double diff;
46 +  int intQuotient;
47 +  
48 +  quotient = dividend / divisor;
49 +
50 +  if (quotient < 0)
51 +    quotient = -quotient;
52 +
53 +  intQuotient = int (quotient + tolerance);
54 +
55 +  diff = fabs(fabs(dividend) - intQuotient  * fabs(divisor));
56 +
57 +  if (diff <= tolerance)
58 +    return true;
59 +  else
60 +    return false;  
61 + }
62 +
63   SimSetup::SimSetup(){
64 +  
65 +  initSuspend = false;
66    isInfoArray = 0;
67    nInfo = 1;
68  
# Line 52 | Line 85 | void SimSetup::setSimInfo(SimInfo* the_info, int theNi
85    info = the_info;
86    nInfo = theNinfo;
87    isInfoArray = 1;
88 +  initSuspend = true;
89   }
90  
91  
# Line 90 | Line 124 | void SimSetup::createSim(void){
124   #endif // is_mpi
125  
126   void SimSetup::createSim(void){
93  int i, j, k, globalAtomIndex;
127  
128    // gather all of the information from the Bass file
129  
# Line 106 | Line 139 | void SimSetup::createSim(void){
139  
140    // initialize the system coordinates
141  
142 <  if (!isInfoArray){
142 >  if ( !initSuspend ){
143      initSystemCoords();
144 +
145 +    if( !(globals->getUseInitTime()) )
146 +      info[0].currentTime = 0.0;
147    }  
148  
149    // make the output filenames
150  
151    makeOutNames();
152 <
117 <  // make the integrator
118 <
119 <  makeIntegrator();
120 <
152 >  
153   #ifdef IS_MPI
154    mpiSim->mpiRefresh();
155   #endif
# Line 125 | Line 157 | void SimSetup::createSim(void){
157    // initialize the Fortran
158  
159    initFortran();
160 +
161 +  //creat constraint manager
162 +  for(int i = 0; i < nInfo; i++)
163 +    info[i].consMan = new ConstraintManager(&info[i]);
164 +
165 +  if (globals->haveMinimizer())
166 +    // make minimizer
167 +    makeMinimizer();
168 +  else
169 +    // make the integrator
170 +    makeIntegrator();
171 +
172   }
173  
174  
175   void SimSetup::makeMolecules(void){
176 <  int k, l;
177 <  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
176 >  int i, j, k;
177 >  int exI, exJ, exK, exL, slI, slJ;
178 >  int tempI, tempJ, tempK, tempL;
179 >  int molI, globalID;
180 >  int stampID, atomOffset, rbOffset, groupOffset;
181    molInit molInfo;
182    DirectionalAtom* dAtom;
183 +  RigidBody* myRB;
184 +  StuntDouble* mySD;
185    LinkedAssign* extras;
186    LinkedAssign* current_extra;
187    AtomStamp* currentAtom;
188    BondStamp* currentBond;
189    BendStamp* currentBend;
190    TorsionStamp* currentTorsion;
191 +  RigidBodyStamp* currentRigidBody;
192 +  CutoffGroupStamp* currentCutoffGroup;
193 +  CutoffGroup* myCutoffGroup;
194 +  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
195 +  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
196  
197    bond_pair* theBonds;
198    bend_set* theBends;
199    torsion_set* theTorsions;
200  
201 +  set<int> skipList;
202  
203 +  double phi, theta, psi;
204 +  char* molName;
205 +  char rbName[100];
206 +
207 +  ConstraintPair* consPair; //constraint pair
208 +  ConstraintElement* consElement1;  //first element of constraint pair
209 +  ConstraintElement* consElement2;  //second element of constraint pair
210 +  int whichRigidBody;
211 +  int consAtomIndex;  //index of constraint atom in rigid body's atom array
212 +  vector<pair<int, int> > jointAtoms;
213 +  double bondLength2;
214    //init the forceField paramters
215  
216    the_ff->readParams();
217  
152
218    // init the atoms
219  
220 <  double ux, uy, uz, u, uSqr;
220 >  int nMembers, nNew, rb1, rb2;
221  
222    for (k = 0; k < nInfo; k++){
223      the_ff->setSimInfo(&(info[k]));
224  
225 + #ifdef IS_MPI
226 +    info[k].globalGroupMembership = new int[mpiSim->getNAtomsGlobal()];
227 +    for (i = 0; i < mpiSim->getNAtomsGlobal(); i++)
228 +      info[k].globalGroupMembership[i] = 0;
229 + #else
230 +    info[k].globalGroupMembership = new int[info[k].n_atoms];
231 +    for (i = 0; i < info[k].n_atoms; i++)
232 +      info[k].globalGroupMembership[i] = 0;
233 + #endif
234 +
235      atomOffset = 0;
236 <    excludeOffset = 0;
236 >    groupOffset = 0;
237 >
238      for (i = 0; i < info[k].n_mol; i++){
239        stampID = info[k].molecules[i].getStampID();
240 +      molName = comp_stamps[stampID]->getID();
241  
242        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
243        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
244        molInfo.nBends = comp_stamps[stampID]->getNBends();
245        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
246 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
246 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
247  
248 +      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
249 +      
250        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
172      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
173      molInfo.myBonds = new Bond * [molInfo.nBonds];
174      molInfo.myBends = new Bend * [molInfo.nBends];
175      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
251  
252 +      if (molInfo.nBonds > 0)
253 +        molInfo.myBonds = new Bond*[molInfo.nBonds];
254 +      else
255 +        molInfo.myBonds = NULL;
256 +
257 +      if (molInfo.nBends > 0)
258 +        molInfo.myBends = new Bend*[molInfo.nBends];
259 +      else
260 +        molInfo.myBends = NULL;
261 +
262 +      if (molInfo.nTorsions > 0)
263 +        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
264 +      else
265 +        molInfo.myTorsions = NULL;
266 +
267        theBonds = new bond_pair[molInfo.nBonds];
268        theBends = new bend_set[molInfo.nBends];
269        theTorsions = new torsion_set[molInfo.nTorsions];
270 <
270 >      
271        // make the Atoms
272  
273        for (j = 0; j < molInfo.nAtoms; j++){
274          currentAtom = comp_stamps[stampID]->getAtom(j);
275 +
276          if (currentAtom->haveOrientation()){
277            dAtom = new DirectionalAtom((j + atomOffset),
278                                        info[k].getConfiguration());
279            info[k].n_oriented++;
280            molInfo.myAtoms[j] = dAtom;
281  
282 <          ux = currentAtom->getOrntX();
283 <          uy = currentAtom->getOrntY();
284 <          uz = currentAtom->getOrntZ();
282 >          // Directional Atoms have standard unit vectors which are oriented
283 >          // in space using the three Euler angles.  We assume the standard
284 >          // unit vector was originally along the z axis below.
285  
286 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
286 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
287 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
288 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
289  
290 <          u = sqrt(uSqr);
291 <          ux = ux / u;
199 <          uy = uy / u;
200 <          uz = uz / u;
201 <
202 <          dAtom->setSUx(ux);
203 <          dAtom->setSUy(uy);
204 <          dAtom->setSUz(uz);
290 >          dAtom->setUnitFrameFromEuler(phi, theta, psi);
291 >            
292          }
293          else{
294 <          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
295 <                                               info[k].getConfiguration());
294 >
295 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
296 >
297          }
210        molInfo.myAtoms[j]->setType(currentAtom->getType());
298  
299 +        molInfo.myAtoms[j]->setType(currentAtom->getType());
300   #ifdef IS_MPI
301 <
214 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
215 <
301 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
302   #endif // is_mpi
303        }
304  
# Line 222 | Line 308 | void SimSetup::makeMolecules(void){
308          theBonds[j].a = currentBond->getA() + atomOffset;
309          theBonds[j].b = currentBond->getB() + atomOffset;
310  
311 <        exI = theBonds[j].a;
312 <        exJ = theBonds[j].b;
311 >        tempI = theBonds[j].a;
312 >        tempJ = theBonds[j].b;
313  
228        // exclude_I must always be the smaller of the pair
229        if (exI > exJ){
230          tempEx = exI;
231          exI = exJ;
232          exJ = tempEx;
233        }
314   #ifdef IS_MPI
315 <        tempEx = exI;
316 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
317 <        tempEx = exJ;
318 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
315 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
316 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
317 > #else
318 >        exI = tempI + 1;
319 >        exJ = tempJ + 1;
320 > #endif
321  
322 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
241 < #else  // isn't MPI
242 <
243 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
244 < #endif  //is_mpi
322 >        info[k].excludes->addPair(exI, exJ);
323        }
246      excludeOffset += molInfo.nBonds;
324  
325        //make the bends
326        for (j = 0; j < molInfo.nBends; j++){
# Line 293 | Line 370 | void SimSetup::makeMolecules(void){
370            }
371          }
372  
373 <        if (!theBends[j].isGhost){
374 <          exI = theBends[j].a;
375 <          exJ = theBends[j].c;
376 <        }
377 <        else{
301 <          exI = theBends[j].a;
302 <          exJ = theBends[j].b;
303 <        }
304 <
305 <        // exclude_I must always be the smaller of the pair
306 <        if (exI > exJ){
307 <          tempEx = exI;
308 <          exI = exJ;
309 <          exJ = tempEx;
310 <        }
373 >        if (theBends[j].isGhost) {
374 >          
375 >          tempI = theBends[j].a;
376 >          tempJ = theBends[j].b;
377 >          
378   #ifdef IS_MPI
379 <        tempEx = exI;
380 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
381 <        tempEx = exJ;
382 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
379 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
380 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
381 > #else
382 >          exI = tempI + 1;
383 >          exJ = tempJ + 1;
384 > #endif          
385 >          info[k].excludes->addPair(exI, exJ);
386  
387 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
388 < #else  // isn't MPI
389 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
390 < #endif  //is_mpi
387 >        } else {
388 >
389 >          tempI = theBends[j].a;
390 >          tempJ = theBends[j].b;
391 >          tempK = theBends[j].c;
392 >          
393 > #ifdef IS_MPI
394 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
395 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
396 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
397 > #else
398 >          exI = tempI + 1;
399 >          exJ = tempJ + 1;
400 >          exK = tempK + 1;
401 > #endif
402 >          
403 >          info[k].excludes->addPair(exI, exK);
404 >          info[k].excludes->addPair(exI, exJ);
405 >          info[k].excludes->addPair(exJ, exK);
406 >        }
407        }
322      excludeOffset += molInfo.nBends;
408  
409        for (j = 0; j < molInfo.nTorsions; j++){
410          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 328 | Line 413 | void SimSetup::makeMolecules(void){
413          theTorsions[j].c = currentTorsion->getC() + atomOffset;
414          theTorsions[j].d = currentTorsion->getD() + atomOffset;
415  
416 <        exI = theTorsions[j].a;
417 <        exJ = theTorsions[j].d;
416 >        tempI = theTorsions[j].a;      
417 >        tempJ = theTorsions[j].b;
418 >        tempK = theTorsions[j].c;
419 >        tempL = theTorsions[j].d;
420  
334        // exclude_I must always be the smaller of the pair
335        if (exI > exJ){
336          tempEx = exI;
337          exI = exJ;
338          exJ = tempEx;
339        }
421   #ifdef IS_MPI
422 <        tempEx = exI;
423 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
424 <        tempEx = exJ;
425 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
422 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
423 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
424 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
425 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
426 > #else
427 >        exI = tempI + 1;
428 >        exJ = tempJ + 1;
429 >        exK = tempK + 1;
430 >        exL = tempL + 1;
431 > #endif
432  
433 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
434 < #else  // isn't MPI
435 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
436 < #endif  //is_mpi
433 >        info[k].excludes->addPair(exI, exJ);
434 >        info[k].excludes->addPair(exI, exK);
435 >        info[k].excludes->addPair(exI, exL);        
436 >        info[k].excludes->addPair(exJ, exK);
437 >        info[k].excludes->addPair(exJ, exL);
438 >        info[k].excludes->addPair(exK, exL);
439        }
351      excludeOffset += molInfo.nTorsions;
440  
441 +      
442 +      molInfo.myRigidBodies.clear();
443 +      
444 +      for (j = 0; j < molInfo.nRigidBodies; j++){
445  
446 <      // send the arrays off to the forceField for init.
446 >        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
447 >        nMembers = currentRigidBody->getNMembers();
448 >
449 >        // Create the Rigid Body:
450 >
451 >        myRB = new RigidBody();
452 >
453 >        sprintf(rbName,"%s_RB_%d", molName, j);
454 >        myRB->setType(rbName);
455 >        
456 >        for (rb1 = 0; rb1 < nMembers; rb1++) {
457 >
458 >          // molI is atom numbering inside this molecule
459 >          molI = currentRigidBody->getMember(rb1);    
460  
461 +          // tempI is atom numbering on local processor
462 +          tempI = molI + atomOffset;
463 +
464 +          // currentAtom is the AtomStamp (which we need for
465 +          // rigid body reference positions)
466 +          currentAtom = comp_stamps[stampID]->getAtom(molI);
467 +
468 +          // When we add to the rigid body, add the atom itself and
469 +          // the stamp info:
470 +
471 +          myRB->addAtom(info[k].atoms[tempI], currentAtom);
472 +          
473 +          // Add this atom to the Skip List for the integrators
474 + #ifdef IS_MPI
475 +          slI = info[k].atoms[tempI]->getGlobalIndex();
476 + #else
477 +          slI = tempI;
478 + #endif
479 +          skipList.insert(slI);
480 +          
481 +        }
482 +        
483 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
484 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
485 +            
486 +            tempI = currentRigidBody->getMember(rb1);
487 +            tempJ = currentRigidBody->getMember(rb2);
488 +            
489 +            // Some explanation is required here.
490 +            // Fortran indexing starts at 1, while c indexing starts at 0
491 +            // Also, in parallel computations, the GlobalIndex is
492 +            // used for the exclude list:
493 +            
494 + #ifdef IS_MPI
495 +            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
496 +            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
497 + #else
498 +            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
499 +            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
500 + #endif
501 +            
502 +            info[k].excludes->addPair(exI, exJ);
503 +            
504 +          }
505 +        }
506 +
507 +        molInfo.myRigidBodies.push_back(myRB);
508 +        info[k].rigidBodies.push_back(myRB);
509 +      }
510 +      
511 +
512 +      //create cutoff group for molecule
513 +
514 +      cutoffAtomSet.clear();
515 +      molInfo.myCutoffGroups.clear();
516 +      
517 +      for (j = 0; j < nCutoffGroups; j++){
518 +
519 +        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
520 +        nMembers = currentCutoffGroup->getNMembers();
521 +
522 +        myCutoffGroup = new CutoffGroup();
523 +        
524 + #ifdef IS_MPI
525 +        myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]);
526 + #else
527 +        myCutoffGroup->setGlobalIndex(groupOffset);
528 + #endif
529 +        
530 +        for (int cg = 0; cg < nMembers; cg++) {
531 +
532 +          // molI is atom numbering inside this molecule
533 +          molI = currentCutoffGroup->getMember(cg);    
534 +
535 +          // tempI is atom numbering on local processor
536 +          tempI = molI + atomOffset;
537 +
538 + #ifdef IS_MPI
539 +          globalID = info[k].atoms[tempI]->getGlobalIndex();
540 +          info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset];
541 + #else
542 +          globalID = info[k].atoms[tempI]->getIndex();
543 +          info[k].globalGroupMembership[globalID] = groupOffset;
544 + #endif                    
545 +          myCutoffGroup->addAtom(info[k].atoms[tempI]);
546 +          cutoffAtomSet.insert(tempI);
547 +        }
548 +        
549 +        molInfo.myCutoffGroups.push_back(myCutoffGroup);
550 +        groupOffset++;
551 +
552 +      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
553 +      
554 +      
555 +      // create a cutoff group for every atom in current molecule which
556 +      // does not belong to cutoffgroup defined at mdl file
557 +      
558 +      for(j = 0; j < molInfo.nAtoms; j++){
559 +        
560 +        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
561 +          myCutoffGroup = new CutoffGroup();
562 +          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
563 +          
564 + #ifdef IS_MPI
565 +          myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]);
566 +          globalID = info[k].atoms[atomOffset + j]->getGlobalIndex();
567 +          info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset];
568 + #else
569 +          myCutoffGroup->setGlobalIndex(groupOffset);
570 +          globalID = info[k].atoms[atomOffset + j]->getIndex();
571 +          info[k].globalGroupMembership[globalID] = groupOffset;
572 + #endif
573 +          molInfo.myCutoffGroups.push_back(myCutoffGroup);
574 +          groupOffset++;
575 +        }          
576 +      }
577 +
578 +      // After this is all set up, scan through the atoms to
579 +      // see if they can be added to the integrableObjects:
580 +
581 +      molInfo.myIntegrableObjects.clear();
582 +      
583 +
584 +      for (j = 0; j < molInfo.nAtoms; j++){
585 +
586 + #ifdef IS_MPI
587 +        slJ = molInfo.myAtoms[j]->getGlobalIndex();
588 + #else
589 +        slJ = j+atomOffset;
590 + #endif
591 +
592 +        // if they aren't on the skip list, then they can be integrated
593 +
594 +        if (skipList.find(slJ) == skipList.end()) {
595 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
596 +          info[k].integrableObjects.push_back(mySD);
597 +          molInfo.myIntegrableObjects.push_back(mySD);
598 +        }
599 +      }
600 +
601 +      // all rigid bodies are integrated:
602 +
603 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
604 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
605 +        info[k].integrableObjects.push_back(mySD);      
606 +        molInfo.myIntegrableObjects.push_back(mySD);
607 +      }
608 +        
609 +      // send the arrays off to the forceField for init.
610 +      
611        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
612        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
613        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
# Line 360 | Line 615 | void SimSetup::makeMolecules(void){
615                                   theTorsions);
616  
617  
618 <      info[k].molecules[i].initialize(molInfo);
618 >      //creat ConstraintPair.
619 >      molInfo.myConstraintPairs.clear();
620 >      
621 >      for (j = 0; j < molInfo.nBonds; j++){
622  
623 +        //if bond is constrained bond, add it into constraint pair
624 +        if(molInfo.myBonds[j]->is_constrained()){
625  
626 +          //if both atoms are in the same rigid body, just skip it
627 +          currentBond = comp_stamps[stampID]->getBond(j);
628 +          
629 +          if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){
630 +
631 +            tempI = currentBond->getA() + atomOffset;
632 +            if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex))
633 +              consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
634 +            else
635 +               consElement1 = new ConstraintAtom(info[k].atoms[tempI]);      
636 +
637 +            tempJ =  currentBond->getB() + atomOffset;
638 +            if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex))
639 +              consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
640 +            else
641 +               consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);    
642 +
643 +            bondLength2 = molInfo.myBonds[j]->get_constraint()->get_dsqr();            
644 +            consPair = new DistanceConstraintPair(consElement1, consElement2, bondLength2);
645 +
646 +            molInfo.myConstraintPairs.push_back(consPair);
647 +          }
648 +        }//end if(molInfo.myBonds[j]->is_constrained())
649 +      }  
650 +      
651 +      //loop over rigid bodies, if two rigid bodies share same joint, creat a JointConstraintPair
652 +      for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){
653 +        for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){
654 +          
655 +          jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2);
656 +
657 +          for(size_t m = 0; m < jointAtoms.size(); m++){          
658 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first);
659 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second);
660 +
661 +            consPair = new JointConstraintPair(consElement1, consElement2);  
662 +            molInfo.myConstraintPairs.push_back(consPair);            
663 +          }
664 +
665 +        }
666 +      }
667 +      
668 +
669 +      info[k].molecules[i].initialize(molInfo);
670 +      
671 +      
672        atomOffset += molInfo.nAtoms;
673        delete[] theBonds;
674        delete[] theBends;
675        delete[] theTorsions;
676      }
677 +
678 +
679 +
680 + #ifdef IS_MPI    
681 +    // Since the globalGroupMembership has been zero filled and we've only
682 +    // poked values into the atoms we know, we can do an Allreduce
683 +    // to get the full globalGroupMembership array (We think).
684 +    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
685 +    // docs said we could.
686 +
687 +    int* ggMjunk = new int[mpiSim->getNAtomsGlobal()];    
688 +
689 +    MPI_Allreduce(info[k].globalGroupMembership,
690 +                  ggMjunk,
691 +                  mpiSim->getNAtomsGlobal(),
692 +                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
693 +
694 +    for (i = 0; i < mpiSim->getNAtomsGlobal(); i++)
695 +      info[k].globalGroupMembership[i] = ggMjunk[i];
696 +
697 +    delete[] ggMjunk;
698 +    
699 + #endif
700 +
701 +
702 +
703    }
704  
705   #ifdef IS_MPI
# Line 375 | Line 707 | void SimSetup::makeMolecules(void){
707    MPIcheckPoint();
708   #endif // is_mpi
709  
378  // clean up the forcefield
379
380  the_ff->calcRcut();
381  the_ff->cleanMe();
710   }
711  
712   void SimSetup::initFromBass(void){
# Line 551 | Line 879 | void SimSetup::gatherInfo(void){
879  
880  
881   void SimSetup::gatherInfo(void){
882 <  int i, j, k;
882 >  int i;
883  
884    ensembleCase = -1;
885    ffCase = -1;
# Line 579 | Line 907 | void SimSetup::gatherInfo(void){
907    else if (!strcasecmp(force_field, "EAM")){
908      ffCase = FF_EAM;
909    }
910 +  else if (!strcasecmp(force_field, "WATER")){
911 +    ffCase = FF_H2O;
912 +  }
913    else{
914      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
915              force_field);
# Line 602 | Line 933 | void SimSetup::gatherInfo(void){
933    else if (!strcasecmp(ensemble, "NPTf")){
934      ensembleCase = NPTf_ENS;
935    }
936 +  else if (!strcasecmp(ensemble, "NPTxyz")){
937 +    ensembleCase = NPTxyz_ENS;
938 +  }
939    else{
940      sprintf(painCave.errMsg,
941 <            "SimSetup Warning. Unrecognized Ensemble -> %s, "
942 <            "reverting to NVE for this simulation.\n",
941 >            "SimSetup Warning. Unrecognized Ensemble -> %s \n"
942 >            "\treverting to NVE for this simulation.\n",
943              ensemble);
944           painCave.isFatal = 0;
945           simError();
# Line 637 | Line 971 | void SimSetup::gatherInfo(void){
971        if (!the_components[i]->haveNMol()){
972          // we have a problem
973          sprintf(painCave.errMsg,
974 <                "SimSetup Error. No global NMol or component NMol"
975 <                " given. Cannot calculate the number of atoms.\n");
974 >                "SimSetup Error. No global NMol or component NMol given.\n"
975 >                "\tCannot calculate the number of atoms.\n");
976          painCave.isFatal = 1;
977          simError();
978        }
# Line 656 | Line 990 | void SimSetup::gatherInfo(void){
990              " Please give nMol in the components.\n");
991      painCave.isFatal = 1;
992      simError();
993 +  }
994 +
995 +  //check whether sample time, status time, thermal time and reset time are divisble by dt
996 +  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
997 +    sprintf(painCave.errMsg,
998 +            "Sample time is not divisible by dt.\n"
999 +            "\tThis will result in samples that are not uniformly\n"
1000 +            "\tdistributed in time.  If this is a problem, change\n"
1001 +            "\tyour sampleTime variable.\n");
1002 +    painCave.isFatal = 0;
1003 +    simError();    
1004 +  }
1005 +
1006 +  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
1007 +    sprintf(painCave.errMsg,
1008 +            "Status time is not divisible by dt.\n"
1009 +            "\tThis will result in status reports that are not uniformly\n"
1010 +            "\tdistributed in time.  If this is a problem, change \n"
1011 +            "\tyour statusTime variable.\n");
1012 +    painCave.isFatal = 0;
1013 +    simError();    
1014    }
1015  
1016 +  if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){
1017 +    sprintf(painCave.errMsg,
1018 +            "Thermal time is not divisible by dt.\n"
1019 +            "\tThis will result in thermalizations that are not uniformly\n"
1020 +            "\tdistributed in time.  If this is a problem, change \n"
1021 +            "\tyour thermalTime variable.\n");
1022 +    painCave.isFatal = 0;
1023 +    simError();    
1024 +  }  
1025 +
1026 +  if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){
1027 +    sprintf(painCave.errMsg,
1028 +            "Reset time is not divisible by dt.\n"
1029 +            "\tThis will result in integrator resets that are not uniformly\n"
1030 +            "\tdistributed in time.  If this is a problem, change\n"
1031 +            "\tyour resetTime variable.\n");
1032 +    painCave.isFatal = 0;
1033 +    simError();    
1034 +  }
1035 +
1036    // set the status, sample, and thermal kick times
1037  
1038    for (i = 0; i < nInfo; i++){
1039      if (globals->haveSampleTime()){
1040        info[i].sampleTime = globals->getSampleTime();
1041        info[i].statusTime = info[i].sampleTime;
667      info[i].thermalTime = info[i].sampleTime;
1042      }
1043      else{
1044        info[i].sampleTime = globals->getRunTime();
1045        info[i].statusTime = info[i].sampleTime;
672      info[i].thermalTime = info[i].sampleTime;
1046      }
1047  
1048      if (globals->haveStatusTime()){
# Line 678 | Line 1051 | void SimSetup::gatherInfo(void){
1051  
1052      if (globals->haveThermalTime()){
1053        info[i].thermalTime = globals->getThermalTime();
1054 +    } else {
1055 +      info[i].thermalTime = globals->getRunTime();
1056      }
1057  
1058      info[i].resetIntegrator = 0;
# Line 687 | Line 1062 | void SimSetup::gatherInfo(void){
1062      }
1063  
1064      // check for the temperature set flag
1065 <
1065 >    
1066      if (globals->haveTempSet())
1067        info[i].setTemp = globals->getTempSet();
1068  
1069 <    // get some of the tricky things that may still be in the globals
1069 >    // check for the extended State init
1070  
1071 <    double boxVector[3];
1072 <    if (globals->haveBox()){
698 <      boxVector[0] = globals->getBox();
699 <      boxVector[1] = globals->getBox();
700 <      boxVector[2] = globals->getBox();
701 <
702 <      info[i].setBox(boxVector);
703 <    }
704 <    else if (globals->haveDensity()){
705 <      double vol;
706 <      vol = (double) tot_nmol / globals->getDensity();
707 <      boxVector[0] = pow(vol, (1.0 / 3.0));
708 <      boxVector[1] = boxVector[0];
709 <      boxVector[2] = boxVector[0];
1071 >    info[i].useInitXSstate = globals->getUseInitXSstate();
1072 >    info[i].orthoTolerance = globals->getOrthoBoxTolerance();
1073  
1074 <      info[i].setBox(boxVector);
1074 >    // check for thermodynamic integration
1075 >    if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) {
1076 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1077 >        info[i].useSolidThermInt = globals->getUseSolidThermInt();
1078 >        info[i].thermIntLambda = globals->getThermIntLambda();
1079 >        info[i].thermIntK = globals->getThermIntK();
1080 >        
1081 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
1082 >        info[i].restraint = myRestraint;
1083 >      }
1084 >      else {
1085 >        sprintf(painCave.errMsg,
1086 >                "SimSetup Error:\n"
1087 >                "\tKeyword useSolidThermInt was set to 'true' but\n"
1088 >                "\tthermodynamicIntegrationLambda (and/or\n"
1089 >                "\tthermodynamicIntegrationK) was not specified.\n"
1090 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1091 >        painCave.isFatal = 1;
1092 >        simError();    
1093 >      }
1094      }
1095 <    else{
1096 <      if (!globals->haveBoxX()){
1097 <        sprintf(painCave.errMsg,
1098 <                "SimSetup error, no periodic BoxX size given.\n");
1099 <        painCave.isFatal = 1;
1100 <        simError();
1095 >    else if(globals->getUseLiquidThermInt()) {
1096 >      if (globals->getUseSolidThermInt()) {
1097 >        sprintf( painCave.errMsg,
1098 >                 "SimSetup Warning: It appears that you have both solid and\n"
1099 >                 "\tliquid thermodynamic integration activated in your .bass\n"
1100 >                 "\tfile. To avoid confusion, specify only one technique in\n"
1101 >                 "\tyour .bass file. Liquid-state thermodynamic integration\n"
1102 >                 "\twill be assumed for the current simulation. If this is not\n"
1103 >                 "\twhat you desire, set useSolidThermInt to 'true' and\n"
1104 >                 "\tuseLiquidThermInt to 'false' in your .bass file.\n");
1105 >        painCave.isFatal = 0;
1106 >        simError();
1107        }
1108 <      boxVector[0] = globals->getBoxX();
1109 <
1110 <      if (!globals->haveBoxY()){
1111 <        sprintf(painCave.errMsg,
724 <                "SimSetup error, no periodic BoxY size given.\n");
725 <        painCave.isFatal = 1;
726 <        simError();
1108 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1109 >        info[i].useLiquidThermInt = globals->getUseLiquidThermInt();
1110 >        info[i].thermIntLambda = globals->getThermIntLambda();
1111 >        info[i].thermIntK = globals->getThermIntK();
1112        }
1113 <      boxVector[1] = globals->getBoxY();
1114 <
1115 <      if (!globals->haveBoxZ()){
1116 <        sprintf(painCave.errMsg,
1117 <                "SimSetup error, no periodic BoxZ size given.\n");
1118 <        painCave.isFatal = 1;
1119 <        simError();
1113 >      else {
1114 >        sprintf(painCave.errMsg,
1115 >                "SimSetup Error:\n"
1116 >                "\tKeyword useLiquidThermInt was set to 'true' but\n"
1117 >                "\tthermodynamicIntegrationLambda (and/or\n"
1118 >                "\tthermodynamicIntegrationK) was not specified.\n"
1119 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1120 >        painCave.isFatal = 1;
1121 >        simError();    
1122        }
736      boxVector[2] = globals->getBoxZ();
737
738      info[i].setBox(boxVector);
1123      }
1124 +    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
1125 +        sprintf(painCave.errMsg,
1126 +                "SimSetup Warning: If you want to use Thermodynamic\n"
1127 +                "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n"
1128 +                "\t'true' in your .bass file.  These keywords are set to\n"
1129 +                "\t'false' by default, so your lambda and/or k values are\n"
1130 +                "\tbeing ignored.\n");
1131 +        painCave.isFatal = 0;
1132 +        simError();  
1133 +    }
1134    }
1135 <
1135 >  
1136    //setup seed for random number generator
1137    int seedValue;
1138  
# Line 778 | Line 1172 | void SimSetup::gatherInfo(void){
1172    for (int i = 0; i < nInfo; i++){
1173      info[i].setSeed(seedValue);
1174    }
1175 <
1175 >  
1176   #ifdef IS_MPI
1177 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
1177 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
1178    MPIcheckPoint();
1179   #endif // is_mpi
1180   }
# Line 789 | Line 1183 | void SimSetup::finalInfoCheck(void){
1183   void SimSetup::finalInfoCheck(void){
1184    int index;
1185    int usesDipoles;
1186 +  int usesCharges;
1187    int i;
1188  
1189    for (i = 0; i < nInfo; i++){
# Line 800 | Line 1195 | void SimSetup::finalInfoCheck(void){
1195        usesDipoles = (info[i].atoms[index])->hasDipole();
1196        index++;
1197      }
1198 <
1198 >    index = 0;
1199 >    usesCharges = 0;
1200 >    while ((index < info[i].n_atoms) && !usesCharges){
1201 >      usesCharges= (info[i].atoms[index])->hasCharge();
1202 >      index++;
1203 >    }
1204   #ifdef IS_MPI
1205      int myUse = usesDipoles;
1206      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1207   #endif //is_mpi
1208  
1209 <    double theEcr, theEst;
1209 >    double theRcut, theRsw;
1210  
1211 +    if (globals->haveRcut()) {
1212 +      theRcut = globals->getRcut();
1213 +
1214 +      if (globals->haveRsw())
1215 +        theRsw = globals->getRsw();
1216 +      else
1217 +        theRsw = theRcut;
1218 +      
1219 +      info[i].setDefaultRcut(theRcut, theRsw);
1220 +
1221 +    } else {
1222 +      
1223 +      the_ff->calcRcut();
1224 +      theRcut = info[i].getRcut();
1225 +
1226 +      if (globals->haveRsw())
1227 +        theRsw = globals->getRsw();
1228 +      else
1229 +        theRsw = theRcut;
1230 +      
1231 +      info[i].setDefaultRcut(theRcut, theRsw);
1232 +    }
1233 +
1234      if (globals->getUseRF()){
1235        info[i].useReactionField = 1;
1236 <
1237 <      if (!globals->haveECR()){
1236 >      
1237 >      if (!globals->haveRcut()){
1238          sprintf(painCave.errMsg,
1239 <                "SimSetup Warning: using default value of 1/2 the smallest "
1240 <                "box length for the electrostaticCutoffRadius.\n"
1241 <                "I hope you have a very fast processor!\n");
1239 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1240 >                "\tOOPSE will use a default value of 15.0 angstroms"
1241 >                "\tfor the cutoffRadius.\n");
1242          painCave.isFatal = 0;
1243          simError();
1244 <        double smallest;
822 <        smallest = info[i].boxL[0];
823 <        if (info[i].boxL[1] <= smallest)
824 <          smallest = info[i].boxL[1];
825 <        if (info[i].boxL[2] <= smallest)
826 <          smallest = info[i].boxL[2];
827 <        theEcr = 0.5 * smallest;
1244 >        theRcut = 15.0;
1245        }
1246        else{
1247 <        theEcr = globals->getECR();
1247 >        theRcut = globals->getRcut();
1248        }
1249  
1250 <      if (!globals->haveEST()){
1250 >      if (!globals->haveRsw()){
1251          sprintf(painCave.errMsg,
1252 <                "SimSetup Warning: using default value of 0.05 * the "
1253 <                "electrostaticCutoffRadius for the electrostaticSkinThickness\n");
1252 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1253 >                "\tOOPSE will use a default value of\n"
1254 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1255          painCave.isFatal = 0;
1256          simError();
1257 <        theEst = 0.05 * theEcr;
1257 >        theRsw = 0.95 * theRcut;
1258        }
1259        else{
1260 <        theEst = globals->getEST();
1260 >        theRsw = globals->getRsw();
1261        }
1262  
1263 <      info[i].setEcr(theEcr, theEst);
1263 >      info[i].setDefaultRcut(theRcut, theRsw);
1264  
1265        if (!globals->haveDielectric()){
1266          sprintf(painCave.errMsg,
1267 <                "SimSetup Error: You are trying to use Reaction Field without"
1268 <                "setting a dielectric constant!\n");
1267 >                "SimSetup Error: No Dielectric constant was set.\n"
1268 >                "\tYou are trying to use Reaction Field without"
1269 >                "\tsetting a dielectric constant!\n");
1270          painCave.isFatal = 1;
1271          simError();
1272        }
1273        info[i].dielectric = globals->getDielectric();
1274      }
1275      else{
1276 <      if (usesDipoles){
1277 <        if (!globals->haveECR()){
1276 >      if (usesDipoles || usesCharges){
1277 >
1278 >        if (!globals->haveRcut()){
1279            sprintf(painCave.errMsg,
1280 <                  "SimSetup Warning: using default value of 1/2 the smallest "
1281 <                  "box length for the electrostaticCutoffRadius.\n"
1282 <                  "I hope you have a very fast processor!\n");
1280 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1281 >                  "\tOOPSE will use a default value of 15.0 angstroms"
1282 >                  "\tfor the cutoffRadius.\n");
1283            painCave.isFatal = 0;
1284            simError();
1285 <          double smallest;
1286 <          smallest = info[i].boxL[0];
867 <          if (info[i].boxL[1] <= smallest)
868 <            smallest = info[i].boxL[1];
869 <          if (info[i].boxL[2] <= smallest)
870 <            smallest = info[i].boxL[2];
871 <          theEcr = 0.5 * smallest;
872 <        }
1285 >          theRcut = 15.0;
1286 >      }
1287          else{
1288 <          theEcr = globals->getECR();
1288 >          theRcut = globals->getRcut();
1289          }
1290 <
1291 <        if (!globals->haveEST()){
1290 >        
1291 >        if (!globals->haveRsw()){
1292            sprintf(painCave.errMsg,
1293 <                  "SimSetup Warning: using default value of 0.05 * the "
1294 <                  "electrostaticCutoffRadius for the "
1295 <                  "electrostaticSkinThickness\n");
1293 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1294 >                  "\tOOPSE will use a default value of\n"
1295 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1296            painCave.isFatal = 0;
1297            simError();
1298 <          theEst = 0.05 * theEcr;
1298 >          theRsw = 0.95 * theRcut;
1299          }
1300          else{
1301 <          theEst = globals->getEST();
1301 >          theRsw = globals->getRsw();
1302          }
1303 <
1304 <        info[i].setEcr(theEcr, theEst);
1303 >        
1304 >        info[i].setDefaultRcut(theRcut, theRsw);
1305 >        
1306        }
1307      }
1308    }
894
1309   #ifdef IS_MPI
1310    strcpy(checkPointMsg, "post processing checks out");
1311    MPIcheckPoint();
1312   #endif // is_mpi
899 }
1313  
1314 +  // clean up the forcefield
1315 +  the_ff->cleanMe();
1316 + }
1317 +  
1318   void SimSetup::initSystemCoords(void){
1319    int i;
1320  
# Line 914 | Line 1331 | void SimSetup::initSystemCoords(void){
1331      if (worldRank == 0){
1332   #endif //is_mpi
1333        inName = globals->getInitialConfig();
917      double* tempDouble = new double[1000000];
1334        fileInit = new InitializeFromFile(inName);
1335   #ifdef IS_MPI
1336      }
# Line 926 | Line 1342 | void SimSetup::initSystemCoords(void){
1342      delete fileInit;
1343    }
1344    else{
1345 < #ifdef IS_MPI
930 <
1345 >    
1346      // no init from bass
1347 <
1347 >    
1348      sprintf(painCave.errMsg,
1349 <            "Cannot intialize a parallel simulation without an initial configuration file.\n");
1350 <    painCave.isFatal;
1349 >            "Cannot intialize a simulation without an initial configuration file.\n");
1350 >    painCave.isFatal = 1;;
1351      simError();
1352 <
938 < #else
939 <
940 <    initFromBass();
941 <
942 <
943 < #endif
1352 >    
1353    }
1354  
1355   #ifdef IS_MPI
# Line 1034 | Line 1443 | void SimSetup::makeOutNames(void){
1443          }
1444        }
1445  
1446 +      strcpy(info[k].rawPotName, inFileName);
1447 +      nameLength = strlen(info[k].rawPotName);
1448 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1449 +      if (!strcmp(endTest, ".bass")){
1450 +        strcpy(endTest, ".raw");
1451 +      }
1452 +      else if (!strcmp(endTest, ".BASS")){
1453 +        strcpy(endTest, ".raw");
1454 +      }
1455 +      else{
1456 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1457 +        if (!strcmp(endTest, ".bss")){
1458 +          strcpy(endTest, ".raw");
1459 +        }
1460 +        else if (!strcmp(endTest, ".mdl")){
1461 +          strcpy(endTest, ".raw");
1462 +        }
1463 +        else{
1464 +          strcat(info[k].rawPotName, ".raw");
1465 +        }
1466 +      }
1467 +
1468   #ifdef IS_MPI
1469  
1470      }
# Line 1094 | Line 1525 | void SimSetup::createFF(void){
1525        the_ff = new EAM_FF();
1526        break;
1527  
1528 +    case FF_H2O:
1529 +      the_ff = new WATER();
1530 +      break;
1531 +
1532      default:
1533        sprintf(painCave.errMsg,
1534                "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1114 | Line 1549 | void SimSetup::compList(void){
1549    LinkedMolStamp* headStamp = new LinkedMolStamp();
1550    LinkedMolStamp* currentStamp = NULL;
1551    comp_stamps = new MoleculeStamp * [n_components];
1552 +  bool haveCutoffGroups;
1553  
1554 +  haveCutoffGroups = false;
1555 +  
1556    // make an array of molecule stamps that match the components used.
1557    // also extract the used stamps out into a separate linked list
1558  
# Line 1149 | Line 1587 | void SimSetup::compList(void){
1587        headStamp->add(currentStamp);
1588        comp_stamps[i] = headStamp->match(id);
1589      }
1590 +
1591 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1592 +      haveCutoffGroups = true;    
1593    }
1594 +    
1595 +  for (i = 0; i < nInfo; i++)
1596 +    info[i].haveCutoffGroups = haveCutoffGroups;
1597  
1598   #ifdef IS_MPI
1599    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1158 | Line 1602 | void SimSetup::calcSysValues(void){
1602   }
1603  
1604   void SimSetup::calcSysValues(void){
1605 <  int i, j, k;
1605 >  int i, j;
1606 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1607  
1608    int* molMembershipArray;
1609 +  CutoffGroupStamp* cg;
1610  
1611    tot_atoms = 0;
1612    tot_bonds = 0;
1613    tot_bends = 0;
1614    tot_torsions = 0;
1615 +  tot_rigid = 0;
1616 +  tot_groups = 0;
1617    for (i = 0; i < n_components; i++){
1618      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1619      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1620      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1621      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1622 <  }
1622 >    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1623  
1624 +    ncutgroups = comp_stamps[i]->getNCutoffGroups();
1625 +    atomsingroups = 0;
1626 +    for (j=0; j < ncutgroups; j++) {
1627 +      cg = comp_stamps[i]->getCutoffGroup(j);
1628 +      atomsingroups += cg->getNMembers();
1629 +    }
1630 +    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
1631 +    tot_groups += components_nmol[i] * ngroupsinstamp;    
1632 +  }
1633 +  
1634    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1635    molMembershipArray = new int[tot_atoms];
1636  
# Line 1183 | Line 1641 | void SimSetup::calcSysValues(void){
1641      info[i].n_torsions = tot_torsions;
1642      info[i].n_SRI = tot_SRI;
1643      info[i].n_mol = tot_nmol;
1644 <
1644 >    info[i].ngroup = tot_groups;
1645      info[i].molMembershipArray = molMembershipArray;
1646    }
1647   }
# Line 1194 | Line 1652 | void SimSetup::mpiMolDivide(void){
1652    int i, j, k;
1653    int localMol, allMol;
1654    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1655 +  int local_rigid, local_groups;
1656 +  vector<int> globalMolIndex;
1657 +  int ncutgroups, atomsingroups, ngroupsinstamp;
1658 +  CutoffGroupStamp* cg;
1659  
1660    mpiSim = new mpiSimulation(info);
1661  
1662 <  globalIndex = mpiSim->divideLabor();
1662 >  mpiSim->divideLabor();
1663 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1664 >  globalGroupIndex = mpiSim->getGlobalGroupIndex();
1665 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1666  
1667    // set up the local variables
1668  
# Line 1210 | Line 1675 | void SimSetup::mpiMolDivide(void){
1675    local_bonds = 0;
1676    local_bends = 0;
1677    local_torsions = 0;
1678 <  globalAtomIndex = 0;
1678 >  local_rigid = 0;
1679 >  local_groups = 0;
1680 >  globalAtomCounter = 0;
1681  
1215
1682    for (i = 0; i < n_components; i++){
1683      for (j = 0; j < components_nmol[i]; j++){
1684        if (mol2proc[allMol] == worldRank){
# Line 1220 | Line 1686 | void SimSetup::mpiMolDivide(void){
1686          local_bonds += comp_stamps[i]->getNBonds();
1687          local_bends += comp_stamps[i]->getNBends();
1688          local_torsions += comp_stamps[i]->getNTorsions();
1689 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1690 +
1691 +        ncutgroups = comp_stamps[i]->getNCutoffGroups();
1692 +        atomsingroups = 0;
1693 +        for (k=0; k < ncutgroups; k++) {
1694 +          cg = comp_stamps[i]->getCutoffGroup(k);
1695 +          atomsingroups += cg->getNMembers();
1696 +        }
1697 +        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +
1698 +          ncutgroups;
1699 +        local_groups += ngroupsinstamp;    
1700 +
1701          localMol++;
1702        }      
1703        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1704 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1705 <        globalAtomIndex++;
1704 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1705 >        globalAtomCounter++;
1706        }
1707  
1708        allMol++;
# Line 1232 | Line 1710 | void SimSetup::mpiMolDivide(void){
1710    }
1711    local_SRI = local_bonds + local_bends + local_torsions;
1712  
1713 <  info[0].n_atoms = mpiSim->getMyNlocal();  
1714 <
1713 >  info[0].n_atoms = mpiSim->getNAtomsLocal();  
1714 >  
1715    if (local_atoms != info[0].n_atoms){
1716      sprintf(painCave.errMsg,
1717 <            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1718 <            " localAtom (%d) are not equal.\n",
1717 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
1718 >            "\tlocalAtom (%d) are not equal.\n",
1719              info[0].n_atoms, local_atoms);
1720      painCave.isFatal = 1;
1721      simError();
1722    }
1723  
1724 +  info[0].ngroup = mpiSim->getNGroupsLocal();  
1725 +  if (local_groups != info[0].ngroup){
1726 +    sprintf(painCave.errMsg,
1727 +            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
1728 +            "\tlocalGroups (%d) are not equal.\n",
1729 +            info[0].ngroup, local_groups);
1730 +    painCave.isFatal = 1;
1731 +    simError();
1732 +  }
1733 +  
1734    info[0].n_bonds = local_bonds;
1735    info[0].n_bends = local_bends;
1736    info[0].n_torsions = local_torsions;
# Line 1257 | Line 1745 | void SimSetup::makeSysArrays(void){
1745  
1746  
1747   void SimSetup::makeSysArrays(void){
1748 <  int i, j, k, l;
1748 >
1749 > #ifndef IS_MPI
1750 >  int k, j;
1751 > #endif // is_mpi
1752 >  int i, l;
1753  
1754    Atom** the_atoms;
1755    Molecule* the_molecules;
1264  Exclude** the_excludes;
1756  
1266
1757    for (l = 0; l < nInfo; l++){
1758      // create the atom and short range interaction arrays
1759  
# Line 1277 | Line 1767 | void SimSetup::makeSysArrays(void){
1767  
1768  
1769      molIndex = 0;
1770 <    for (i = 0; i < mpiSim->getTotNmol(); i++){
1770 >    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
1771        if (mol2proc[i] == worldRank){
1772          the_molecules[molIndex].setStampID(molCompType[i]);
1773          the_molecules[molIndex].setMyIndex(molIndex);
# Line 1289 | Line 1779 | void SimSetup::makeSysArrays(void){
1779   #else // is_mpi
1780  
1781      molIndex = 0;
1782 <    globalAtomIndex = 0;
1782 >    globalAtomCounter = 0;
1783      for (i = 0; i < n_components; i++){
1784        for (j = 0; j < components_nmol[i]; j++){
1785          the_molecules[molIndex].setStampID(i);
1786          the_molecules[molIndex].setMyIndex(molIndex);
1787          the_molecules[molIndex].setGlobalIndex(molIndex);
1788          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1789 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1790 <          globalAtomIndex++;
1789 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1790 >          globalAtomCounter++;
1791          }
1792          molIndex++;
1793        }
# Line 1306 | Line 1796 | void SimSetup::makeSysArrays(void){
1796  
1797   #endif // is_mpi
1798  
1799 <
1800 <    if (info[l].n_SRI){
1801 <      Exclude::createArray(info[l].n_SRI);
1312 <      the_excludes = new Exclude * [info[l].n_SRI];
1313 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1314 <        the_excludes[ex] = new Exclude(ex);
1315 <      }
1316 <      info[l].globalExcludes = new int;
1317 <      info[l].n_exclude = info[l].n_SRI;
1318 <    }
1319 <    else{
1320 <      Exclude::createArray(1);
1321 <      the_excludes = new Exclude * ;
1322 <      the_excludes[0] = new Exclude(0);
1323 <      the_excludes[0]->setPair(0, 0);
1324 <      info[l].globalExcludes = new int;
1325 <      info[l].globalExcludes[0] = 0;
1326 <      info[l].n_exclude = 0;
1327 <    }
1328 <
1799 >    info[l].globalExcludes = new int;
1800 >    info[l].globalExcludes[0] = 0;
1801 >    
1802      // set the arrays into the SimInfo object
1803  
1804      info[l].atoms = the_atoms;
1805      info[l].molecules = the_molecules;
1806      info[l].nGlobalExcludes = 0;
1807 <    info[l].excludes = the_excludes;
1335 <
1807 >    
1808      the_ff->setSimInfo(info);
1809    }
1810   }
# Line 1340 | Line 1812 | void SimSetup::makeIntegrator(void){
1812   void SimSetup::makeIntegrator(void){
1813    int k;
1814  
1815 +  NVE<RealIntegrator>* myNVE = NULL;
1816    NVT<RealIntegrator>* myNVT = NULL;
1817    NPTi<NPT<RealIntegrator> >* myNPTi = NULL;
1818    NPTf<NPT<RealIntegrator> >* myNPTf = NULL;
1819 +  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL;
1820    
1821    for (k = 0; k < nInfo; k++){
1822      switch (ensembleCase){
1823        case NVE_ENS:
1824          if (globals->haveZconstraints()){
1825            setupZConstraint(info[k]);
1826 <          new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff);
1826 >          myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff);
1827          }
1828 <        else
1829 <          new NVE<RealIntegrator>(&(info[k]), the_ff);
1828 >        else{
1829 >          myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff);
1830 >        }
1831 >        
1832 >        info->the_integrator = myNVE;
1833          break;
1834  
1835        case NVT_ENS:
# Line 1370 | Line 1847 | void SimSetup::makeIntegrator(void){
1847          else{
1848            sprintf(painCave.errMsg,
1849                    "SimSetup error: If you use the NVT\n"
1850 <                  "    ensemble, you must set tauThermostat.\n");
1850 >                  "\tensemble, you must set tauThermostat.\n");
1851            painCave.isFatal = 1;
1852            simError();
1853          }
1854 +
1855 +        info->the_integrator = myNVT;
1856          break;
1857  
1858        case NPTi_ENS:
# Line 1391 | Line 1870 | void SimSetup::makeIntegrator(void){
1870          else{
1871            sprintf(painCave.errMsg,
1872                    "SimSetup error: If you use a constant pressure\n"
1873 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1873 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1874            painCave.isFatal = 1;
1875            simError();
1876          }
# Line 1401 | Line 1880 | void SimSetup::makeIntegrator(void){
1880          else{
1881            sprintf(painCave.errMsg,
1882                    "SimSetup error: If you use an NPT\n"
1883 <                  "    ensemble, you must set tauThermostat.\n");
1883 >                  "\tensemble, you must set tauThermostat.\n");
1884            painCave.isFatal = 1;
1885            simError();
1886          }
# Line 1411 | Line 1890 | void SimSetup::makeIntegrator(void){
1890          else{
1891            sprintf(painCave.errMsg,
1892                    "SimSetup error: If you use an NPT\n"
1893 <                  "    ensemble, you must set tauBarostat.\n");
1893 >                  "\tensemble, you must set tauBarostat.\n");
1894            painCave.isFatal = 1;
1895            simError();
1896          }
1897 +
1898 +        info->the_integrator = myNPTi;
1899          break;
1900  
1901        case NPTf_ENS:
# Line 1432 | Line 1913 | void SimSetup::makeIntegrator(void){
1913          else{
1914            sprintf(painCave.errMsg,
1915                    "SimSetup error: If you use a constant pressure\n"
1916 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1916 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1917            painCave.isFatal = 1;
1918            simError();
1919          }    
1920  
1921          if (globals->haveTauThermostat())
1922            myNPTf->setTauThermostat(globals->getTauThermostat());
1923 +
1924          else{
1925            sprintf(painCave.errMsg,
1926                    "SimSetup error: If you use an NPT\n"
1927 <                  "    ensemble, you must set tauThermostat.\n");
1927 >                  "\tensemble, you must set tauThermostat.\n");
1928            painCave.isFatal = 1;
1929            simError();
1930          }
1931  
1932          if (globals->haveTauBarostat())
1933            myNPTf->setTauBarostat(globals->getTauBarostat());
1934 +
1935          else{
1936            sprintf(painCave.errMsg,
1937                    "SimSetup error: If you use an NPT\n"
1938 <                  "    ensemble, you must set tauBarostat.\n");
1938 >                  "\tensemble, you must set tauBarostat.\n");
1939 >          painCave.isFatal = 1;
1940 >          simError();
1941 >        }
1942 >
1943 >        info->the_integrator = myNPTf;
1944 >        break;
1945 >
1946 >      case NPTxyz_ENS:
1947 >        if (globals->haveZconstraints()){
1948 >          setupZConstraint(info[k]);
1949 >          myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1950 >        }
1951 >        else
1952 >          myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff);
1953 >
1954 >        myNPTxyz->setTargetTemp(globals->getTargetTemp());
1955 >
1956 >        if (globals->haveTargetPressure())
1957 >          myNPTxyz->setTargetPressure(globals->getTargetPressure());
1958 >        else{
1959 >          sprintf(painCave.errMsg,
1960 >                  "SimSetup error: If you use a constant pressure\n"
1961 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1962            painCave.isFatal = 1;
1963            simError();
1964 +        }    
1965 +
1966 +        if (globals->haveTauThermostat())
1967 +          myNPTxyz->setTauThermostat(globals->getTauThermostat());
1968 +        else{
1969 +          sprintf(painCave.errMsg,
1970 +                  "SimSetup error: If you use an NPT\n"
1971 +                  "\tensemble, you must set tauThermostat.\n");
1972 +          painCave.isFatal = 1;
1973 +          simError();
1974          }
1975 +
1976 +        if (globals->haveTauBarostat())
1977 +          myNPTxyz->setTauBarostat(globals->getTauBarostat());
1978 +        else{
1979 +          sprintf(painCave.errMsg,
1980 +                  "SimSetup error: If you use an NPT\n"
1981 +                  "\tensemble, you must set tauBarostat.\n");
1982 +          painCave.isFatal = 1;
1983 +          simError();
1984 +        }
1985 +
1986 +        info->the_integrator = myNPTxyz;
1987          break;
1988  
1989        default:
# Line 1503 | Line 2031 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2031    }
2032    else{
2033      sprintf(painCave.errMsg,
2034 <            "ZConstraint error: If you use an ZConstraint\n"
2035 <            " , you must set sample time.\n");
2034 >            "ZConstraint error: If you use a ZConstraint,\n"
2035 >            "\tyou must set zconsTime.\n");
2036      painCave.isFatal = 1;
2037      simError();
2038    }
# Line 1519 | Line 2047 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2047    else{
2048      double defaultZConsTol = 0.01;
2049      sprintf(painCave.errMsg,
2050 <            "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
2051 <            " , default value %f is used.\n",
2050 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
2051 >            "\tOOPSE will use a default value of %f.\n"
2052 >            "\tTo set the tolerance, use the zconsTol variable.\n",
2053              defaultZConsTol);
2054      painCave.isFatal = 0;
2055      simError();      
# Line 1538 | Line 2067 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2067    }
2068    else{
2069      sprintf(painCave.errMsg,
2070 <            "ZConstraint Warning: User does not set force Subtraction policy, "
2071 <            "PolicyByMass is used\n");
2070 >            "ZConstraint Warning: No force subtraction policy was set.\n"
2071 >            "\tOOPSE will use PolicyByMass.\n"
2072 >            "\tTo set the policy, use the zconsForcePolicy variable.\n");
2073      painCave.isFatal = 0;
2074      simError();
2075      zconsForcePolicy->setData("BYMASS");
# Line 1547 | Line 2077 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2077  
2078    theInfo.addProperty(zconsForcePolicy);
2079  
2080 +  //set zcons gap
2081 +  DoubleData* zconsGap = new DoubleData();
2082 +  zconsGap->setID(ZCONSGAP_ID);
2083 +
2084 +  if (globals->haveZConsGap()){
2085 +    zconsGap->setData(globals->getZconsGap());
2086 +    theInfo.addProperty(zconsGap);  
2087 +  }
2088 +
2089 +  //set zcons fixtime
2090 +  DoubleData* zconsFixtime = new DoubleData();
2091 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
2092 +
2093 +  if (globals->haveZConsFixTime()){
2094 +    zconsFixtime->setData(globals->getZconsFixtime());
2095 +    theInfo.addProperty(zconsFixtime);  
2096 +  }
2097 +
2098 +  //set zconsUsingSMD
2099 +  IntData* zconsUsingSMD = new IntData();
2100 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
2101 +
2102 +  if (globals->haveZConsUsingSMD()){
2103 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
2104 +    theInfo.addProperty(zconsUsingSMD);  
2105 +  }
2106 +
2107    //Determine the name of ouput file and add it into SimInfo's property list
2108    //Be careful, do not use inFileName, since it is a pointer which
2109    //point to a string at master node, and slave nodes do not contain that string
# Line 1576 | Line 2133 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2133      tempParaItem.zPos = zconStamp[i]->getZpos();
2134      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
2135      tempParaItem.kRatio = zconStamp[i]->getKratio();
2136 <
2136 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
2137 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
2138      zconsParaData->addItem(tempParaItem);
2139    }
2140  
2141    //check the uniqueness of index  
2142    if(!zconsParaData->isIndexUnique()){
2143      sprintf(painCave.errMsg,
2144 <            "ZConstraint Error: molIndex is not unique\n");
2144 >            "ZConstraint Error: molIndex is not unique!\n");
2145      painCave.isFatal = 1;
2146      simError();
2147    }
# Line 1594 | Line 2152 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2152    //push data into siminfo, therefore, we can retrieve later
2153    theInfo.addProperty(zconsParaData);
2154   }
2155 +
2156 + void SimSetup::makeMinimizer(){
2157 +
2158 +  OOPSEMinimizer* myOOPSEMinimizer;
2159 +  MinimizerParameterSet* param;
2160 +  char minimizerName[100];
2161 +  
2162 +  for (int i = 0; i < nInfo; i++){
2163 +    
2164 +    //prepare parameter set for minimizer
2165 +    param = new MinimizerParameterSet();
2166 +    param->setDefaultParameter();
2167 +
2168 +    if (globals->haveMinimizer()){
2169 +      param->setFTol(globals->getMinFTol());
2170 +    }
2171 +
2172 +    if (globals->haveMinGTol()){
2173 +      param->setGTol(globals->getMinGTol());
2174 +    }
2175 +
2176 +    if (globals->haveMinMaxIter()){
2177 +      param->setMaxIteration(globals->getMinMaxIter());
2178 +    }
2179 +
2180 +    if (globals->haveMinWriteFrq()){
2181 +      param->setMaxIteration(globals->getMinMaxIter());
2182 +    }
2183 +
2184 +    if (globals->haveMinWriteFrq()){
2185 +      param->setWriteFrq(globals->getMinWriteFrq());
2186 +    }
2187 +    
2188 +    if (globals->haveMinStepSize()){
2189 +      param->setStepSize(globals->getMinStepSize());
2190 +    }
2191 +
2192 +    if (globals->haveMinLSMaxIter()){
2193 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
2194 +    }    
2195 +
2196 +    if (globals->haveMinLSTol()){
2197 +      param->setLineSearchTol(globals->getMinLSTol());
2198 +    }    
2199 +
2200 +    strcpy(minimizerName, globals->getMinimizer());
2201 +
2202 +    if (!strcasecmp(minimizerName, "CG")){
2203 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
2204 +    }
2205 +    else if (!strcasecmp(minimizerName, "SD")){
2206 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
2207 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
2208 +    }
2209 +    else{
2210 +          sprintf(painCave.errMsg,
2211 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
2212 +          painCave.isFatal = 0;
2213 +          simError();
2214 +
2215 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
2216 +    }
2217 +     info[i].the_integrator = myOOPSEMinimizer;
2218 +
2219 +     //store the minimizer into simInfo
2220 +     info[i].the_minimizer = myOOPSEMinimizer;
2221 +     info[i].has_minimizer = true;
2222 +  }
2223 +
2224 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines