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 1066 by tim, Tue Feb 24 16:36:33 2004 UTC vs.
Revision 1204 by gezelter, Thu May 27 19:26:42 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 "ConjugateMinimizer.hpp"
14   #include "OOPSEMinimizer.hpp"
15  
# Line 146 | Line 147 | void SimSetup::createSim(void){
147    // make the output filenames
148  
149    makeOutNames();
149 –
150 –  if (globals->haveMinimizer())
151 –    // make minimizer
152 –    makeMinimizer();
153 –  else
154 –    // make the integrator
155 –    makeIntegrator();
150    
151   #ifdef IS_MPI
152    mpiSim->mpiRefresh();
# Line 161 | Line 155 | void SimSetup::createSim(void){
155    // initialize the Fortran
156  
157    initFortran();
158 +
159 +  if (globals->haveMinimizer())
160 +    // make minimizer
161 +    makeMinimizer();
162 +  else
163 +    // make the integrator
164 +    makeIntegrator();
165 +
166   }
167  
168  
169   void SimSetup::makeMolecules(void){
170 <  int k;
171 <  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
170 >  int i, j, k;
171 >  int exI, exJ, exK, exL, slI, slJ;
172 >  int tempI, tempJ, tempK, tempL;
173 >  int molI;
174 >  int stampID, atomOffset, rbOffset;
175    molInit molInfo;
176    DirectionalAtom* dAtom;
177 +  RigidBody* myRB;
178 +  StuntDouble* mySD;
179    LinkedAssign* extras;
180    LinkedAssign* current_extra;
181    AtomStamp* currentAtom;
182    BondStamp* currentBond;
183    BendStamp* currentBend;
184    TorsionStamp* currentTorsion;
185 +  RigidBodyStamp* currentRigidBody;
186 +  CutoffGroupStamp* currentCutoffGroup;
187 +  CutoffGroup* myCutoffGroup;
188 +  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
189 +  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
190  
191    bond_pair* theBonds;
192    bend_set* theBends;
193    torsion_set* theTorsions;
194  
195 +  set<int> skipList;
196 +
197 +  double phi, theta, psi;
198 +  char* molName;
199 +  char rbName[100];
200 +
201    //init the forceField paramters
202  
203    the_ff->readParams();
204  
187 –
205    // init the atoms
206  
207 <  double phi, theta, psi;
191 <  double sux, suy, suz;
192 <  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
193 <  double ux, uy, uz, u, uSqr;
207 >  int nMembers, nNew, rb1, rb2;
208  
209    for (k = 0; k < nInfo; k++){
210      the_ff->setSimInfo(&(info[k]));
211  
212      atomOffset = 0;
213 <    excludeOffset = 0;
213 >
214      for (i = 0; i < info[k].n_mol; i++){
215        stampID = info[k].molecules[i].getStampID();
216 +      molName = comp_stamps[stampID]->getID();
217  
218        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
219        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
220        molInfo.nBends = comp_stamps[stampID]->getNBends();
221        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
222 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
222 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
223  
224 +      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
225 +      
226        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
210 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
211 –      molInfo.myBonds = new Bond * [molInfo.nBonds];
212 –      molInfo.myBends = new Bend * [molInfo.nBends];
213 –      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
227  
228 +      if (molInfo.nBonds > 0)
229 +        molInfo.myBonds = new Bond*[molInfo.nBonds];
230 +      else
231 +        molInfo.myBonds = NULL;
232 +
233 +      if (molInfo.nBends > 0)
234 +        molInfo.myBends = new Bend*[molInfo.nBends];
235 +      else
236 +        molInfo.myBends = NULL;
237 +
238 +      if (molInfo.nTorsions > 0)
239 +        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
240 +      else
241 +        molInfo.myTorsions = NULL;
242 +
243        theBonds = new bond_pair[molInfo.nBonds];
244        theBends = new bend_set[molInfo.nBends];
245        theTorsions = new torsion_set[molInfo.nTorsions];
246 <
246 >      
247        // make the Atoms
248  
249        for (j = 0; j < molInfo.nAtoms; j++){
250          currentAtom = comp_stamps[stampID]->getAtom(j);
251 +
252          if (currentAtom->haveOrientation()){
253            dAtom = new DirectionalAtom((j + atomOffset),
254                                        info[k].getConfiguration());
# Line 233 | Line 262 | void SimSetup::makeMolecules(void){
262            phi = currentAtom->getEulerPhi() * M_PI / 180.0;
263            theta = currentAtom->getEulerTheta() * M_PI / 180.0;
264            psi = currentAtom->getEulerPsi()* M_PI / 180.0;
265 +
266 +          dAtom->setUnitFrameFromEuler(phi, theta, psi);
267              
268 <          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
269 <          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
239 <          Axz = sin(theta) * sin(psi);
240 <          
241 <          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
242 <          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
243 <          Ayz = sin(theta) * cos(psi);
244 <          
245 <          Azx = sin(phi) * sin(theta);
246 <          Azy = -cos(phi) * sin(theta);
247 <          Azz = cos(theta);
268 >        }
269 >        else{
270  
271 <          sux = 0.0;
250 <          suy = 0.0;
251 <          suz = 1.0;
271 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
272  
273 <          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
254 <          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
255 <          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
256 <
257 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
258 <
259 <          u = sqrt(uSqr);
260 <          ux = ux / u;
261 <          uy = uy / u;
262 <          uz = uz / u;
273 >        }
274  
264 –          dAtom->setSUx(ux);
265 –          dAtom->setSUy(uy);
266 –          dAtom->setSUz(uz);
267 –        }
268 –        else{
269 –          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
270 –                                               info[k].getConfiguration());
271 –        }
275          molInfo.myAtoms[j]->setType(currentAtom->getType());
273 –
276   #ifdef IS_MPI
277  
278 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
278 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
279  
280   #endif // is_mpi
281        }
# Line 284 | Line 286 | void SimSetup::makeMolecules(void){
286          theBonds[j].a = currentBond->getA() + atomOffset;
287          theBonds[j].b = currentBond->getB() + atomOffset;
288  
289 <        exI = theBonds[j].a;
290 <        exJ = theBonds[j].b;
289 >        tempI = theBonds[j].a;
290 >        tempJ = theBonds[j].b;
291  
290 –        // exclude_I must always be the smaller of the pair
291 –        if (exI > exJ){
292 –          tempEx = exI;
293 –          exI = exJ;
294 –          exJ = tempEx;
295 –        }
292   #ifdef IS_MPI
293 <        tempEx = exI;
294 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
295 <        tempEx = exJ;
296 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
293 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
294 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
295 > #else
296 >        exI = tempI + 1;
297 >        exJ = tempJ + 1;
298 > #endif
299  
300 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
303 < #else  // isn't MPI
304 <
305 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
306 < #endif  //is_mpi
300 >        info[k].excludes->addPair(exI, exJ);
301        }
308 –      excludeOffset += molInfo.nBonds;
302  
303        //make the bends
304        for (j = 0; j < molInfo.nBends; j++){
# Line 355 | Line 348 | void SimSetup::makeMolecules(void){
348            }
349          }
350  
351 <        if (!theBends[j].isGhost){
352 <          exI = theBends[j].a;
353 <          exJ = theBends[j].c;
354 <        }
355 <        else{
363 <          exI = theBends[j].a;
364 <          exJ = theBends[j].b;
365 <        }
366 <
367 <        // exclude_I must always be the smaller of the pair
368 <        if (exI > exJ){
369 <          tempEx = exI;
370 <          exI = exJ;
371 <          exJ = tempEx;
372 <        }
351 >        if (theBends[j].isGhost) {
352 >          
353 >          tempI = theBends[j].a;
354 >          tempJ = theBends[j].b;
355 >          
356   #ifdef IS_MPI
357 <        tempEx = exI;
358 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
359 <        tempEx = exJ;
360 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
357 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
358 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
359 > #else
360 >          exI = tempI + 1;
361 >          exJ = tempJ + 1;
362 > #endif          
363 >          info[k].excludes->addPair(exI, exJ);
364  
365 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
366 < #else  // isn't MPI
367 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
368 < #endif  //is_mpi
365 >        } else {
366 >
367 >          tempI = theBends[j].a;
368 >          tempJ = theBends[j].b;
369 >          tempK = theBends[j].c;
370 >          
371 > #ifdef IS_MPI
372 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
373 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
374 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
375 > #else
376 >          exI = tempI + 1;
377 >          exJ = tempJ + 1;
378 >          exK = tempK + 1;
379 > #endif
380 >          
381 >          info[k].excludes->addPair(exI, exK);
382 >          info[k].excludes->addPair(exI, exJ);
383 >          info[k].excludes->addPair(exJ, exK);
384 >        }
385        }
384 –      excludeOffset += molInfo.nBends;
386  
387        for (j = 0; j < molInfo.nTorsions; j++){
388          currentTorsion = comp_stamps[stampID]->getTorsion(j);
# Line 390 | Line 391 | void SimSetup::makeMolecules(void){
391          theTorsions[j].c = currentTorsion->getC() + atomOffset;
392          theTorsions[j].d = currentTorsion->getD() + atomOffset;
393  
394 <        exI = theTorsions[j].a;
395 <        exJ = theTorsions[j].d;
394 >        tempI = theTorsions[j].a;      
395 >        tempJ = theTorsions[j].b;
396 >        tempK = theTorsions[j].c;
397 >        tempL = theTorsions[j].d;
398  
399 <        // exclude_I must always be the smaller of the pair
400 <        if (exI > exJ){
401 <          tempEx = exI;
402 <          exI = exJ;
403 <          exJ = tempEx;
399 > #ifdef IS_MPI
400 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
401 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
402 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
403 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
404 > #else
405 >        exI = tempI + 1;
406 >        exJ = tempJ + 1;
407 >        exK = tempK + 1;
408 >        exL = tempL + 1;
409 > #endif
410 >
411 >        info[k].excludes->addPair(exI, exJ);
412 >        info[k].excludes->addPair(exI, exK);
413 >        info[k].excludes->addPair(exI, exL);        
414 >        info[k].excludes->addPair(exJ, exK);
415 >        info[k].excludes->addPair(exJ, exL);
416 >        info[k].excludes->addPair(exK, exL);
417 >      }
418 >
419 >      
420 >      molInfo.myRigidBodies.clear();
421 >      
422 >      for (j = 0; j < molInfo.nRigidBodies; j++){
423 >
424 >        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
425 >        nMembers = currentRigidBody->getNMembers();
426 >
427 >        // Create the Rigid Body:
428 >
429 >        myRB = new RigidBody();
430 >
431 >        sprintf(rbName,"%s_RB_%d", molName, j);
432 >        myRB->setType(rbName);
433 >        
434 >        for (rb1 = 0; rb1 < nMembers; rb1++) {
435 >
436 >          // molI is atom numbering inside this molecule
437 >          molI = currentRigidBody->getMember(rb1);    
438 >
439 >          // tempI is atom numbering on local processor
440 >          tempI = molI + atomOffset;
441 >
442 >          // currentAtom is the AtomStamp (which we need for
443 >          // rigid body reference positions)
444 >          currentAtom = comp_stamps[stampID]->getAtom(molI);
445 >
446 >          // When we add to the rigid body, add the atom itself and
447 >          // the stamp info:
448 >
449 >          myRB->addAtom(info[k].atoms[tempI], currentAtom);
450 >          
451 >          // Add this atom to the Skip List for the integrators
452 > #ifdef IS_MPI
453 >          slI = info[k].atoms[tempI]->getGlobalIndex();
454 > #else
455 >          slI = tempI;
456 > #endif
457 >          skipList.insert(slI);
458 >          
459          }
460 +        
461 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
462 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
463 +            
464 +            tempI = currentRigidBody->getMember(rb1);
465 +            tempJ = currentRigidBody->getMember(rb2);
466 +            
467 +            // Some explanation is required here.
468 +            // Fortran indexing starts at 1, while c indexing starts at 0
469 +            // Also, in parallel computations, the GlobalIndex is
470 +            // used for the exclude list:
471 +            
472   #ifdef IS_MPI
473 <        tempEx = exI;
474 <        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
475 <        tempEx = exJ;
476 <        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
473 >            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
474 >            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
475 > #else
476 >            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
477 >            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
478 > #endif
479 >            
480 >            info[k].excludes->addPair(exI, exJ);
481 >            
482 >          }
483 >        }
484  
485 <        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
486 < #else  // isn't MPI
410 <        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
411 < #endif  //is_mpi
485 >        molInfo.myRigidBodies.push_back(myRB);
486 >        info[k].rigidBodies.push_back(myRB);
487        }
488 <      excludeOffset += molInfo.nTorsions;
488 >      
489  
490 +      //create cutoff group for molecule
491  
492 <      // send the arrays off to the forceField for init.
492 >      cutoffAtomSet.clear();
493 >      molInfo.myCutoffGroups.clear();
494 >      
495 >      for (j = 0; j < nCutoffGroups; j++){
496 >
497 >        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
498 >        nMembers = currentCutoffGroup->getNMembers();
499 >
500 >        myCutoffGroup = new CutoffGroup();
501 >        
502 >        for (int cg = 0; cg < nMembers; cg++) {
503 >
504 >          // molI is atom numbering inside this molecule
505 >          molI = currentCutoffGroup->getMember(cg);    
506 >
507 >          // tempI is atom numbering on local processor
508 >          tempI = molI + atomOffset;
509 >          
510 >          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
511 >
512 >          cutoffAtomSet.insert(tempI);
513 >        }
514 >
515 >        molInfo.myCutoffGroups.push_back(myCutoffGroup);
516 >      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
517 >
518 >      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
519 >
520 >      for(j = 0; j < molInfo.nAtoms; j++){
521 >
522 >        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
523 >          myCutoffGroup = new CutoffGroup();
524 >          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
525 >          molInfo.myCutoffGroups.push_back(myCutoffGroup);
526 >        }
527 >          
528 >      }
529 >
530 >              
531 >
532 >
533 >      // After this is all set up, scan through the atoms to
534 >      // see if they can be added to the integrableObjects:
535 >
536 >      molInfo.myIntegrableObjects.clear();
537 >      
538 >
539 >      for (j = 0; j < molInfo.nAtoms; j++){
540 >
541 > #ifdef IS_MPI
542 >        slJ = molInfo.myAtoms[j]->getGlobalIndex();
543 > #else
544 >        slJ = j+atomOffset;
545 > #endif
546  
547 +        // if they aren't on the skip list, then they can be integrated
548 +
549 +        if (skipList.find(slJ) == skipList.end()) {
550 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
551 +          info[k].integrableObjects.push_back(mySD);
552 +          molInfo.myIntegrableObjects.push_back(mySD);
553 +        }
554 +      }
555 +
556 +      // all rigid bodies are integrated:
557 +
558 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
559 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
560 +        info[k].integrableObjects.push_back(mySD);      
561 +        molInfo.myIntegrableObjects.push_back(mySD);
562 +      }
563 +    
564 +      
565 +      // send the arrays off to the forceField for init.
566 +      
567        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
568        the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
569        the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
570        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
571                                   theTorsions);
572  
424 –
573        info[k].molecules[i].initialize(molInfo);
574  
575  
# Line 429 | Line 577 | void SimSetup::makeMolecules(void){
577        delete[] theBonds;
578        delete[] theBends;
579        delete[] theTorsions;
580 <    }
580 >    }    
581    }
582  
583   #ifdef IS_MPI
584    sprintf(checkPointMsg, "all molecules initialized succesfully");
585    MPIcheckPoint();
586   #endif // is_mpi
439 –
440 –  // clean up the forcefield
587  
442 –  the_ff->calcRcut();
443 –  the_ff->cleanMe();
588   }
589  
590   void SimSetup::initFromBass(void){
# Line 727 | Line 871 | void SimSetup::gatherInfo(void){
871    }
872  
873    //check whether sample time, status time, thermal time and reset time are divisble by dt
874 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
874 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
875      sprintf(painCave.errMsg,
876              "Sample time is not divisible by dt.\n"
877              "\tThis will result in samples that are not uniformly\n"
# Line 737 | Line 881 | void SimSetup::gatherInfo(void){
881      simError();    
882    }
883  
884 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
884 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
885      sprintf(painCave.errMsg,
886              "Status time is not divisible by dt.\n"
887              "\tThis will result in status reports that are not uniformly\n"
# Line 773 | Line 917 | void SimSetup::gatherInfo(void){
917      if (globals->haveSampleTime()){
918        info[i].sampleTime = globals->getSampleTime();
919        info[i].statusTime = info[i].sampleTime;
776 –      info[i].thermalTime = info[i].sampleTime;
920      }
921      else{
922        info[i].sampleTime = globals->getRunTime();
923        info[i].statusTime = info[i].sampleTime;
781 –      info[i].thermalTime = info[i].sampleTime;
924      }
925  
926      if (globals->haveStatusTime()){
# Line 787 | Line 929 | void SimSetup::gatherInfo(void){
929  
930      if (globals->haveThermalTime()){
931        info[i].thermalTime = globals->getThermalTime();
932 +    } else {
933 +      info[i].thermalTime = globals->getRunTime();
934      }
935  
936      info[i].resetIntegrator = 0;
# Line 804 | Line 948 | void SimSetup::gatherInfo(void){
948  
949      info[i].useInitXSstate = globals->getUseInitXSstate();
950      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
951 <    
951 >
952 >    // check for thermodynamic integration
953 >    if (globals->getUseThermInt()) {
954 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
955 >        info[i].useThermInt = globals->getUseThermInt();
956 >        info[i].thermIntLambda = globals->getThermIntLambda();
957 >        info[i].thermIntK = globals->getThermIntK();
958 >        
959 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
960 >        info[i].restraint = myRestraint;
961 >      }
962 >      else {
963 >        sprintf(painCave.errMsg,
964 >                "SimSetup Error:\n"
965 >                "\tKeyword useThermInt was set to 'true' but\n"
966 >                "\tthermodynamicIntegrationLambda (and/or\n"
967 >                "\tthermodynamicIntegrationK) was not specified.\n"
968 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
969 >        painCave.isFatal = 1;
970 >        simError();    
971 >      }
972 >    }
973 >    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
974 >        sprintf(painCave.errMsg,
975 >                "SimSetup Warning: If you want to use Thermodynamic\n"
976 >                "\tIntegration, set useThermInt to 'true' in your .bass file.\n"
977 >                "\tThe useThermInt keyword is 'false' by default, so your\n"
978 >                "\tlambda and/or k values are being ignored.\n");
979 >        painCave.isFatal = 0;
980 >        simError();  
981 >    }
982    }
983    
984    //setup seed for random number generator
# Line 857 | Line 1031 | void SimSetup::finalInfoCheck(void){
1031   void SimSetup::finalInfoCheck(void){
1032    int index;
1033    int usesDipoles;
1034 +  int usesCharges;
1035    int i;
1036  
1037    for (i = 0; i < nInfo; i++){
# Line 868 | Line 1043 | void SimSetup::finalInfoCheck(void){
1043        usesDipoles = (info[i].atoms[index])->hasDipole();
1044        index++;
1045      }
1046 <
1046 >    index = 0;
1047 >    usesCharges = 0;
1048 >    while ((index < info[i].n_atoms) && !usesCharges){
1049 >      usesCharges= (info[i].atoms[index])->hasCharge();
1050 >      index++;
1051 >    }
1052   #ifdef IS_MPI
1053      int myUse = usesDipoles;
1054      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1055   #endif //is_mpi
1056  
1057 <    double theEcr, theEst;
1057 >    double theRcut, theRsw;
1058  
1059 <    if (globals->getUseRF()){
1060 <      info[i].useReactionField = 1;
1059 >    if (globals->haveRcut()) {
1060 >      theRcut = globals->getRcut();
1061  
1062 <      if (!globals->haveECR()){
1063 <        sprintf(painCave.errMsg,
1064 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1062 >      if (globals->haveRsw())
1063 >        theRsw = globals->getRsw();
1064 >      else
1065 >        theRsw = theRcut;
1066 >      
1067 >      info[i].setDefaultRcut(theRcut, theRsw);
1068 >
1069 >    } else {
1070 >      
1071 >      the_ff->calcRcut();
1072 >      theRcut = info[i].getRcut();
1073 >
1074 >      if (globals->haveRsw())
1075 >        theRsw = globals->getRsw();
1076 >      else
1077 >        theRsw = theRcut;
1078 >      
1079 >      info[i].setDefaultRcut(theRcut, theRsw);
1080 >    }
1081 >
1082 >    if (globals->getUseRF()){
1083 >      info[i].useReactionField = 1;
1084 >      
1085 >      if (!globals->haveRcut()){
1086 >        sprintf(painCave.errMsg,
1087 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1088                  "\tOOPSE will use a default value of 15.0 angstroms"
1089 <                "\tfor the electrostaticCutoffRadius.\n");
1089 >                "\tfor the cutoffRadius.\n");
1090          painCave.isFatal = 0;
1091          simError();
1092 <        theEcr = 15.0;
1092 >        theRcut = 15.0;
1093        }
1094        else{
1095 <        theEcr = globals->getECR();
1095 >        theRcut = globals->getRcut();
1096        }
1097  
1098 <      if (!globals->haveEST()){
1098 >      if (!globals->haveRsw()){
1099          sprintf(painCave.errMsg,
1100 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1100 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1101                  "\tOOPSE will use a default value of\n"
1102 <                "\t0.05 * electrostaticCutoffRadius\n"
900 <                "\tfor the electrostaticSkinThickness\n");
1102 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1103          painCave.isFatal = 0;
1104          simError();
1105 <        theEst = 0.05 * theEcr;
1105 >        theRsw = 0.95 * theRcut;
1106        }
1107        else{
1108 <        theEst = globals->getEST();
1108 >        theRsw = globals->getRsw();
1109        }
1110  
1111 <      info[i].setDefaultEcr(theEcr, theEst);
1111 >      info[i].setDefaultRcut(theRcut, theRsw);
1112  
1113        if (!globals->haveDielectric()){
1114          sprintf(painCave.errMsg,
# Line 919 | Line 1121 | void SimSetup::finalInfoCheck(void){
1121        info[i].dielectric = globals->getDielectric();
1122      }
1123      else{
1124 <      if (usesDipoles){
1125 <        if (!globals->haveECR()){
1124 >      if (usesDipoles || usesCharges){
1125 >
1126 >        if (!globals->haveRcut()){
1127            sprintf(painCave.errMsg,
1128 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1128 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1129                    "\tOOPSE will use a default value of 15.0 angstroms"
1130 <                  "\tfor the electrostaticCutoffRadius.\n");
1131 <          painCave.isFatal = 0;
1132 <          simError();
1133 <          theEcr = 15.0;
1134 <        }
1130 >                  "\tfor the cutoffRadius.\n");
1131 >          painCave.isFatal = 0;
1132 >          simError();
1133 >          theRcut = 15.0;
1134 >      }
1135          else{
1136 <          theEcr = globals->getECR();
1136 >          theRcut = globals->getRcut();
1137          }
1138 <        
1139 <        if (!globals->haveEST()){
1138 >        
1139 >        if (!globals->haveRsw()){
1140            sprintf(painCave.errMsg,
1141 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1141 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1142                    "\tOOPSE will use a default value of\n"
1143 <                  "\t0.05 * electrostaticCutoffRadius\n"
941 <                  "\tfor the electrostaticSkinThickness\n");
1143 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1144            painCave.isFatal = 0;
1145            simError();
1146 <          theEst = 0.05 * theEcr;
1146 >          theRsw = 0.95 * theRcut;
1147          }
1148          else{
1149 <          theEst = globals->getEST();
1149 >          theRsw = globals->getRsw();
1150          }
1151 +        
1152 +        info[i].setDefaultRcut(theRcut, theRsw);
1153          
950 –        info[i].setDefaultEcr(theEcr, theEst);
1154        }
1155      }
1156    }
# Line 955 | Line 1158 | void SimSetup::finalInfoCheck(void){
1158    strcpy(checkPointMsg, "post processing checks out");
1159    MPIcheckPoint();
1160   #endif // is_mpi
1161 +
1162 +  // clean up the forcefield
1163 +  the_ff->cleanMe();
1164   }
1165    
1166   void SimSetup::initSystemCoords(void){
# Line 1082 | Line 1288 | void SimSetup::makeOutNames(void){
1288          }
1289          else{
1290            strcat(info[k].statusName, ".stat");
1291 +        }
1292 +      }
1293 +
1294 +      strcpy(info[k].rawPotName, inFileName);
1295 +      nameLength = strlen(info[k].rawPotName);
1296 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1297 +      if (!strcmp(endTest, ".bass")){
1298 +        strcpy(endTest, ".raw");
1299 +      }
1300 +      else if (!strcmp(endTest, ".BASS")){
1301 +        strcpy(endTest, ".raw");
1302 +      }
1303 +      else{
1304 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1305 +        if (!strcmp(endTest, ".bss")){
1306 +          strcpy(endTest, ".raw");
1307 +        }
1308 +        else if (!strcmp(endTest, ".mdl")){
1309 +          strcpy(endTest, ".raw");
1310 +        }
1311 +        else{
1312 +          strcat(info[k].rawPotName, ".raw");
1313          }
1314        }
1315  
# Line 1169 | Line 1397 | void SimSetup::compList(void){
1397    LinkedMolStamp* headStamp = new LinkedMolStamp();
1398    LinkedMolStamp* currentStamp = NULL;
1399    comp_stamps = new MoleculeStamp * [n_components];
1400 +  bool haveCutoffGroups;
1401  
1402 +  haveCutoffGroups = false;
1403 +  
1404    // make an array of molecule stamps that match the components used.
1405    // also extract the used stamps out into a separate linked list
1406  
# Line 1204 | Line 1435 | void SimSetup::compList(void){
1435        headStamp->add(currentStamp);
1436        comp_stamps[i] = headStamp->match(id);
1437      }
1438 +
1439 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1440 +      haveCutoffGroups = true;    
1441    }
1442 +    
1443 +  for (i = 0; i < nInfo; i++)
1444 +    info[i].haveCutoffGroups = haveCutoffGroups;
1445  
1446   #ifdef IS_MPI
1447    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1213 | Line 1450 | void SimSetup::calcSysValues(void){
1450   }
1451  
1452   void SimSetup::calcSysValues(void){
1453 <  int i;
1453 >  int i, j;
1454 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1455  
1456    int* molMembershipArray;
1457 +  CutoffGroupStamp* cg;
1458  
1459    tot_atoms = 0;
1460    tot_bonds = 0;
1461    tot_bends = 0;
1462    tot_torsions = 0;
1463 +  tot_rigid = 0;
1464 +  tot_groups = 0;
1465    for (i = 0; i < n_components; i++){
1466      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1467      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1468      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1469      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1470 <  }
1470 >    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1471  
1472 +    ncutgroups = comp_stamps[i]->getNCutoffGroups();
1473 +    atomsingroups = 0;
1474 +    for (j=0; j < ncutgroups; j++) {
1475 +      cg = comp_stamps[i]->getCutoffGroup(j);
1476 +      atomsingroups += cg->getNMembers();
1477 +    }
1478 +    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
1479 +    tot_groups += components_nmol[i] * ngroupsinstamp;    
1480 +  }
1481 +  
1482    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1483    molMembershipArray = new int[tot_atoms];
1484  
# Line 1238 | Line 1489 | void SimSetup::calcSysValues(void){
1489      info[i].n_torsions = tot_torsions;
1490      info[i].n_SRI = tot_SRI;
1491      info[i].n_mol = tot_nmol;
1492 <
1492 >    info[i].ngroup = tot_groups;
1493      info[i].molMembershipArray = molMembershipArray;
1494    }
1495   }
# Line 1249 | Line 1500 | void SimSetup::mpiMolDivide(void){
1500    int i, j, k;
1501    int localMol, allMol;
1502    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1503 +  int local_rigid, local_groups;
1504 +  vector<int> globalMolIndex;
1505 +  int ncutgroups, atomsingroups, ngroupsinstamp;
1506 +  CutoffGroupStamp* cg;
1507  
1508    mpiSim = new mpiSimulation(info);
1509  
1510 <  globalIndex = mpiSim->divideLabor();
1510 >  mpiSim->divideLabor();
1511 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1512 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1513  
1514    // set up the local variables
1515  
# Line 1265 | Line 1522 | void SimSetup::mpiMolDivide(void){
1522    local_bonds = 0;
1523    local_bends = 0;
1524    local_torsions = 0;
1525 <  globalAtomIndex = 0;
1525 >  local_rigid = 0;
1526 >  local_groups = 0;
1527 >  globalAtomCounter = 0;
1528  
1270 –
1529    for (i = 0; i < n_components; i++){
1530      for (j = 0; j < components_nmol[i]; j++){
1531        if (mol2proc[allMol] == worldRank){
# Line 1275 | Line 1533 | void SimSetup::mpiMolDivide(void){
1533          local_bonds += comp_stamps[i]->getNBonds();
1534          local_bends += comp_stamps[i]->getNBends();
1535          local_torsions += comp_stamps[i]->getNTorsions();
1536 +        local_rigid += comp_stamps[i]->getNRigidBodies();
1537 +
1538 +        ncutgroups = comp_stamps[i]->getNCutoffGroups();
1539 +        atomsingroups = 0;
1540 +        for (k=0; k < ncutgroups; k++) {
1541 +          cg = comp_stamps[i]->getCutoffGroup(k);
1542 +          atomsingroups += cg->getNMembers();
1543 +        }
1544 +        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +
1545 +          ncutgroups;
1546 +        local_groups += ngroupsinstamp;    
1547 +
1548          localMol++;
1549        }      
1550        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1551 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1552 <        globalAtomIndex++;
1551 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1552 >        globalAtomCounter++;
1553        }
1554  
1555        allMol++;
# Line 1287 | Line 1557 | void SimSetup::mpiMolDivide(void){
1557    }
1558    local_SRI = local_bonds + local_bends + local_torsions;
1559  
1560 <  info[0].n_atoms = mpiSim->getMyNlocal();  
1561 <
1560 >  info[0].n_atoms = mpiSim->getNAtomsLocal();  
1561 >  
1562    if (local_atoms != info[0].n_atoms){
1563      sprintf(painCave.errMsg,
1564              "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
# Line 1298 | Line 1568 | void SimSetup::mpiMolDivide(void){
1568      simError();
1569    }
1570  
1571 +  info[0].ngroup = mpiSim->getNGroupsLocal();  
1572 +  if (local_groups != info[0].ngroup){
1573 +    sprintf(painCave.errMsg,
1574 +            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
1575 +            "\tlocalGroups (%d) are not equal.\n",
1576 +            info[0].ngroup, local_groups);
1577 +    painCave.isFatal = 1;
1578 +    simError();
1579 +  }
1580 +  
1581    info[0].n_bonds = local_bonds;
1582    info[0].n_bends = local_bends;
1583    info[0].n_torsions = local_torsions;
# Line 1320 | Line 1600 | void SimSetup::makeSysArrays(void){
1600  
1601    Atom** the_atoms;
1602    Molecule* the_molecules;
1323 –  Exclude** the_excludes;
1603  
1325 –
1604    for (l = 0; l < nInfo; l++){
1605      // create the atom and short range interaction arrays
1606  
# Line 1336 | Line 1614 | void SimSetup::makeSysArrays(void){
1614  
1615  
1616      molIndex = 0;
1617 <    for (i = 0; i < mpiSim->getTotNmol(); i++){
1617 >    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
1618        if (mol2proc[i] == worldRank){
1619          the_molecules[molIndex].setStampID(molCompType[i]);
1620          the_molecules[molIndex].setMyIndex(molIndex);
# Line 1348 | Line 1626 | void SimSetup::makeSysArrays(void){
1626   #else // is_mpi
1627  
1628      molIndex = 0;
1629 <    globalAtomIndex = 0;
1629 >    globalAtomCounter = 0;
1630      for (i = 0; i < n_components; i++){
1631        for (j = 0; j < components_nmol[i]; j++){
1632          the_molecules[molIndex].setStampID(i);
1633          the_molecules[molIndex].setMyIndex(molIndex);
1634          the_molecules[molIndex].setGlobalIndex(molIndex);
1635          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1636 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1637 <          globalAtomIndex++;
1636 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1637 >          globalAtomCounter++;
1638          }
1639          molIndex++;
1640        }
# Line 1365 | Line 1643 | void SimSetup::makeSysArrays(void){
1643  
1644   #endif // is_mpi
1645  
1646 <
1647 <    if (info[l].n_SRI){
1648 <      Exclude::createArray(info[l].n_SRI);
1371 <      the_excludes = new Exclude * [info[l].n_SRI];
1372 <      for (int ex = 0; ex < info[l].n_SRI; ex++){
1373 <        the_excludes[ex] = new Exclude(ex);
1374 <      }
1375 <      info[l].globalExcludes = new int;
1376 <      info[l].n_exclude = info[l].n_SRI;
1377 <    }
1378 <    else{
1379 <      Exclude::createArray(1);
1380 <      the_excludes = new Exclude * ;
1381 <      the_excludes[0] = new Exclude(0);
1382 <      the_excludes[0]->setPair(0, 0);
1383 <      info[l].globalExcludes = new int;
1384 <      info[l].globalExcludes[0] = 0;
1385 <      info[l].n_exclude = 0;
1386 <    }
1387 <
1646 >    info[l].globalExcludes = new int;
1647 >    info[l].globalExcludes[0] = 0;
1648 >    
1649      // set the arrays into the SimInfo object
1650  
1651      info[l].atoms = the_atoms;
1652      info[l].molecules = the_molecules;
1653      info[l].nGlobalExcludes = 0;
1654 <    info[l].excludes = the_excludes;
1394 <
1654 >    
1655      the_ff->setSimInfo(info);
1656    }
1657   }
# Line 1664 | Line 1924 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1924  
1925    theInfo.addProperty(zconsForcePolicy);
1926  
1927 +  //set zcons gap
1928 +  DoubleData* zconsGap = new DoubleData();
1929 +  zconsGap->setID(ZCONSGAP_ID);
1930 +
1931 +  if (globals->haveZConsGap()){
1932 +    zconsGap->setData(globals->getZconsGap());
1933 +    theInfo.addProperty(zconsGap);  
1934 +  }
1935 +
1936 +  //set zcons fixtime
1937 +  DoubleData* zconsFixtime = new DoubleData();
1938 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
1939 +
1940 +  if (globals->haveZConsFixTime()){
1941 +    zconsFixtime->setData(globals->getZconsFixtime());
1942 +    theInfo.addProperty(zconsFixtime);  
1943 +  }
1944 +
1945 +  //set zconsUsingSMD
1946 +  IntData* zconsUsingSMD = new IntData();
1947 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1948 +
1949 +  if (globals->haveZConsUsingSMD()){
1950 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1951 +    theInfo.addProperty(zconsUsingSMD);  
1952 +  }
1953 +
1954    //Determine the name of ouput file and add it into SimInfo's property list
1955    //Be careful, do not use inFileName, since it is a pointer which
1956    //point to a string at master node, and slave nodes do not contain that string
# Line 1693 | Line 1980 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1980      tempParaItem.zPos = zconStamp[i]->getZpos();
1981      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1982      tempParaItem.kRatio = zconStamp[i]->getKratio();
1983 <
1983 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1984 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1985      zconsParaData->addItem(tempParaItem);
1986    }
1987  

Diff Legend

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