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 965 by gezelter, Mon Jan 19 21:17:39 2004 UTC vs.
Revision 1212 by chrisfen, Tue Jun 1 17:15:43 2004 UTC

# Line 9 | Line 9
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  
17   #ifdef IS_MPI
18   #include "mpiBASS.h"
# Line 24 | Line 28
28   #define NPTxyz_ENS     4
29  
30  
31 < #define FF_DUFF 0
32 < #define FF_LJ   1
33 < #define FF_EAM  2
31 > #define FF_DUFF  0
32 > #define FF_LJ    1
33 > #define FF_EAM   2
34 > #define FF_H2O   3
35  
36   using namespace std;
37  
# Line 143 | Line 148 | void SimSetup::createSim(void){
148    // make the output filenames
149  
150    makeOutNames();
151 <
147 <  // make the integrator
148 <
149 <  makeIntegrator();
150 <
151 >  
152   #ifdef IS_MPI
153    mpiSim->mpiRefresh();
154   #endif
# Line 155 | Line 156 | void SimSetup::createSim(void){
156    // initialize the Fortran
157  
158    initFortran();
159 +
160 +  if (globals->haveMinimizer())
161 +    // make minimizer
162 +    makeMinimizer();
163 +  else
164 +    // make the integrator
165 +    makeIntegrator();
166 +
167   }
168  
169  
170   void SimSetup::makeMolecules(void){
171 <  int k;
172 <  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
171 >  int i, j, k;
172 >  int exI, exJ, exK, exL, slI, slJ;
173 >  int tempI, tempJ, tempK, tempL;
174 >  int molI;
175 >  int stampID, atomOffset, rbOffset;
176    molInit molInfo;
177    DirectionalAtom* dAtom;
178 +  RigidBody* myRB;
179 +  StuntDouble* mySD;
180    LinkedAssign* extras;
181    LinkedAssign* current_extra;
182    AtomStamp* currentAtom;
183    BondStamp* currentBond;
184    BendStamp* currentBend;
185    TorsionStamp* currentTorsion;
186 +  RigidBodyStamp* currentRigidBody;
187 +  CutoffGroupStamp* currentCutoffGroup;
188 +  CutoffGroup* myCutoffGroup;
189 +  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
190 +  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
191  
192    bond_pair* theBonds;
193    bend_set* theBends;
194    torsion_set* theTorsions;
195  
196 +  set<int> skipList;
197  
198 +  double phi, theta, psi;
199 +  char* molName;
200 +  char rbName[100];
201 +
202 +  //ConstraintPair* consPair; //constraint pair
203 +  //ConstraintElement* consElement1;  //first element of constraint pair
204 +  //ConstraintElement* consElement2;  //second element of constraint pair
205 +  //int whichRigidBody;
206 +  //int consAtomIndex;  //index of constraint atom in rigid body's atom array
207 +  //vector<pair<int, int> > jointAtoms;
208    //init the forceField paramters
209  
210    the_ff->readParams();
211  
182 –
212    // init the atoms
213  
214 <  double ux, uy, uz, u, uSqr;
214 >  int nMembers, nNew, rb1, rb2;
215  
216    for (k = 0; k < nInfo; k++){
217      the_ff->setSimInfo(&(info[k]));
218  
219      atomOffset = 0;
220 <    excludeOffset = 0;
220 >    groupOffset = 0;
221 >
222      for (i = 0; i < info[k].n_mol; i++){
223        stampID = info[k].molecules[i].getStampID();
224 +      molName = comp_stamps[stampID]->getID();
225  
226        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
227        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
228        molInfo.nBends = comp_stamps[stampID]->getNBends();
229        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
230 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
230 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
231  
232 +      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
233 +      
234        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
202 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
203 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
204 –      molInfo.myBends = new Bend * [molInfo.nBends];
205 –      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
235  
236 +      if (molInfo.nBonds > 0)
237 +        molInfo.myBonds = new Bond*[molInfo.nBonds];
238 +      else
239 +        molInfo.myBonds = NULL;
240 +
241 +      if (molInfo.nBends > 0)
242 +        molInfo.myBends = new Bend*[molInfo.nBends];
243 +      else
244 +        molInfo.myBends = NULL;
245 +
246 +      if (molInfo.nTorsions > 0)
247 +        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
248 +      else
249 +        molInfo.myTorsions = NULL;
250 +
251        theBonds = new bond_pair[molInfo.nBonds];
252        theBends = new bend_set[molInfo.nBends];
253        theTorsions = new torsion_set[molInfo.nTorsions];
254 <
254 >      
255        // make the Atoms
256  
257        for (j = 0; j < molInfo.nAtoms; j++){
258          currentAtom = comp_stamps[stampID]->getAtom(j);
259 +
260          if (currentAtom->haveOrientation()){
261            dAtom = new DirectionalAtom((j + atomOffset),
262                                        info[k].getConfiguration());
263            info[k].n_oriented++;
264            molInfo.myAtoms[j] = dAtom;
265  
266 <          ux = currentAtom->getOrntX();
267 <          uy = currentAtom->getOrntY();
268 <          uz = currentAtom->getOrntZ();
266 >          // Directional Atoms have standard unit vectors which are oriented
267 >          // in space using the three Euler angles.  We assume the standard
268 >          // unit vector was originally along the z axis below.
269  
270 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
270 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
271 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
272 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
273  
274 <          u = sqrt(uSqr);
275 <          ux = ux / u;
229 <          uy = uy / u;
230 <          uz = uz / u;
231 <
232 <          dAtom->setSUx(ux);
233 <          dAtom->setSUy(uy);
234 <          dAtom->setSUz(uz);
274 >          dAtom->setUnitFrameFromEuler(phi, theta, psi);
275 >            
276          }
277          else{
278 <          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
279 <                                               info[k].getConfiguration());
278 >
279 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
280 >
281          }
240 –        molInfo.myAtoms[j]->setType(currentAtom->getType());
282  
283 +        molInfo.myAtoms[j]->setType(currentAtom->getType());
284   #ifdef IS_MPI
285  
286 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
286 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
287  
288   #endif // is_mpi
289        }
# Line 252 | Line 294 | void SimSetup::makeMolecules(void){
294          theBonds[j].a = currentBond->getA() + atomOffset;
295          theBonds[j].b = currentBond->getB() + atomOffset;
296  
297 <        exI = theBonds[j].a;
298 <        exJ = theBonds[j].b;
297 >        tempI = theBonds[j].a;
298 >        tempJ = theBonds[j].b;
299  
300 <        // exclude_I must always be the smaller of the pair
301 <        if (exI > exJ){
302 <          tempEx = exI;
303 <          exI = exJ;
304 <          exJ = tempEx;
305 <        }
306 < #ifdef IS_MPI
265 <        tempEx = exI;
266 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
267 <        tempEx = exJ;
268 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
300 > #ifdef IS_MPI
301 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
302 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
303 > #else
304 >        exI = tempI + 1;
305 >        exJ = tempJ + 1;
306 > #endif
307  
308 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
271 < #else  // isn't MPI
272 <
273 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
274 < #endif  //is_mpi
308 >        info[k].excludes->addPair(exI, exJ);
309        }
276 –      excludeOffset += molInfo.nBonds;
310  
311        //make the bends
312        for (j = 0; j < molInfo.nBends; j++){
# Line 323 | Line 356 | void SimSetup::makeMolecules(void){
356            }
357          }
358  
359 <        if (!theBends[j].isGhost){
360 <          exI = theBends[j].a;
361 <          exJ = theBends[j].c;
362 <        }
363 <        else{
331 <          exI = theBends[j].a;
332 <          exJ = theBends[j].b;
333 <        }
334 <
335 <        // exclude_I must always be the smaller of the pair
336 <        if (exI > exJ){
337 <          tempEx = exI;
338 <          exI = exJ;
339 <          exJ = tempEx;
340 <        }
359 >        if (theBends[j].isGhost) {
360 >          
361 >          tempI = theBends[j].a;
362 >          tempJ = theBends[j].b;
363 >          
364   #ifdef IS_MPI
365 <        tempEx = exI;
366 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
367 <        tempEx = exJ;
368 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
365 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
366 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
367 > #else
368 >          exI = tempI + 1;
369 >          exJ = tempJ + 1;
370 > #endif          
371 >          info[k].excludes->addPair(exI, exJ);
372  
373 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
374 < #else  // isn't MPI
375 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
376 < #endif  //is_mpi
373 >        } else {
374 >
375 >          tempI = theBends[j].a;
376 >          tempJ = theBends[j].b;
377 >          tempK = theBends[j].c;
378 >          
379 > #ifdef IS_MPI
380 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
381 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
382 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
383 > #else
384 >          exI = tempI + 1;
385 >          exJ = tempJ + 1;
386 >          exK = tempK + 1;
387 > #endif
388 >          
389 >          info[k].excludes->addPair(exI, exK);
390 >          info[k].excludes->addPair(exI, exJ);
391 >          info[k].excludes->addPair(exJ, exK);
392 >        }
393        }
352 –      excludeOffset += molInfo.nBends;
394  
395        for (j = 0; j < molInfo.nTorsions; j++){
396          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 358 | Line 399 | void SimSetup::makeMolecules(void){
399          theTorsions[j].c = currentTorsion->getC() + atomOffset;
400          theTorsions[j].d = currentTorsion->getD() + atomOffset;
401  
402 <        exI = theTorsions[j].a;
403 <        exJ = theTorsions[j].d;
402 >        tempI = theTorsions[j].a;      
403 >        tempJ = theTorsions[j].b;
404 >        tempK = theTorsions[j].c;
405 >        tempL = theTorsions[j].d;
406  
364 –        // exclude_I must always be the smaller of the pair
365 –        if (exI > exJ){
366 –          tempEx = exI;
367 –          exI = exJ;
368 –          exJ = tempEx;
369 –        }
407   #ifdef IS_MPI
408 <        tempEx = exI;
409 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
410 <        tempEx = exJ;
411 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
412 <
413 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
414 < #else  // isn't MPI
415 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
416 < #endif  //is_mpi
417 <      }
381 <      excludeOffset += molInfo.nTorsions;
408 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
409 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
410 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
411 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
412 > #else
413 >        exI = tempI + 1;
414 >        exJ = tempJ + 1;
415 >        exK = tempK + 1;
416 >        exL = tempL + 1;
417 > #endif
418  
419 +        info[k].excludes->addPair(exI, exJ);
420 +        info[k].excludes->addPair(exI, exK);
421 +        info[k].excludes->addPair(exI, exL);        
422 +        info[k].excludes->addPair(exJ, exK);
423 +        info[k].excludes->addPair(exJ, exL);
424 +        info[k].excludes->addPair(exK, exL);
425 +      }
426  
427 <      // send the arrays off to the forceField for init.
427 >      
428 >      molInfo.myRigidBodies.clear();
429 >      
430 >      for (j = 0; j < molInfo.nRigidBodies; j++){
431  
432 +        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
433 +        nMembers = currentRigidBody->getNMembers();
434 +
435 +        // Create the Rigid Body:
436 +
437 +        myRB = new RigidBody();
438 +
439 +        sprintf(rbName,"%s_RB_%d", molName, j);
440 +        myRB->setType(rbName);
441 +        
442 +        for (rb1 = 0; rb1 < nMembers; rb1++) {
443 +
444 +          // molI is atom numbering inside this molecule
445 +          molI = currentRigidBody->getMember(rb1);    
446 +
447 +          // tempI is atom numbering on local processor
448 +          tempI = molI + atomOffset;
449 +
450 +          // currentAtom is the AtomStamp (which we need for
451 +          // rigid body reference positions)
452 +          currentAtom = comp_stamps[stampID]->getAtom(molI);
453 +
454 +          // When we add to the rigid body, add the atom itself and
455 +          // the stamp info:
456 +
457 +          myRB->addAtom(info[k].atoms[tempI], currentAtom);
458 +          
459 +          // Add this atom to the Skip List for the integrators
460 + #ifdef IS_MPI
461 +          slI = info[k].atoms[tempI]->getGlobalIndex();
462 + #else
463 +          slI = tempI;
464 + #endif
465 +          skipList.insert(slI);
466 +          
467 +        }
468 +        
469 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
470 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
471 +            
472 +            tempI = currentRigidBody->getMember(rb1);
473 +            tempJ = currentRigidBody->getMember(rb2);
474 +            
475 +            // Some explanation is required here.
476 +            // Fortran indexing starts at 1, while c indexing starts at 0
477 +            // Also, in parallel computations, the GlobalIndex is
478 +            // used for the exclude list:
479 +            
480 + #ifdef IS_MPI
481 +            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
482 +            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
483 + #else
484 +            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
485 +            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
486 + #endif
487 +            
488 +            info[k].excludes->addPair(exI, exJ);
489 +            
490 +          }
491 +        }
492 +
493 +        molInfo.myRigidBodies.push_back(myRB);
494 +        info[k].rigidBodies.push_back(myRB);
495 +      }
496 +      
497 +
498 +      //create cutoff group for molecule
499 +
500 +      cutoffAtomSet.clear();
501 +      molInfo.myCutoffGroups.clear();
502 +      
503 +      for (j = 0; j < nCutoffGroups; j++){
504 +
505 +        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
506 +        nMembers = currentCutoffGroup->getNMembers();
507 +
508 +        myCutoffGroup = new CutoffGroup();
509 +        myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
510 +        
511 +        for (int cg = 0; cg < nMembers; cg++) {
512 +
513 +          // molI is atom numbering inside this molecule
514 +          molI = currentCutoffGroup->getMember(cg);    
515 +
516 +          // tempI is atom numbering on local processor
517 +          tempI = molI + atomOffset;
518 +
519 + #ifdef IS_MPI
520 +          globalID = info[k].atoms[tempI]->getGlobalIndex()
521 + #else
522 +          globalID = info[k].atoms[tempI]->getIndex();
523 + #endif
524 +
525 +          globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
526 +
527 +          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
528 +
529 +          cutoffAtomSet.insert(tempI);
530 +        }
531 +      
532 +        molInfo.myCutoffGroups.push_back(myCutoffGroup);
533 +        groupOffset++;
534 +
535 +      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
536 +
537 +      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
538 +
539 +      for(j = 0; j < molInfo.nAtoms; j++){
540 +
541 +        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
542 +          myCutoffGroup = new CutoffGroup();
543 +          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
544 +          myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
545 + #ifdef IS_MPI
546 +          globalID = info[k].atoms[atomOffset + j]->getGlobalIndex()
547 + #else
548 +          globalID = info[k].atoms[atomOffset + j]->getIndex();
549 + #endif
550 +          globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
551 +          molInfo.myCutoffGroups.push_back(myCutoffGroup);
552 +          groupOffset++;
553 +        }
554 +          
555 +      }
556 +
557 +      // After this is all set up, scan through the atoms to
558 +      // see if they can be added to the integrableObjects:
559 +
560 +      molInfo.myIntegrableObjects.clear();
561 +      
562 +
563 +      for (j = 0; j < molInfo.nAtoms; j++){
564 +
565 + #ifdef IS_MPI
566 +        slJ = molInfo.myAtoms[j]->getGlobalIndex();
567 + #else
568 +        slJ = j+atomOffset;
569 + #endif
570 +
571 +        // if they aren't on the skip list, then they can be integrated
572 +
573 +        if (skipList.find(slJ) == skipList.end()) {
574 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
575 +          info[k].integrableObjects.push_back(mySD);
576 +          molInfo.myIntegrableObjects.push_back(mySD);
577 +        }
578 +      }
579 +
580 +      // all rigid bodies are integrated:
581 +
582 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
583 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
584 +        info[k].integrableObjects.push_back(mySD);      
585 +        molInfo.myIntegrableObjects.push_back(mySD);
586 +      }
587 +
588 +
589 +    /*
590 +
591 +      //creat ConstraintPair.
592 +      molInfo.myConstraintPair.clear();
593 +      
594 +      for (j = 0; j < molInfo.nBonds; j++){
595 +
596 +        //if both atoms are in the same rigid body, just skip it
597 +        currentBond = comp_stamps[stampID]->getBond(j);
598 +        if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){
599 +
600 +          tempI = currentBond->getA() + atomOffset;
601 +          if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex))
602 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
603 +          else
604 +             consElement1 = new ConstraintAtom(info[k].atoms[tempI]);      
605 +
606 +          tempJ =  currentBond->getB() + atomOffset;
607 +          if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex))
608 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
609 +          else
610 +             consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);    
611 +
612 +          consPair = new DistanceConstraintPair(consElement1, consElement2);
613 +          molInfo.myConstraintPairs.push_back(consPair);
614 +        }
615 +      }  
616 +      
617 +      //loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair
618 +      for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){
619 +        for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){
620 +          
621 +          jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2);
622 +
623 +          for(size_t m = 0; m < jointAtoms.size(); m++){          
624 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first);
625 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second);
626 +
627 +            consPair = new JointConstraintPair(consElement1, consElement2);  
628 +            molInfo.myConstraintPairs.push_back(consPair);            
629 +          }
630 +
631 +        }
632 +      }
633 +      
634 + */      
635 +      // send the arrays off to the forceField for init.
636 +      
637        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
638        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
639        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
640        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
641                                   theTorsions);
642  
392 –
643        info[k].molecules[i].initialize(molInfo);
644  
645  
# Line 397 | Line 647 | void SimSetup::makeMolecules(void){
647        delete[] theBonds;
648        delete[] theBends;
649        delete[] theTorsions;
650 <    }
650 >    }    
651    }
652  
653   #ifdef IS_MPI
# Line 405 | Line 655 | void SimSetup::makeMolecules(void){
655    MPIcheckPoint();
656   #endif // is_mpi
657  
408 –  // clean up the forcefield
409 –
410 –  the_ff->calcRcut();
411 –  the_ff->cleanMe();
658   }
659  
660   void SimSetup::initFromBass(void){
# Line 609 | Line 855 | void SimSetup::gatherInfo(void){
855    else if (!strcasecmp(force_field, "EAM")){
856      ffCase = FF_EAM;
857    }
858 +  else if (!strcasecmp(force_field, "WATER")){
859 +    ffCase = FF_H2O;
860 +  }
861    else{
862      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
863              force_field);
# Line 692 | Line 941 | void SimSetup::gatherInfo(void){
941    }
942  
943    //check whether sample time, status time, thermal time and reset time are divisble by dt
944 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
944 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
945      sprintf(painCave.errMsg,
946              "Sample time is not divisible by dt.\n"
947              "\tThis will result in samples that are not uniformly\n"
# Line 702 | Line 951 | void SimSetup::gatherInfo(void){
951      simError();    
952    }
953  
954 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
954 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
955      sprintf(painCave.errMsg,
956              "Status time is not divisible by dt.\n"
957              "\tThis will result in status reports that are not uniformly\n"
# Line 738 | Line 987 | void SimSetup::gatherInfo(void){
987      if (globals->haveSampleTime()){
988        info[i].sampleTime = globals->getSampleTime();
989        info[i].statusTime = info[i].sampleTime;
741 –      info[i].thermalTime = info[i].sampleTime;
990      }
991      else{
992        info[i].sampleTime = globals->getRunTime();
993        info[i].statusTime = info[i].sampleTime;
746 –      info[i].thermalTime = info[i].sampleTime;
994      }
995  
996      if (globals->haveStatusTime()){
# Line 752 | Line 999 | void SimSetup::gatherInfo(void){
999  
1000      if (globals->haveThermalTime()){
1001        info[i].thermalTime = globals->getThermalTime();
1002 +    } else {
1003 +      info[i].thermalTime = globals->getRunTime();
1004      }
1005  
1006      info[i].resetIntegrator = 0;
# Line 769 | Line 1018 | void SimSetup::gatherInfo(void){
1018  
1019      info[i].useInitXSstate = globals->getUseInitXSstate();
1020      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
1021 <    
1021 >
1022 >    // check for thermodynamic integration
1023 >    if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) {
1024 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1025 >        info[i].useSolidThermInt = globals->getUseSolidThermInt();
1026 >        info[i].thermIntLambda = globals->getThermIntLambda();
1027 >        info[i].thermIntK = globals->getThermIntK();
1028 >        
1029 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
1030 >        info[i].restraint = myRestraint;
1031 >      }
1032 >      else {
1033 >        sprintf(painCave.errMsg,
1034 >                "SimSetup Error:\n"
1035 >                "\tKeyword useSolidThermInt was set to 'true' but\n"
1036 >                "\tthermodynamicIntegrationLambda (and/or\n"
1037 >                "\tthermodynamicIntegrationK) was not specified.\n"
1038 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1039 >        painCave.isFatal = 1;
1040 >        simError();    
1041 >      }
1042 >    }
1043 >    else if(globals->getUseLiquidThermInt()) {
1044 >      if (globals->getUseSolidThermInt()) {
1045 >        sprintf( painCave.errMsg,
1046 >                 "SimSetup Warning: It appears that you have both solid and\n"
1047 >                 "\tliquid thermodynamic integration activated in your .bass\n"
1048 >                 "\tfile. To avoid confusion, specify only one technique in\n"
1049 >                 "\tyour .bass file. Liquid-state thermodynamic integration\n"
1050 >                 "\twill be assumed for the current simulation. If this is not\n"
1051 >                 "\twhat you desire, set useSolidThermInt to 'true' and\n"
1052 >                 "\tuseLiquidThermInt to 'false' in your .bass file.\n");
1053 >        painCave.isFatal = 0;
1054 >        simError();
1055 >      }
1056 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1057 >        info[i].useLiquidThermInt = globals->getUseLiquidThermInt();
1058 >        info[i].thermIntLambda = globals->getThermIntLambda();
1059 >        info[i].thermIntK = globals->getThermIntK();
1060 >      }
1061 >      else {
1062 >        sprintf(painCave.errMsg,
1063 >                "SimSetup Error:\n"
1064 >                "\tKeyword useLiquidThermInt was set to 'true' but\n"
1065 >                "\tthermodynamicIntegrationLambda (and/or\n"
1066 >                "\tthermodynamicIntegrationK) was not specified.\n"
1067 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1068 >        painCave.isFatal = 1;
1069 >        simError();    
1070 >      }
1071 >    }
1072 >    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
1073 >        sprintf(painCave.errMsg,
1074 >                "SimSetup Warning: If you want to use Thermodynamic\n"
1075 >                "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n"
1076 >                "\t'true' in your .bass file.  These keywords are set to\n"
1077 >                "\t'false' by default, so your lambda and/or k values are\n"
1078 >                "\tbeing ignored.\n");
1079 >        painCave.isFatal = 0;
1080 >        simError();  
1081 >    }
1082    }
1083    
1084    //setup seed for random number generator
# Line 811 | Line 1120 | void SimSetup::gatherInfo(void){
1120    for (int i = 0; i < nInfo; i++){
1121      info[i].setSeed(seedValue);
1122    }
1123 <
1123 >  
1124   #ifdef IS_MPI
1125 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
1125 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
1126    MPIcheckPoint();
1127   #endif // is_mpi
1128   }
# Line 822 | Line 1131 | void SimSetup::finalInfoCheck(void){
1131   void SimSetup::finalInfoCheck(void){
1132    int index;
1133    int usesDipoles;
1134 +  int usesCharges;
1135    int i;
1136  
1137    for (i = 0; i < nInfo; i++){
# Line 833 | Line 1143 | void SimSetup::finalInfoCheck(void){
1143        usesDipoles = (info[i].atoms[index])->hasDipole();
1144        index++;
1145      }
1146 <
1146 >    index = 0;
1147 >    usesCharges = 0;
1148 >    while ((index < info[i].n_atoms) && !usesCharges){
1149 >      usesCharges= (info[i].atoms[index])->hasCharge();
1150 >      index++;
1151 >    }
1152   #ifdef IS_MPI
1153      int myUse = usesDipoles;
1154      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1155   #endif //is_mpi
1156  
1157 <    double theEcr, theEst;
1157 >    double theRcut, theRsw;
1158  
1159 +    if (globals->haveRcut()) {
1160 +      theRcut = globals->getRcut();
1161 +
1162 +      if (globals->haveRsw())
1163 +        theRsw = globals->getRsw();
1164 +      else
1165 +        theRsw = theRcut;
1166 +      
1167 +      info[i].setDefaultRcut(theRcut, theRsw);
1168 +
1169 +    } else {
1170 +      
1171 +      the_ff->calcRcut();
1172 +      theRcut = info[i].getRcut();
1173 +
1174 +      if (globals->haveRsw())
1175 +        theRsw = globals->getRsw();
1176 +      else
1177 +        theRsw = theRcut;
1178 +      
1179 +      info[i].setDefaultRcut(theRcut, theRsw);
1180 +    }
1181 +
1182      if (globals->getUseRF()){
1183        info[i].useReactionField = 1;
1184 <
1185 <      if (!globals->haveECR()){
1186 <        sprintf(painCave.errMsg,
1187 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1184 >      
1185 >      if (!globals->haveRcut()){
1186 >        sprintf(painCave.errMsg,
1187 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1188                  "\tOOPSE will use a default value of 15.0 angstroms"
1189 <                "\tfor the electrostaticCutoffRadius.\n");
1189 >                "\tfor the cutoffRadius.\n");
1190          painCave.isFatal = 0;
1191          simError();
1192 <        theEcr = 15.0;
1192 >        theRcut = 15.0;
1193        }
1194        else{
1195 <        theEcr = globals->getECR();
1195 >        theRcut = globals->getRcut();
1196        }
1197  
1198 <      if (!globals->haveEST()){
1198 >      if (!globals->haveRsw()){
1199          sprintf(painCave.errMsg,
1200 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1200 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1201                  "\tOOPSE will use a default value of\n"
1202 <                "\t0.05 * electrostaticCutoffRadius\n"
865 <                "\tfor the electrostaticSkinThickness\n");
1202 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1203          painCave.isFatal = 0;
1204          simError();
1205 <        theEst = 0.05 * theEcr;
1205 >        theRsw = 0.95 * theRcut;
1206        }
1207        else{
1208 <        theEst = globals->getEST();
1208 >        theRsw = globals->getRsw();
1209        }
1210  
1211 <      info[i].setDefaultEcr(theEcr, theEst);
1211 >      info[i].setDefaultRcut(theRcut, theRsw);
1212  
1213        if (!globals->haveDielectric()){
1214          sprintf(painCave.errMsg,
# Line 884 | Line 1221 | void SimSetup::finalInfoCheck(void){
1221        info[i].dielectric = globals->getDielectric();
1222      }
1223      else{
1224 <      if (usesDipoles){
1225 <        if (!globals->haveECR()){
1224 >      if (usesDipoles || usesCharges){
1225 >
1226 >        if (!globals->haveRcut()){
1227            sprintf(painCave.errMsg,
1228 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1228 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1229                    "\tOOPSE will use a default value of 15.0 angstroms"
1230 <                  "\tfor the electrostaticCutoffRadius.\n");
1231 <          painCave.isFatal = 0;
1232 <          simError();
1233 <          theEcr = 15.0;
1234 <        }
1230 >                  "\tfor the cutoffRadius.\n");
1231 >          painCave.isFatal = 0;
1232 >          simError();
1233 >          theRcut = 15.0;
1234 >      }
1235          else{
1236 <          theEcr = globals->getECR();
1236 >          theRcut = globals->getRcut();
1237          }
1238 <        
1239 <        if (!globals->haveEST()){
1238 >        
1239 >        if (!globals->haveRsw()){
1240            sprintf(painCave.errMsg,
1241 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1241 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1242                    "\tOOPSE will use a default value of\n"
1243 <                  "\t0.05 * electrostaticCutoffRadius\n"
906 <                  "\tfor the electrostaticSkinThickness\n");
1243 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1244            painCave.isFatal = 0;
1245            simError();
1246 <          theEst = 0.05 * theEcr;
1246 >          theRsw = 0.95 * theRcut;
1247          }
1248          else{
1249 <          theEst = globals->getEST();
1249 >          theRsw = globals->getRsw();
1250          }
1251 +        
1252 +        info[i].setDefaultRcut(theRcut, theRsw);
1253          
915 –        info[i].setDefaultEcr(theEcr, theEst);
1254        }
1255      }
1256    }
# Line 920 | Line 1258 | void SimSetup::finalInfoCheck(void){
1258    strcpy(checkPointMsg, "post processing checks out");
1259    MPIcheckPoint();
1260   #endif // is_mpi
1261 +
1262 +  // clean up the forcefield
1263 +  the_ff->cleanMe();
1264   }
1265    
1266   void SimSetup::initSystemCoords(void){
# Line 1050 | Line 1391 | void SimSetup::makeOutNames(void){
1391          }
1392        }
1393  
1394 +      strcpy(info[k].rawPotName, inFileName);
1395 +      nameLength = strlen(info[k].rawPotName);
1396 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1397 +      if (!strcmp(endTest, ".bass")){
1398 +        strcpy(endTest, ".raw");
1399 +      }
1400 +      else if (!strcmp(endTest, ".BASS")){
1401 +        strcpy(endTest, ".raw");
1402 +      }
1403 +      else{
1404 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1405 +        if (!strcmp(endTest, ".bss")){
1406 +          strcpy(endTest, ".raw");
1407 +        }
1408 +        else if (!strcmp(endTest, ".mdl")){
1409 +          strcpy(endTest, ".raw");
1410 +        }
1411 +        else{
1412 +          strcat(info[k].rawPotName, ".raw");
1413 +        }
1414 +      }
1415 +
1416   #ifdef IS_MPI
1417  
1418      }
# Line 1110 | Line 1473 | void SimSetup::createFF(void){
1473        the_ff = new EAM_FF();
1474        break;
1475  
1476 +    case FF_H2O:
1477 +      the_ff = new WATER();
1478 +      break;
1479 +
1480      default:
1481        sprintf(painCave.errMsg,
1482                "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1130 | Line 1497 | void SimSetup::compList(void){
1497    LinkedMolStamp* headStamp = new LinkedMolStamp();
1498    LinkedMolStamp* currentStamp = NULL;
1499    comp_stamps = new MoleculeStamp * [n_components];
1500 +  bool haveCutoffGroups;
1501  
1502 +  haveCutoffGroups = false;
1503 +  
1504    // make an array of molecule stamps that match the components used.
1505    // also extract the used stamps out into a separate linked list
1506  
# Line 1165 | Line 1535 | void SimSetup::compList(void){
1535        headStamp->add(currentStamp);
1536        comp_stamps[i] = headStamp->match(id);
1537      }
1538 +
1539 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1540 +      haveCutoffGroups = true;    
1541    }
1542 +    
1543 +  for (i = 0; i < nInfo; i++)
1544 +    info[i].haveCutoffGroups = haveCutoffGroups;
1545  
1546   #ifdef IS_MPI
1547    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1174 | Line 1550 | void SimSetup::calcSysValues(void){
1550   }
1551  
1552   void SimSetup::calcSysValues(void){
1553 <  int i;
1553 >  int i, j;
1554 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1555  
1556    int* molMembershipArray;
1557 +  CutoffGroupStamp* cg;
1558  
1559    tot_atoms = 0;
1560    tot_bonds = 0;
1561    tot_bends = 0;
1562    tot_torsions = 0;
1563 +  tot_rigid = 0;
1564 +  tot_groups = 0;
1565    for (i = 0; i < n_components; i++){
1566      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1567      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1568      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1569      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1570 <  }
1570 >    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1571  
1572 +    ncutgroups = comp_stamps[i]->getNCutoffGroups();
1573 +    atomsingroups = 0;
1574 +    for (j=0; j < ncutgroups; j++) {
1575 +      cg = comp_stamps[i]->getCutoffGroup(j);
1576 +      atomsingroups += cg->getNMembers();
1577 +    }
1578 +    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
1579 +    tot_groups += components_nmol[i] * ngroupsinstamp;    
1580 +  }
1581 +  
1582    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1583    molMembershipArray = new int[tot_atoms];
1584  
# Line 1199 | Line 1589 | void SimSetup::calcSysValues(void){
1589      info[i].n_torsions = tot_torsions;
1590      info[i].n_SRI = tot_SRI;
1591      info[i].n_mol = tot_nmol;
1592 <
1592 >    info[i].ngroup = tot_groups;
1593      info[i].molMembershipArray = molMembershipArray;
1594    }
1595   }
# Line 1210 | Line 1600 | void SimSetup::mpiMolDivide(void){
1600    int i, j, k;
1601    int localMol, allMol;
1602    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1603 +  int local_rigid, local_groups;
1604 +  vector<int> globalMolIndex;
1605 +  int ncutgroups, atomsingroups, ngroupsinstamp;
1606 +  CutoffGroupStamp* cg;
1607  
1608    mpiSim = new mpiSimulation(info);
1609  
1610 <  globalIndex = mpiSim->divideLabor();
1610 >  mpiSim->divideLabor();
1611 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1612 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1613  
1614    // set up the local variables
1615  
# Line 1226 | Line 1622 | void SimSetup::mpiMolDivide(void){
1622    local_bonds = 0;
1623    local_bends = 0;
1624    local_torsions = 0;
1625 <  globalAtomIndex = 0;
1625 >  local_rigid = 0;
1626 >  local_groups = 0;
1627 >  globalAtomCounter = 0;
1628  
1231 –
1629    for (i = 0; i < n_components; i++){
1630      for (j = 0; j < components_nmol[i]; j++){
1631        if (mol2proc[allMol] == worldRank){
# Line 1236 | Line 1633 | void SimSetup::mpiMolDivide(void){
1633          local_bonds += comp_stamps[i]->getNBonds();
1634          local_bends += comp_stamps[i]->getNBends();
1635          local_torsions += comp_stamps[i]->getNTorsions();
1636 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1637 +
1638 +        ncutgroups = comp_stamps[i]->getNCutoffGroups();
1639 +        atomsingroups = 0;
1640 +        for (k=0; k < ncutgroups; k++) {
1641 +          cg = comp_stamps[i]->getCutoffGroup(k);
1642 +          atomsingroups += cg->getNMembers();
1643 +        }
1644 +        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +
1645 +          ncutgroups;
1646 +        local_groups += ngroupsinstamp;    
1647 +
1648          localMol++;
1649        }      
1650        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1651 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1652 <        globalAtomIndex++;
1651 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1652 >        globalAtomCounter++;
1653        }
1654  
1655        allMol++;
# Line 1248 | Line 1657 | void SimSetup::mpiMolDivide(void){
1657    }
1658    local_SRI = local_bonds + local_bends + local_torsions;
1659  
1660 <  info[0].n_atoms = mpiSim->getMyNlocal();  
1661 <
1660 >  info[0].n_atoms = mpiSim->getNAtomsLocal();  
1661 >  
1662    if (local_atoms != info[0].n_atoms){
1663      sprintf(painCave.errMsg,
1664              "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
# Line 1259 | Line 1668 | void SimSetup::mpiMolDivide(void){
1668      simError();
1669    }
1670  
1671 +  info[0].ngroup = mpiSim->getNGroupsLocal();  
1672 +  if (local_groups != info[0].ngroup){
1673 +    sprintf(painCave.errMsg,
1674 +            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
1675 +            "\tlocalGroups (%d) are not equal.\n",
1676 +            info[0].ngroup, local_groups);
1677 +    painCave.isFatal = 1;
1678 +    simError();
1679 +  }
1680 +  
1681    info[0].n_bonds = local_bonds;
1682    info[0].n_bends = local_bends;
1683    info[0].n_torsions = local_torsions;
# Line 1281 | Line 1700 | void SimSetup::makeSysArrays(void){
1700  
1701    Atom** the_atoms;
1702    Molecule* the_molecules;
1284 –  Exclude** the_excludes;
1703  
1286 –
1704    for (l = 0; l < nInfo; l++){
1705      // create the atom and short range interaction arrays
1706  
# Line 1297 | Line 1714 | void SimSetup::makeSysArrays(void){
1714  
1715  
1716      molIndex = 0;
1717 <    for (i = 0; i < mpiSim->getTotNmol(); i++){
1717 >    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
1718        if (mol2proc[i] == worldRank){
1719          the_molecules[molIndex].setStampID(molCompType[i]);
1720          the_molecules[molIndex].setMyIndex(molIndex);
# Line 1309 | Line 1726 | void SimSetup::makeSysArrays(void){
1726   #else // is_mpi
1727  
1728      molIndex = 0;
1729 <    globalAtomIndex = 0;
1729 >    globalAtomCounter = 0;
1730      for (i = 0; i < n_components; i++){
1731        for (j = 0; j < components_nmol[i]; j++){
1732          the_molecules[molIndex].setStampID(i);
1733          the_molecules[molIndex].setMyIndex(molIndex);
1734          the_molecules[molIndex].setGlobalIndex(molIndex);
1735          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1736 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1737 <          globalAtomIndex++;
1736 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1737 >          globalAtomCounter++;
1738          }
1739          molIndex++;
1740        }
# Line 1326 | Line 1743 | void SimSetup::makeSysArrays(void){
1743  
1744   #endif // is_mpi
1745  
1746 <
1747 <    if (info[l].n_SRI){
1748 <      Exclude::createArray(info[l].n_SRI);
1332 <      the_excludes = new Exclude * [info[l].n_SRI];
1333 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1334 <        the_excludes[ex] = new Exclude(ex);
1335 <      }
1336 <      info[l].globalExcludes = new int;
1337 <      info[l].n_exclude = info[l].n_SRI;
1338 <    }
1339 <    else{
1340 <      Exclude::createArray(1);
1341 <      the_excludes = new Exclude * ;
1342 <      the_excludes[0] = new Exclude(0);
1343 <      the_excludes[0]->setPair(0, 0);
1344 <      info[l].globalExcludes = new int;
1345 <      info[l].globalExcludes[0] = 0;
1346 <      info[l].n_exclude = 0;
1347 <    }
1348 <
1746 >    info[l].globalExcludes = new int;
1747 >    info[l].globalExcludes[0] = 0;
1748 >    
1749      // set the arrays into the SimInfo object
1750  
1751      info[l].atoms = the_atoms;
1752      info[l].molecules = the_molecules;
1753      info[l].nGlobalExcludes = 0;
1754 <    info[l].excludes = the_excludes;
1355 <
1754 >    
1755      the_ff->setSimInfo(info);
1756    }
1757   }
# Line 1624 | Line 2023 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2023    }
2024  
2025    theInfo.addProperty(zconsForcePolicy);
2026 +
2027 +  //set zcons gap
2028 +  DoubleData* zconsGap = new DoubleData();
2029 +  zconsGap->setID(ZCONSGAP_ID);
2030 +
2031 +  if (globals->haveZConsGap()){
2032 +    zconsGap->setData(globals->getZconsGap());
2033 +    theInfo.addProperty(zconsGap);  
2034 +  }
2035  
2036 +  //set zcons fixtime
2037 +  DoubleData* zconsFixtime = new DoubleData();
2038 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
2039 +
2040 +  if (globals->haveZConsFixTime()){
2041 +    zconsFixtime->setData(globals->getZconsFixtime());
2042 +    theInfo.addProperty(zconsFixtime);  
2043 +  }
2044 +
2045 +  //set zconsUsingSMD
2046 +  IntData* zconsUsingSMD = new IntData();
2047 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
2048 +
2049 +  if (globals->haveZConsUsingSMD()){
2050 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
2051 +    theInfo.addProperty(zconsUsingSMD);  
2052 +  }
2053 +
2054    //Determine the name of ouput file and add it into SimInfo's property list
2055    //Be careful, do not use inFileName, since it is a pointer which
2056    //point to a string at master node, and slave nodes do not contain that string
# Line 1654 | Line 2080 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2080      tempParaItem.zPos = zconStamp[i]->getZpos();
2081      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
2082      tempParaItem.kRatio = zconStamp[i]->getKratio();
2083 <
2083 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
2084 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
2085      zconsParaData->addItem(tempParaItem);
2086    }
2087  
# Line 1671 | Line 2098 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2098    
2099    //push data into siminfo, therefore, we can retrieve later
2100    theInfo.addProperty(zconsParaData);
2101 + }
2102 +
2103 + void SimSetup::makeMinimizer(){
2104 +
2105 +  OOPSEMinimizer* myOOPSEMinimizer;
2106 +  MinimizerParameterSet* param;
2107 +  char minimizerName[100];
2108 +  
2109 +  for (int i = 0; i < nInfo; i++){
2110 +    
2111 +    //prepare parameter set for minimizer
2112 +    param = new MinimizerParameterSet();
2113 +    param->setDefaultParameter();
2114 +
2115 +    if (globals->haveMinimizer()){
2116 +      param->setFTol(globals->getMinFTol());
2117 +    }
2118 +
2119 +    if (globals->haveMinGTol()){
2120 +      param->setGTol(globals->getMinGTol());
2121 +    }
2122 +
2123 +    if (globals->haveMinMaxIter()){
2124 +      param->setMaxIteration(globals->getMinMaxIter());
2125 +    }
2126 +
2127 +    if (globals->haveMinWriteFrq()){
2128 +      param->setMaxIteration(globals->getMinMaxIter());
2129 +    }
2130 +
2131 +    if (globals->haveMinWriteFrq()){
2132 +      param->setWriteFrq(globals->getMinWriteFrq());
2133 +    }
2134 +    
2135 +    if (globals->haveMinStepSize()){
2136 +      param->setStepSize(globals->getMinStepSize());
2137 +    }
2138 +
2139 +    if (globals->haveMinLSMaxIter()){
2140 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
2141 +    }    
2142 +
2143 +    if (globals->haveMinLSTol()){
2144 +      param->setLineSearchTol(globals->getMinLSTol());
2145 +    }    
2146 +
2147 +    strcpy(minimizerName, globals->getMinimizer());
2148 +
2149 +    if (!strcasecmp(minimizerName, "CG")){
2150 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
2151 +    }
2152 +    else if (!strcasecmp(minimizerName, "SD")){
2153 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
2154 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
2155 +    }
2156 +    else{
2157 +          sprintf(painCave.errMsg,
2158 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
2159 +          painCave.isFatal = 0;
2160 +          simError();
2161 +
2162 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
2163 +    }
2164 +     info[i].the_integrator = myOOPSEMinimizer;
2165 +
2166 +     //store the minimizer into simInfo
2167 +     info[i].the_minimizer = myOOPSEMinimizer;
2168 +     info[i].has_minimizer = true;
2169 +  }
2170 +
2171   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines