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 999 by chrisfen, Fri Jan 30 15:01:09 2004 UTC vs.
Revision 1261 by gezelter, Fri Jun 11 14:14:10 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 + #include "ConstraintManager.hpp"
17  
18   #ifdef IS_MPI
19   #include "mpiBASS.h"
# Line 27 | Line 32
32   #define FF_DUFF  0
33   #define FF_LJ    1
34   #define FF_EAM   2
35 < #define FF_H2O 3
35 > #define FF_H2O   3
36  
37   using namespace std;
38  
# Line 144 | Line 149 | void SimSetup::createSim(void){
149    // make the output filenames
150  
151    makeOutNames();
152 <
148 <  // make the integrator
149 <
150 <  makeIntegrator();
151 <
152 >  
153   #ifdef IS_MPI
154    mpiSim->mpiRefresh();
155   #endif
# Line 156 | 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;
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  
182 –
218    // init the atoms
219  
220 <  double phi, theta, psi;
186 <  double sux, suy, suz;
187 <  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
188 <  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]);
205 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
206 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
207 –      molInfo.myBends = new Bend * [molInfo.nBends];
208 –      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());
# Line 228 | Line 286 | void SimSetup::makeMolecules(void){
286            phi = currentAtom->getEulerPhi() * M_PI / 180.0;
287            theta = currentAtom->getEulerTheta() * M_PI / 180.0;
288            psi = currentAtom->getEulerPsi()* M_PI / 180.0;
231 –            
232 –          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
233 –          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
234 –          Axz = sin(theta) * sin(psi);
235 –          
236 –          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
237 –          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
238 –          Ayz = sin(theta) * cos(psi);
239 –          
240 –          Azx = sin(phi) * sin(theta);
241 –          Azy = -cos(phi) * sin(theta);
242 –          Azz = cos(theta);
289  
290 <          sux = 0.0;
291 <          suy = 0.0;
246 <          suz = 1.0;
247 <
248 <          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
249 <          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
250 <          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
251 <
252 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
253 <
254 <          u = sqrt(uSqr);
255 <          ux = ux / u;
256 <          uy = uy / u;
257 <          uz = uz / u;
258 <
259 <          dAtom->setSUx(ux);
260 <          dAtom->setSUy(uy);
261 <          dAtom->setSUz(uz);
290 >          dAtom->setUnitFrameFromEuler(phi, theta, psi);
291 >            
292          }
293          else{
264 –          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
265 –                                               info[k].getConfiguration());
266 –        }
267 –        molInfo.myAtoms[j]->setType(currentAtom->getType());
294  
295 < #ifdef IS_MPI
295 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
296  
297 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
297 >        }
298  
299 +        molInfo.myAtoms[j]->setType(currentAtom->getType());
300 + #ifdef IS_MPI
301 +        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
302   #endif // is_mpi
303        }
304  
# Line 279 | 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  
285 –        // exclude_I must always be the smaller of the pair
286 –        if (exI > exJ){
287 –          tempEx = exI;
288 –          exI = exJ;
289 –          exJ = tempEx;
290 –        }
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);
298 < #else  // isn't MPI
299 <
300 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
301 < #endif  //is_mpi
322 >        info[k].excludes->addPair(exI, exJ);
323        }
303 –      excludeOffset += molInfo.nBonds;
324  
325        //make the bends
326        for (j = 0; j < molInfo.nBends; j++){
# Line 350 | 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{
358 <          exI = theBends[j].a;
359 <          exJ = theBends[j].b;
360 <        }
361 <
362 <        // exclude_I must always be the smaller of the pair
363 <        if (exI > exJ){
364 <          tempEx = exI;
365 <          exI = exJ;
366 <          exJ = tempEx;
367 <        }
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        }
379 –      excludeOffset += molInfo.nBends;
408  
409        for (j = 0; j < molInfo.nTorsions; j++){
410          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 385 | 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;
418 <
419 <        // exclude_I must always be the smaller of the pair
392 <        if (exI > exJ){
393 <          tempEx = exI;
394 <          exI = exJ;
395 <          exJ = tempEx;
396 <        }
397 < #ifdef IS_MPI
398 <        tempEx = exI;
399 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
400 <        tempEx = exJ;
401 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
402 <
403 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
404 < #else  // isn't MPI
405 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
406 < #endif  //is_mpi
407 <      }
408 <      excludeOffset += molInfo.nTorsions;
416 >        tempI = theTorsions[j].a;      
417 >        tempJ = theTorsions[j].b;
418 >        tempK = theTorsions[j].c;
419 >        tempL = theTorsions[j].d;
420  
421 + #ifdef IS_MPI
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 <      // send the arrays off to the forceField for init.
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 >      }
440  
441 +      
442 +      molInfo.myRigidBodies.clear();
443 +      
444 +      for (j = 0; j < molInfo.nRigidBodies; j++){
445 +
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 417 | 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 432 | Line 707 | void SimSetup::makeMolecules(void){
707    MPIcheckPoint();
708   #endif // is_mpi
709  
435 –  // clean up the forcefield
436 –
437 –  the_ff->calcRcut();
438 –  the_ff->cleanMe();
710   }
711  
712   void SimSetup::initFromBass(void){
# Line 645 | Line 916 | void SimSetup::gatherInfo(void){
916           painCave.isFatal = 1;
917           simError();
918    }
919 <
920 <    // get the ensemble
919 >  if (globals->haveForceFieldVariant()) {
920 >    strcpy(forcefield_variant, globals->getForceFieldVariant());
921 >    has_forcefield_variant = 1;
922 >  }
923 >  
924 >  // get the ensemble
925  
926    strcpy(ensemble, globals->getEnsemble());
927  
# Line 722 | Line 997 | void SimSetup::gatherInfo(void){
997    }
998  
999    //check whether sample time, status time, thermal time and reset time are divisble by dt
1000 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
1000 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
1001      sprintf(painCave.errMsg,
1002              "Sample time is not divisible by dt.\n"
1003              "\tThis will result in samples that are not uniformly\n"
# Line 732 | Line 1007 | void SimSetup::gatherInfo(void){
1007      simError();    
1008    }
1009  
1010 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
1010 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
1011      sprintf(painCave.errMsg,
1012              "Status time is not divisible by dt.\n"
1013              "\tThis will result in status reports that are not uniformly\n"
# Line 768 | Line 1043 | void SimSetup::gatherInfo(void){
1043      if (globals->haveSampleTime()){
1044        info[i].sampleTime = globals->getSampleTime();
1045        info[i].statusTime = info[i].sampleTime;
771 –      info[i].thermalTime = info[i].sampleTime;
1046      }
1047      else{
1048        info[i].sampleTime = globals->getRunTime();
1049        info[i].statusTime = info[i].sampleTime;
776 –      info[i].thermalTime = info[i].sampleTime;
1050      }
1051  
1052      if (globals->haveStatusTime()){
# Line 782 | Line 1055 | void SimSetup::gatherInfo(void){
1055  
1056      if (globals->haveThermalTime()){
1057        info[i].thermalTime = globals->getThermalTime();
1058 +    } else {
1059 +      info[i].thermalTime = globals->getRunTime();
1060      }
1061  
1062      info[i].resetIntegrator = 0;
# Line 799 | Line 1074 | void SimSetup::gatherInfo(void){
1074  
1075      info[i].useInitXSstate = globals->getUseInitXSstate();
1076      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
1077 <    
1077 >
1078 >    // check for thermodynamic integration
1079 >    if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) {
1080 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1081 >        info[i].useSolidThermInt = globals->getUseSolidThermInt();
1082 >        info[i].thermIntLambda = globals->getThermIntLambda();
1083 >        info[i].thermIntK = globals->getThermIntK();
1084 >        
1085 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
1086 >        info[i].restraint = myRestraint;
1087 >      }
1088 >      else {
1089 >        sprintf(painCave.errMsg,
1090 >                "SimSetup Error:\n"
1091 >                "\tKeyword useSolidThermInt was set to 'true' but\n"
1092 >                "\tthermodynamicIntegrationLambda (and/or\n"
1093 >                "\tthermodynamicIntegrationK) was not specified.\n"
1094 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1095 >        painCave.isFatal = 1;
1096 >        simError();    
1097 >      }
1098 >    }
1099 >    else if(globals->getUseLiquidThermInt()) {
1100 >      if (globals->getUseSolidThermInt()) {
1101 >        sprintf( painCave.errMsg,
1102 >                 "SimSetup Warning: It appears that you have both solid and\n"
1103 >                 "\tliquid thermodynamic integration activated in your .bass\n"
1104 >                 "\tfile. To avoid confusion, specify only one technique in\n"
1105 >                 "\tyour .bass file. Liquid-state thermodynamic integration\n"
1106 >                 "\twill be assumed for the current simulation. If this is not\n"
1107 >                 "\twhat you desire, set useSolidThermInt to 'true' and\n"
1108 >                 "\tuseLiquidThermInt to 'false' in your .bass file.\n");
1109 >        painCave.isFatal = 0;
1110 >        simError();
1111 >      }
1112 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1113 >        info[i].useLiquidThermInt = globals->getUseLiquidThermInt();
1114 >        info[i].thermIntLambda = globals->getThermIntLambda();
1115 >        info[i].thermIntK = globals->getThermIntK();
1116 >      }
1117 >      else {
1118 >        sprintf(painCave.errMsg,
1119 >                "SimSetup Error:\n"
1120 >                "\tKeyword useLiquidThermInt was set to 'true' but\n"
1121 >                "\tthermodynamicIntegrationLambda (and/or\n"
1122 >                "\tthermodynamicIntegrationK) was not specified.\n"
1123 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1124 >        painCave.isFatal = 1;
1125 >        simError();    
1126 >      }
1127 >    }
1128 >    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
1129 >        sprintf(painCave.errMsg,
1130 >                "SimSetup Warning: If you want to use Thermodynamic\n"
1131 >                "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n"
1132 >                "\t'true' in your .bass file.  These keywords are set to\n"
1133 >                "\t'false' by default, so your lambda and/or k values are\n"
1134 >                "\tbeing ignored.\n");
1135 >        painCave.isFatal = 0;
1136 >        simError();  
1137 >    }
1138    }
1139    
1140    //setup seed for random number generator
# Line 841 | Line 1176 | void SimSetup::gatherInfo(void){
1176    for (int i = 0; i < nInfo; i++){
1177      info[i].setSeed(seedValue);
1178    }
1179 <
1179 >  
1180   #ifdef IS_MPI
1181    strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
1182    MPIcheckPoint();
# Line 852 | Line 1187 | void SimSetup::finalInfoCheck(void){
1187   void SimSetup::finalInfoCheck(void){
1188    int index;
1189    int usesDipoles;
1190 +  int usesCharges;
1191    int i;
1192  
1193    for (i = 0; i < nInfo; i++){
# Line 863 | Line 1199 | void SimSetup::finalInfoCheck(void){
1199        usesDipoles = (info[i].atoms[index])->hasDipole();
1200        index++;
1201      }
1202 <
1202 >    index = 0;
1203 >    usesCharges = 0;
1204 >    while ((index < info[i].n_atoms) && !usesCharges){
1205 >      usesCharges= (info[i].atoms[index])->hasCharge();
1206 >      index++;
1207 >    }
1208   #ifdef IS_MPI
1209      int myUse = usesDipoles;
1210      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1211   #endif //is_mpi
1212  
1213 <    double theEcr, theEst;
1213 >    double theRcut, theRsw;
1214 >
1215 >    if (globals->haveRcut()) {
1216 >      theRcut = globals->getRcut();
1217 >
1218 >      if (globals->haveRsw())
1219 >        theRsw = globals->getRsw();
1220 >      else
1221 >        theRsw = theRcut;
1222 >      
1223 >      info[i].setDefaultRcut(theRcut, theRsw);
1224  
1225 +    } else {
1226 +      
1227 +      the_ff->calcRcut();
1228 +      theRcut = info[i].getRcut();
1229 +
1230 +      if (globals->haveRsw())
1231 +        theRsw = globals->getRsw();
1232 +      else
1233 +        theRsw = theRcut;
1234 +      
1235 +      info[i].setDefaultRcut(theRcut, theRsw);
1236 +    }
1237 +
1238      if (globals->getUseRF()){
1239        info[i].useReactionField = 1;
1240 <
1241 <      if (!globals->haveECR()){
1240 >      
1241 >      if (!globals->haveRcut()){
1242          sprintf(painCave.errMsg,
1243 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1243 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1244                  "\tOOPSE will use a default value of 15.0 angstroms"
1245 <                "\tfor the electrostaticCutoffRadius.\n");
1245 >                "\tfor the cutoffRadius.\n");
1246          painCave.isFatal = 0;
1247          simError();
1248 <        theEcr = 15.0;
1248 >        theRcut = 15.0;
1249        }
1250        else{
1251 <        theEcr = globals->getECR();
1251 >        theRcut = globals->getRcut();
1252        }
1253  
1254 <      if (!globals->haveEST()){
1254 >      if (!globals->haveRsw()){
1255          sprintf(painCave.errMsg,
1256 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1256 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1257                  "\tOOPSE will use a default value of\n"
1258 <                "\t0.05 * electrostaticCutoffRadius\n"
895 <                "\tfor the electrostaticSkinThickness\n");
1258 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1259          painCave.isFatal = 0;
1260          simError();
1261 <        theEst = 0.05 * theEcr;
1261 >        theRsw = 0.95 * theRcut;
1262        }
1263        else{
1264 <        theEst = globals->getEST();
1264 >        theRsw = globals->getRsw();
1265        }
1266  
1267 <      info[i].setDefaultEcr(theEcr, theEst);
1267 >      info[i].setDefaultRcut(theRcut, theRsw);
1268  
1269        if (!globals->haveDielectric()){
1270          sprintf(painCave.errMsg,
# Line 914 | Line 1277 | void SimSetup::finalInfoCheck(void){
1277        info[i].dielectric = globals->getDielectric();
1278      }
1279      else{
1280 <      if (usesDipoles){
1281 <        if (!globals->haveECR()){
1280 >      if (usesDipoles || usesCharges){
1281 >
1282 >        if (!globals->haveRcut()){
1283            sprintf(painCave.errMsg,
1284 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1284 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1285                    "\tOOPSE will use a default value of 15.0 angstroms"
1286 <                  "\tfor the electrostaticCutoffRadius.\n");
1287 <          painCave.isFatal = 0;
1288 <          simError();
1289 <          theEcr = 15.0;
1290 <        }
1286 >                  "\tfor the cutoffRadius.\n");
1287 >          painCave.isFatal = 0;
1288 >          simError();
1289 >          theRcut = 15.0;
1290 >      }
1291          else{
1292 <          theEcr = globals->getECR();
1292 >          theRcut = globals->getRcut();
1293          }
1294 <        
1295 <        if (!globals->haveEST()){
1294 >        
1295 >        if (!globals->haveRsw()){
1296            sprintf(painCave.errMsg,
1297 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1297 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1298                    "\tOOPSE will use a default value of\n"
1299 <                  "\t0.05 * electrostaticCutoffRadius\n"
936 <                  "\tfor the electrostaticSkinThickness\n");
1299 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1300            painCave.isFatal = 0;
1301            simError();
1302 <          theEst = 0.05 * theEcr;
1302 >          theRsw = 0.95 * theRcut;
1303          }
1304          else{
1305 <          theEst = globals->getEST();
1305 >          theRsw = globals->getRsw();
1306          }
1307 +        
1308 +        info[i].setDefaultRcut(theRcut, theRsw);
1309          
945 –        info[i].setDefaultEcr(theEcr, theEst);
1310        }
1311      }
1312    }
# Line 950 | Line 1314 | void SimSetup::finalInfoCheck(void){
1314    strcpy(checkPointMsg, "post processing checks out");
1315    MPIcheckPoint();
1316   #endif // is_mpi
1317 +
1318 +  // clean up the forcefield
1319 +  the_ff->cleanMe();
1320   }
1321    
1322   void SimSetup::initSystemCoords(void){
# Line 1077 | Line 1444 | void SimSetup::makeOutNames(void){
1444          }
1445          else{
1446            strcat(info[k].statusName, ".stat");
1447 +        }
1448 +      }
1449 +
1450 +      strcpy(info[k].rawPotName, inFileName);
1451 +      nameLength = strlen(info[k].rawPotName);
1452 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1453 +      if (!strcmp(endTest, ".bass")){
1454 +        strcpy(endTest, ".raw");
1455 +      }
1456 +      else if (!strcmp(endTest, ".BASS")){
1457 +        strcpy(endTest, ".raw");
1458 +      }
1459 +      else{
1460 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1461 +        if (!strcmp(endTest, ".bss")){
1462 +          strcpy(endTest, ".raw");
1463          }
1464 +        else if (!strcmp(endTest, ".mdl")){
1465 +          strcpy(endTest, ".raw");
1466 +        }
1467 +        else{
1468 +          strcat(info[k].rawPotName, ".raw");
1469 +        }
1470        }
1471  
1472   #ifdef IS_MPI
# Line 1129 | Line 1518 | void SimSetup::createFF(void){
1518   void SimSetup::createFF(void){
1519    switch (ffCase){
1520      case FF_DUFF:
1521 <      the_ff = new DUFF();
1521 >        the_ff = new DUFF();
1522        break;
1523  
1524      case FF_LJ:
# Line 1137 | Line 1526 | void SimSetup::createFF(void){
1526        break;
1527  
1528      case FF_EAM:
1529 <      the_ff = new EAM_FF();
1529 >      if (has_forcefield_variant)
1530 >        the_ff = new EAM_FF(forcefield_variant);
1531 >      else
1532 >        the_ff = new EAM_FF();
1533        break;
1534  
1535      case FF_H2O:
# Line 1151 | Line 1543 | void SimSetup::createFF(void){
1543        simError();
1544    }
1545  
1546 +
1547   #ifdef IS_MPI
1548    strcpy(checkPointMsg, "ForceField creation successful");
1549    MPIcheckPoint();
# Line 1164 | Line 1557 | void SimSetup::compList(void){
1557    LinkedMolStamp* headStamp = new LinkedMolStamp();
1558    LinkedMolStamp* currentStamp = NULL;
1559    comp_stamps = new MoleculeStamp * [n_components];
1560 +  bool haveCutoffGroups;
1561  
1562 +  haveCutoffGroups = false;
1563 +  
1564    // make an array of molecule stamps that match the components used.
1565    // also extract the used stamps out into a separate linked list
1566  
# Line 1199 | Line 1595 | void SimSetup::compList(void){
1595        headStamp->add(currentStamp);
1596        comp_stamps[i] = headStamp->match(id);
1597      }
1598 +
1599 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1600 +      haveCutoffGroups = true;    
1601    }
1602 +    
1603 +  for (i = 0; i < nInfo; i++)
1604 +    info[i].haveCutoffGroups = haveCutoffGroups;
1605  
1606   #ifdef IS_MPI
1607    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1208 | Line 1610 | void SimSetup::calcSysValues(void){
1610   }
1611  
1612   void SimSetup::calcSysValues(void){
1613 <  int i;
1613 >  int i, j;
1614 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1615  
1616    int* molMembershipArray;
1617 +  CutoffGroupStamp* cg;
1618  
1619    tot_atoms = 0;
1620    tot_bonds = 0;
1621    tot_bends = 0;
1622    tot_torsions = 0;
1623 +  tot_rigid = 0;
1624 +  tot_groups = 0;
1625    for (i = 0; i < n_components; i++){
1626      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1627      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1628      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1629      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1630 <  }
1630 >    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1631  
1632 +    ncutgroups = comp_stamps[i]->getNCutoffGroups();
1633 +    atomsingroups = 0;
1634 +    for (j=0; j < ncutgroups; j++) {
1635 +      cg = comp_stamps[i]->getCutoffGroup(j);
1636 +      atomsingroups += cg->getNMembers();
1637 +    }
1638 +    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
1639 +    tot_groups += components_nmol[i] * ngroupsinstamp;    
1640 +  }
1641 +  
1642    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1643    molMembershipArray = new int[tot_atoms];
1644  
# Line 1233 | Line 1649 | void SimSetup::calcSysValues(void){
1649      info[i].n_torsions = tot_torsions;
1650      info[i].n_SRI = tot_SRI;
1651      info[i].n_mol = tot_nmol;
1652 <
1652 >    info[i].ngroup = tot_groups;
1653      info[i].molMembershipArray = molMembershipArray;
1654    }
1655   }
# Line 1244 | Line 1660 | void SimSetup::mpiMolDivide(void){
1660    int i, j, k;
1661    int localMol, allMol;
1662    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1663 +  int local_rigid, local_groups;
1664 +  vector<int> globalMolIndex;
1665 +  int ncutgroups, atomsingroups, ngroupsinstamp;
1666 +  CutoffGroupStamp* cg;
1667  
1668    mpiSim = new mpiSimulation(info);
1669  
1670 <  globalIndex = mpiSim->divideLabor();
1670 >  mpiSim->divideLabor();
1671 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1672 >  globalGroupIndex = mpiSim->getGlobalGroupIndex();
1673 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1674  
1675    // set up the local variables
1676  
# Line 1260 | Line 1683 | void SimSetup::mpiMolDivide(void){
1683    local_bonds = 0;
1684    local_bends = 0;
1685    local_torsions = 0;
1686 <  globalAtomIndex = 0;
1687 <
1686 >  local_rigid = 0;
1687 >  local_groups = 0;
1688 >  globalAtomCounter = 0;
1689  
1690    for (i = 0; i < n_components; i++){
1691      for (j = 0; j < components_nmol[i]; j++){
# Line 1270 | Line 1694 | void SimSetup::mpiMolDivide(void){
1694          local_bonds += comp_stamps[i]->getNBonds();
1695          local_bends += comp_stamps[i]->getNBends();
1696          local_torsions += comp_stamps[i]->getNTorsions();
1697 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1698 +
1699 +        ncutgroups = comp_stamps[i]->getNCutoffGroups();
1700 +        atomsingroups = 0;
1701 +        for (k=0; k < ncutgroups; k++) {
1702 +          cg = comp_stamps[i]->getCutoffGroup(k);
1703 +          atomsingroups += cg->getNMembers();
1704 +        }
1705 +        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +
1706 +          ncutgroups;
1707 +        local_groups += ngroupsinstamp;    
1708 +
1709          localMol++;
1710        }      
1711        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1712 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1713 <        globalAtomIndex++;
1712 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1713 >        globalAtomCounter++;
1714        }
1715  
1716        allMol++;
# Line 1282 | Line 1718 | void SimSetup::mpiMolDivide(void){
1718    }
1719    local_SRI = local_bonds + local_bends + local_torsions;
1720  
1721 <  info[0].n_atoms = mpiSim->getMyNlocal();  
1722 <
1721 >  info[0].n_atoms = mpiSim->getNAtomsLocal();  
1722 >  
1723    if (local_atoms != info[0].n_atoms){
1724      sprintf(painCave.errMsg,
1725              "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
# Line 1293 | Line 1729 | void SimSetup::mpiMolDivide(void){
1729      simError();
1730    }
1731  
1732 +  info[0].ngroup = mpiSim->getNGroupsLocal();  
1733 +  if (local_groups != info[0].ngroup){
1734 +    sprintf(painCave.errMsg,
1735 +            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
1736 +            "\tlocalGroups (%d) are not equal.\n",
1737 +            info[0].ngroup, local_groups);
1738 +    painCave.isFatal = 1;
1739 +    simError();
1740 +  }
1741 +  
1742    info[0].n_bonds = local_bonds;
1743    info[0].n_bends = local_bends;
1744    info[0].n_torsions = local_torsions;
# Line 1315 | Line 1761 | void SimSetup::makeSysArrays(void){
1761  
1762    Atom** the_atoms;
1763    Molecule* the_molecules;
1318 –  Exclude** the_excludes;
1764  
1320 –
1765    for (l = 0; l < nInfo; l++){
1766      // create the atom and short range interaction arrays
1767  
# Line 1331 | Line 1775 | void SimSetup::makeSysArrays(void){
1775  
1776  
1777      molIndex = 0;
1778 <    for (i = 0; i < mpiSim->getTotNmol(); i++){
1778 >    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
1779        if (mol2proc[i] == worldRank){
1780          the_molecules[molIndex].setStampID(molCompType[i]);
1781          the_molecules[molIndex].setMyIndex(molIndex);
# Line 1343 | Line 1787 | void SimSetup::makeSysArrays(void){
1787   #else // is_mpi
1788  
1789      molIndex = 0;
1790 <    globalAtomIndex = 0;
1790 >    globalAtomCounter = 0;
1791      for (i = 0; i < n_components; i++){
1792        for (j = 0; j < components_nmol[i]; j++){
1793          the_molecules[molIndex].setStampID(i);
1794          the_molecules[molIndex].setMyIndex(molIndex);
1795          the_molecules[molIndex].setGlobalIndex(molIndex);
1796          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1797 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1798 <          globalAtomIndex++;
1797 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1798 >          globalAtomCounter++;
1799          }
1800          molIndex++;
1801        }
# Line 1360 | Line 1804 | void SimSetup::makeSysArrays(void){
1804  
1805   #endif // is_mpi
1806  
1807 <
1808 <    if (info[l].n_SRI){
1809 <      Exclude::createArray(info[l].n_SRI);
1366 <      the_excludes = new Exclude * [info[l].n_SRI];
1367 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1368 <        the_excludes[ex] = new Exclude(ex);
1369 <      }
1370 <      info[l].globalExcludes = new int;
1371 <      info[l].n_exclude = info[l].n_SRI;
1372 <    }
1373 <    else{
1374 <      Exclude::createArray(1);
1375 <      the_excludes = new Exclude * ;
1376 <      the_excludes[0] = new Exclude(0);
1377 <      the_excludes[0]->setPair(0, 0);
1378 <      info[l].globalExcludes = new int;
1379 <      info[l].globalExcludes[0] = 0;
1380 <      info[l].n_exclude = 0;
1381 <    }
1382 <
1807 >    info[l].globalExcludes = new int;
1808 >    info[l].globalExcludes[0] = 0;
1809 >    
1810      // set the arrays into the SimInfo object
1811  
1812      info[l].atoms = the_atoms;
1813      info[l].molecules = the_molecules;
1814      info[l].nGlobalExcludes = 0;
1815 <    info[l].excludes = the_excludes;
1389 <
1815 >    
1816      the_ff->setSimInfo(info);
1817    }
1818   }
# Line 1659 | Line 2085 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2085  
2086    theInfo.addProperty(zconsForcePolicy);
2087  
2088 +  //set zcons gap
2089 +  DoubleData* zconsGap = new DoubleData();
2090 +  zconsGap->setID(ZCONSGAP_ID);
2091 +
2092 +  if (globals->haveZConsGap()){
2093 +    zconsGap->setData(globals->getZconsGap());
2094 +    theInfo.addProperty(zconsGap);  
2095 +  }
2096 +
2097 +  //set zcons fixtime
2098 +  DoubleData* zconsFixtime = new DoubleData();
2099 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
2100 +
2101 +  if (globals->haveZConsFixTime()){
2102 +    zconsFixtime->setData(globals->getZconsFixtime());
2103 +    theInfo.addProperty(zconsFixtime);  
2104 +  }
2105 +
2106 +  //set zconsUsingSMD
2107 +  IntData* zconsUsingSMD = new IntData();
2108 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
2109 +
2110 +  if (globals->haveZConsUsingSMD()){
2111 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
2112 +    theInfo.addProperty(zconsUsingSMD);  
2113 +  }
2114 +
2115    //Determine the name of ouput file and add it into SimInfo's property list
2116    //Be careful, do not use inFileName, since it is a pointer which
2117    //point to a string at master node, and slave nodes do not contain that string
# Line 1688 | Line 2141 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2141      tempParaItem.zPos = zconStamp[i]->getZpos();
2142      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
2143      tempParaItem.kRatio = zconStamp[i]->getKratio();
2144 <
2144 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
2145 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
2146      zconsParaData->addItem(tempParaItem);
2147    }
2148  
# Line 1706 | Line 2160 | void SimSetup::setupZConstraint(SimInfo& theInfo){
2160    //push data into siminfo, therefore, we can retrieve later
2161    theInfo.addProperty(zconsParaData);
2162   }
2163 +
2164 + void SimSetup::makeMinimizer(){
2165 +
2166 +  OOPSEMinimizer* myOOPSEMinimizer;
2167 +  MinimizerParameterSet* param;
2168 +  char minimizerName[100];
2169 +  
2170 +  for (int i = 0; i < nInfo; i++){
2171 +    
2172 +    //prepare parameter set for minimizer
2173 +    param = new MinimizerParameterSet();
2174 +    param->setDefaultParameter();
2175 +
2176 +    if (globals->haveMinimizer()){
2177 +      param->setFTol(globals->getMinFTol());
2178 +    }
2179 +
2180 +    if (globals->haveMinGTol()){
2181 +      param->setGTol(globals->getMinGTol());
2182 +    }
2183 +
2184 +    if (globals->haveMinMaxIter()){
2185 +      param->setMaxIteration(globals->getMinMaxIter());
2186 +    }
2187 +
2188 +    if (globals->haveMinWriteFrq()){
2189 +      param->setMaxIteration(globals->getMinMaxIter());
2190 +    }
2191 +
2192 +    if (globals->haveMinWriteFrq()){
2193 +      param->setWriteFrq(globals->getMinWriteFrq());
2194 +    }
2195 +    
2196 +    if (globals->haveMinStepSize()){
2197 +      param->setStepSize(globals->getMinStepSize());
2198 +    }
2199 +
2200 +    if (globals->haveMinLSMaxIter()){
2201 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
2202 +    }    
2203 +
2204 +    if (globals->haveMinLSTol()){
2205 +      param->setLineSearchTol(globals->getMinLSTol());
2206 +    }    
2207 +
2208 +    strcpy(minimizerName, globals->getMinimizer());
2209 +
2210 +    if (!strcasecmp(minimizerName, "CG")){
2211 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
2212 +    }
2213 +    else if (!strcasecmp(minimizerName, "SD")){
2214 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
2215 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
2216 +    }
2217 +    else{
2218 +          sprintf(painCave.errMsg,
2219 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
2220 +          painCave.isFatal = 0;
2221 +          simError();
2222 +
2223 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
2224 +    }
2225 +     info[i].the_integrator = myOOPSEMinimizer;
2226 +
2227 +     //store the minimizer into simInfo
2228 +     info[i].the_minimizer = myOOPSEMinimizer;
2229 +     info[i].has_minimizer = true;
2230 +  }
2231 +
2232 + }

Diff Legend

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