ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/SimInfo.cpp
(Generate patch)

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 749 by tim, Wed Nov 16 23:10:02 2005 UTC

# Line 48 | Line 48
48  
49   #include <algorithm>
50   #include <set>
51 + #include <map>
52  
53   #include "brains/SimInfo.hpp"
54   #include "math/Vector3.hpp"
55   #include "primitives/Molecule.hpp"
56 + #include "UseTheForce/fCutoffPolicy.h"
57 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
58 + #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
59 + #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
60   #include "UseTheForce/doForces_interface.h"
61 + #include "UseTheForce/DarkSide/electrostatic_interface.h"
62   #include "UseTheForce/notifyCutoffs_interface.h"
63 + #include "UseTheForce/DarkSide/switcheroo_interface.h"
64   #include "utils/MemoryUtils.hpp"
65   #include "utils/simError.h"
66   #include "selection/SelectionManager.hpp"
# Line 64 | Line 71 | namespace oopse {
71   #endif
72  
73   namespace oopse {
74 +  std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) {
75 +    std::map<int, std::set<int> >::iterator i = container.find(index);
76 +    std::set<int> result;
77 +    if (i != container.end()) {
78 +        result = i->second;
79 +    }
80  
81 +    return result;
82 +  }
83 +  
84    SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
85                     ForceField* ff, Globals* simParams) :
86      stamps_(stamps), forceField_(ff), simParams_(simParams),
# Line 80 | Line 96 | namespace oopse {
96        MoleculeStamp* molStamp;
97        int nMolWithSameStamp;
98        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
99 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
99 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
100        CutoffGroupStamp* cgStamp;    
101        RigidBodyStamp* rbStamp;
102        int nRigidAtoms = 0;
# Line 105 | Line 121 | namespace oopse {
121          }
122  
123          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
124 +
125          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
126  
127          //calculate atoms in rigid bodies
# Line 121 | Line 138 | namespace oopse {
138          
139        }
140  
141 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
142 <      //therefore the total number of cutoff groups in the system is equal to
143 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
144 <      //file plus the number of cutoff groups defined in meta-data file
141 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
142 >      //group therefore the total number of cutoff groups in the system is
143 >      //equal to the total number of atoms minus number of atoms belong to
144 >      //cutoff group defined in meta-data file plus the number of cutoff
145 >      //groups defined in meta-data file
146        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
147  
148 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
149 <      //therefore the total number of  integrable objects in the system is equal to
150 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
151 <      //file plus the number of  rigid bodies defined in meta-data file
152 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
153 <
148 >      //every free atom (atom does not belong to rigid bodies) is an
149 >      //integrable object therefore the total number of integrable objects
150 >      //in the system is equal to the total number of atoms minus number of
151 >      //atoms belong to rigid body defined in meta-data file plus the number
152 >      //of rigid bodies defined in meta-data file
153 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
154 >                                                + nGlobalRigidBodies_;
155 >  
156        nGlobalMols_ = molStampIds_.size();
157  
158   #ifdef IS_MPI    
# Line 335 | Line 355 | namespace oopse {
355      int b;
356      int c;
357      int d;
358 +
359 +    std::map<int, std::set<int> > atomGroups;
360 +
361 +    Molecule::RigidBodyIterator rbIter;
362 +    RigidBody* rb;
363 +    Molecule::IntegrableObjectIterator ii;
364 +    StuntDouble* integrableObject;
365      
366 +    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
367 +           integrableObject = mol->nextIntegrableObject(ii)) {
368 +
369 +      if (integrableObject->isRigidBody()) {
370 +          rb = static_cast<RigidBody*>(integrableObject);
371 +          std::vector<Atom*> atoms = rb->getAtoms();
372 +          std::set<int> rigidAtoms;
373 +          for (int i = 0; i < atoms.size(); ++i) {
374 +            rigidAtoms.insert(atoms[i]->getGlobalIndex());
375 +          }
376 +          for (int i = 0; i < atoms.size(); ++i) {
377 +            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
378 +          }      
379 +      } else {
380 +        std::set<int> oneAtomSet;
381 +        oneAtomSet.insert(integrableObject->getGlobalIndex());
382 +        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
383 +      }
384 +    }  
385 +
386 +    
387 +    
388      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
389        a = bond->getAtomA()->getGlobalIndex();
390        b = bond->getAtomB()->getGlobalIndex();        
# Line 346 | Line 395 | namespace oopse {
395        a = bend->getAtomA()->getGlobalIndex();
396        b = bend->getAtomB()->getGlobalIndex();        
397        c = bend->getAtomC()->getGlobalIndex();
398 +      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
399 +      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
400 +      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
401  
402 <      exclude_.addPair(a, b);
403 <      exclude_.addPair(a, c);
404 <      exclude_.addPair(b, c);        
402 >      exclude_.addPairs(rigidSetA, rigidSetB);
403 >      exclude_.addPairs(rigidSetA, rigidSetC);
404 >      exclude_.addPairs(rigidSetB, rigidSetC);
405 >      
406 >      //exclude_.addPair(a, b);
407 >      //exclude_.addPair(a, c);
408 >      //exclude_.addPair(b, c);        
409      }
410  
411      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
# Line 357 | Line 413 | namespace oopse {
413        b = torsion->getAtomB()->getGlobalIndex();        
414        c = torsion->getAtomC()->getGlobalIndex();        
415        d = torsion->getAtomD()->getGlobalIndex();        
416 +      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
417 +      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
418 +      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
419 +      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
420  
421 +      exclude_.addPairs(rigidSetA, rigidSetB);
422 +      exclude_.addPairs(rigidSetA, rigidSetC);
423 +      exclude_.addPairs(rigidSetA, rigidSetD);
424 +      exclude_.addPairs(rigidSetB, rigidSetC);
425 +      exclude_.addPairs(rigidSetB, rigidSetD);
426 +      exclude_.addPairs(rigidSetC, rigidSetD);
427 +
428 +      /*
429 +      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
430 +      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
431 +      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
432 +      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
433 +      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
434 +      exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
435 +        
436 +      
437        exclude_.addPair(a, b);
438        exclude_.addPair(a, c);
439        exclude_.addPair(a, d);
440        exclude_.addPair(b, c);
441        exclude_.addPair(b, d);
442        exclude_.addPair(c, d);        
443 +      */
444      }
445  
369 –    Molecule::RigidBodyIterator rbIter;
370 –    RigidBody* rb;
446      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
447        std::vector<Atom*> atoms = rb->getAtoms();
448        for (int i = 0; i < atoms.size() -1 ; ++i) {
# Line 392 | Line 467 | namespace oopse {
467      int b;
468      int c;
469      int d;
470 +
471 +    std::map<int, std::set<int> > atomGroups;
472 +
473 +    Molecule::RigidBodyIterator rbIter;
474 +    RigidBody* rb;
475 +    Molecule::IntegrableObjectIterator ii;
476 +    StuntDouble* integrableObject;
477      
478 +    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
479 +           integrableObject = mol->nextIntegrableObject(ii)) {
480 +
481 +      if (integrableObject->isRigidBody()) {
482 +          rb = static_cast<RigidBody*>(integrableObject);
483 +          std::vector<Atom*> atoms = rb->getAtoms();
484 +          std::set<int> rigidAtoms;
485 +          for (int i = 0; i < atoms.size(); ++i) {
486 +            rigidAtoms.insert(atoms[i]->getGlobalIndex());
487 +          }
488 +          for (int i = 0; i < atoms.size(); ++i) {
489 +            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
490 +          }      
491 +      } else {
492 +        std::set<int> oneAtomSet;
493 +        oneAtomSet.insert(integrableObject->getGlobalIndex());
494 +        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
495 +      }
496 +    }  
497 +
498 +    
499      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
500        a = bond->getAtomA()->getGlobalIndex();
501        b = bond->getAtomB()->getGlobalIndex();        
# Line 404 | Line 507 | namespace oopse {
507        b = bend->getAtomB()->getGlobalIndex();        
508        c = bend->getAtomC()->getGlobalIndex();
509  
510 <      exclude_.removePair(a, b);
511 <      exclude_.removePair(a, c);
512 <      exclude_.removePair(b, c);        
510 >      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
511 >      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
512 >      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
513 >
514 >      exclude_.removePairs(rigidSetA, rigidSetB);
515 >      exclude_.removePairs(rigidSetA, rigidSetC);
516 >      exclude_.removePairs(rigidSetB, rigidSetC);
517 >      
518 >      //exclude_.removePair(a, b);
519 >      //exclude_.removePair(a, c);
520 >      //exclude_.removePair(b, c);        
521      }
522  
523      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
# Line 415 | Line 526 | namespace oopse {
526        c = torsion->getAtomC()->getGlobalIndex();        
527        d = torsion->getAtomD()->getGlobalIndex();        
528  
529 +      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
530 +      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
531 +      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
532 +      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
533 +
534 +      exclude_.removePairs(rigidSetA, rigidSetB);
535 +      exclude_.removePairs(rigidSetA, rigidSetC);
536 +      exclude_.removePairs(rigidSetA, rigidSetD);
537 +      exclude_.removePairs(rigidSetB, rigidSetC);
538 +      exclude_.removePairs(rigidSetB, rigidSetD);
539 +      exclude_.removePairs(rigidSetC, rigidSetD);
540 +
541 +      /*
542 +      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
543 +      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
544 +      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
545 +      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
546 +      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
547 +      exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
548 +
549 +      
550        exclude_.removePair(a, b);
551        exclude_.removePair(a, c);
552        exclude_.removePair(a, d);
553        exclude_.removePair(b, c);
554        exclude_.removePair(b, d);
555        exclude_.removePair(c, d);        
556 +      */
557      }
558  
426 –    Molecule::RigidBodyIterator rbIter;
427 –    RigidBody* rb;
559      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
560        std::vector<Atom*> atoms = rb->getAtoms();
561        for (int i = 0; i < atoms.size() -1 ; ++i) {
# Line 462 | Line 593 | namespace oopse {
593      //setup fortran force field
594      /** @deprecate */    
595      int isError = 0;
596 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
596 >    
597 >    setupElectrostaticSummationMethod( isError );
598 >    setupSwitchingFunction();
599 >
600      if(isError){
601        sprintf( painCave.errMsg,
602                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 506 | Line 640 | namespace oopse {
640      int useLennardJones = 0;
641      int useElectrostatic = 0;
642      int useEAM = 0;
643 +    int useSC = 0;
644      int useCharge = 0;
645      int useDirectional = 0;
646      int useDipole = 0;
647      int useGayBerne = 0;
648      int useSticky = 0;
649 +    int useStickyPower = 0;
650      int useShape = 0;
651      int useFLARB = 0; //it is not in AtomType yet
652      int useDirectionalAtom = 0;    
653      int useElectrostatics = 0;
654      //usePBC and useRF are from simParams
655 <    int usePBC = simParams_->getPBC();
656 <    int useRF = simParams_->getUseRF();
655 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
656 >    int useRF;
657 >    int useSF;
658 >    std::string myMethod;
659 >
660 >    // set the useRF logical
661 >    useRF = 0;
662 >    useSF = 0;
663 >
664 >
665 >    if (simParams_->haveElectrostaticSummationMethod()) {
666 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
667 >      toUpper(myMethod);
668 >      if (myMethod == "REACTION_FIELD") {
669 >        useRF=1;
670 >      } else {
671 >        if (myMethod == "SHIFTED_FORCE") {
672 >          useSF = 1;
673 >        }
674 >      }
675 >    }
676  
677      //loop over all of the atom types
678      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
679        useLennardJones |= (*i)->isLennardJones();
680        useElectrostatic |= (*i)->isElectrostatic();
681        useEAM |= (*i)->isEAM();
682 +      useSC |= (*i)->isSC();
683        useCharge |= (*i)->isCharge();
684        useDirectional |= (*i)->isDirectional();
685        useDipole |= (*i)->isDipole();
686        useGayBerne |= (*i)->isGayBerne();
687        useSticky |= (*i)->isSticky();
688 +      useStickyPower |= (*i)->isStickyPower();
689        useShape |= (*i)->isShape();
690      }
691  
692 <    if (useSticky || useDipole || useGayBerne || useShape) {
692 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
693        useDirectionalAtom = 1;
694      }
695  
# Line 564 | Line 721 | namespace oopse {
721      temp = useSticky;
722      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
723  
724 +    temp = useStickyPower;
725 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
726 +    
727      temp = useGayBerne;
728      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
729  
730      temp = useEAM;
731      MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
732  
733 +    temp = useSC;
734 +    MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
735 +    
736      temp = useShape;
737      MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
738  
# Line 578 | Line 741 | namespace oopse {
741  
742      temp = useRF;
743      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
744 <    
744 >
745 >    temp = useSF;
746 >    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
747 >
748   #endif
749  
750      fInfo_.SIM_uses_PBC = usePBC;    
# Line 588 | Line 754 | namespace oopse {
754      fInfo_.SIM_uses_Charges = useCharge;
755      fInfo_.SIM_uses_Dipoles = useDipole;
756      fInfo_.SIM_uses_Sticky = useSticky;
757 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
758      fInfo_.SIM_uses_GayBerne = useGayBerne;
759      fInfo_.SIM_uses_EAM = useEAM;
760 +    fInfo_.SIM_uses_SC = useSC;
761      fInfo_.SIM_uses_Shapes = useShape;
762      fInfo_.SIM_uses_FLARB = useFLARB;
763      fInfo_.SIM_uses_RF = useRF;
764 +    fInfo_.SIM_uses_SF = useSF;
765  
766 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
767 <
766 >    if( myMethod == "REACTION_FIELD") {
767 >      
768        if (simParams_->haveDielectric()) {
769          fInfo_.dielect = simParams_->getDielectric();
770        } else {
# Line 605 | Line 774 | namespace oopse {
774                  "\tsetting a dielectric constant!\n");
775          painCave.isFatal = 1;
776          simError();
777 <      }
609 <        
610 <    } else {
611 <      fInfo_.dielect = 0.0;
777 >      }      
778      }
779  
780    }
# Line 644 | Line 810 | namespace oopse {
810  
811          totalMass = cg->getMass();
812          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
813 <          mfact.push_back(atom->getMass()/totalMass);
813 >          // Check for massless groups - set mfact to 1 if true
814 >          if (totalMass != 0)
815 >            mfact.push_back(atom->getMass()/totalMass);
816 >          else
817 >            mfact.push_back( 1.0 );
818          }
819  
820        }      
# Line 780 | Line 950 | namespace oopse {
950      
951      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
952          
953 <      if (!simParams_->haveRcut()){
953 >      if (!simParams_->haveCutoffRadius()){
954          sprintf(painCave.errMsg,
955                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
956                  "\tOOPSE will use a default value of 15.0 angstroms"
# Line 789 | Line 959 | namespace oopse {
959          simError();
960          rcut = 15.0;
961        } else{
962 <        rcut = simParams_->getRcut();
962 >        rcut = simParams_->getCutoffRadius();
963        }
964  
965 <      if (!simParams_->haveRsw()){
965 >      if (!simParams_->haveSwitchingRadius()){
966          sprintf(painCave.errMsg,
967                  "SimCreator Warning: No value was set for switchingRadius.\n"
968                  "\tOOPSE will use a default value of\n"
969 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
969 >                "\t0.85 * cutoffRadius for the switchingRadius\n");
970          painCave.isFatal = 0;
971          simError();
972 <        rsw = 0.95 * rcut;
972 >        rsw = 0.85 * rcut;
973        } else{
974 <        rsw = simParams_->getRsw();
974 >        rsw = simParams_->getSwitchingRadius();
975        }
976  
977      } else {
978        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
979        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
980          
981 <      if (simParams_->haveRcut()) {
982 <        rcut = simParams_->getRcut();
981 >      if (simParams_->haveCutoffRadius()) {
982 >        rcut = simParams_->getCutoffRadius();
983        } else {
984          //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
985          rcut = calcMaxCutoffRadius();
986        }
987  
988 <      if (simParams_->haveRsw()) {
989 <        rsw  = simParams_->getRsw();
988 >      if (simParams_->haveSwitchingRadius()) {
989 >        rsw  = simParams_->getSwitchingRadius();
990        } else {
991          rsw = rcut;
992        }
# Line 824 | Line 994 | namespace oopse {
994      }
995    }
996  
997 <  void SimInfo::setupCutoff() {
997 >  void SimInfo::setupCutoff() {    
998      getCutoff(rcut_, rsw_);    
999      double rnblist = rcut_ + 1; // skin of neighbor list
1000  
1001      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
1002 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
1002 >    
1003 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
1004 >    if (simParams_->haveCutoffPolicy()) {
1005 >      std::string myPolicy = simParams_->getCutoffPolicy();
1006 >      toUpper(myPolicy);
1007 >      if (myPolicy == "MIX") {
1008 >        cp = MIX_CUTOFF_POLICY;
1009 >      } else {
1010 >        if (myPolicy == "MAX") {
1011 >          cp = MAX_CUTOFF_POLICY;
1012 >        } else {
1013 >          if (myPolicy == "TRADITIONAL") {            
1014 >            cp = TRADITIONAL_CUTOFF_POLICY;
1015 >          } else {
1016 >            // throw error        
1017 >            sprintf( painCave.errMsg,
1018 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
1019 >            painCave.isFatal = 1;
1020 >            simError();
1021 >          }    
1022 >        }          
1023 >      }
1024 >    }
1025 >
1026 >
1027 >    if (simParams_->haveSkinThickness()) {
1028 >      double skinThickness = simParams_->getSkinThickness();
1029 >    }
1030 >
1031 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
1032 >    // also send cutoff notification to electrostatics
1033 >    setElectrostaticCutoffRadius(&rcut_, &rsw_);
1034    }
1035  
1036 +  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
1037 +    
1038 +    int errorOut;
1039 +    int esm =  NONE;
1040 +    int sm = UNDAMPED;
1041 +    double alphaVal;
1042 +    double dielectric;
1043 +
1044 +    errorOut = isError;
1045 +    alphaVal = simParams_->getDampingAlpha();
1046 +    dielectric = simParams_->getDielectric();
1047 +
1048 +    if (simParams_->haveElectrostaticSummationMethod()) {
1049 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
1050 +      toUpper(myMethod);
1051 +      if (myMethod == "NONE") {
1052 +        esm = NONE;
1053 +      } else {
1054 +        if (myMethod == "SWITCHING_FUNCTION") {
1055 +          esm = SWITCHING_FUNCTION;
1056 +        } else {
1057 +          if (myMethod == "SHIFTED_POTENTIAL") {
1058 +            esm = SHIFTED_POTENTIAL;
1059 +          } else {
1060 +            if (myMethod == "SHIFTED_FORCE") {            
1061 +              esm = SHIFTED_FORCE;
1062 +            } else {
1063 +              if (myMethod == "REACTION_FIELD") {            
1064 +                esm = REACTION_FIELD;
1065 +              } else {
1066 +                // throw error        
1067 +                sprintf( painCave.errMsg,
1068 +                         "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() );
1069 +                painCave.isFatal = 1;
1070 +                simError();
1071 +              }    
1072 +            }          
1073 +          }
1074 +        }
1075 +      }
1076 +    }
1077 +    
1078 +    if (simParams_->haveElectrostaticScreeningMethod()) {
1079 +      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
1080 +      toUpper(myScreen);
1081 +      if (myScreen == "UNDAMPED") {
1082 +        sm = UNDAMPED;
1083 +      } else {
1084 +        if (myScreen == "DAMPED") {
1085 +          sm = DAMPED;
1086 +          if (!simParams_->haveDampingAlpha()) {
1087 +            //throw error
1088 +            sprintf( painCave.errMsg,
1089 +                     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal);
1090 +            painCave.isFatal = 0;
1091 +            simError();
1092 +          }
1093 +        } else {
1094 +          // throw error        
1095 +          sprintf( painCave.errMsg,
1096 +                   "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() );
1097 +          painCave.isFatal = 1;
1098 +          simError();
1099 +        }
1100 +      }
1101 +    }
1102 +    
1103 +    // let's pass some summation method variables to fortran
1104 +    setElectrostaticSummationMethod( &esm );
1105 +    setScreeningMethod( &sm );
1106 +    setDampingAlpha( &alphaVal );
1107 +    setReactionFieldDielectric( &dielectric );
1108 +    initFortranFF( &esm, &errorOut );
1109 +  }
1110 +
1111 +  void SimInfo::setupSwitchingFunction() {    
1112 +    int ft = CUBIC;
1113 +
1114 +    if (simParams_->haveSwitchingFunctionType()) {
1115 +      std::string funcType = simParams_->getSwitchingFunctionType();
1116 +      toUpper(funcType);
1117 +      if (funcType == "CUBIC") {
1118 +        ft = CUBIC;
1119 +      } else {
1120 +        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1121 +          ft = FIFTH_ORDER_POLY;
1122 +        } else {
1123 +          // throw error        
1124 +          sprintf( painCave.errMsg,
1125 +                   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1126 +          painCave.isFatal = 1;
1127 +          simError();
1128 +        }          
1129 +      }
1130 +    }
1131 +
1132 +    // send switching function notification to switcheroo
1133 +    setFunctionType(&ft);
1134 +
1135 +  }
1136 +
1137    void SimInfo::addProperty(GenericData* genData) {
1138      properties_.addProperty(genData);  
1139    }
# Line 939 | Line 1241 | namespace oopse {
1241  
1242      return o;
1243    }
1244 +  
1245 +  
1246 +   /*
1247 +   Returns center of mass and center of mass velocity in one function call.
1248 +   */
1249 +  
1250 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1251 +      SimInfo::MoleculeIterator i;
1252 +      Molecule* mol;
1253 +      
1254 +    
1255 +      double totalMass = 0.0;
1256 +    
1257  
1258 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1259 +         double mass = mol->getMass();
1260 +         totalMass += mass;
1261 +         com += mass * mol->getCom();
1262 +         comVel += mass * mol->getComVel();          
1263 +      }  
1264 +      
1265 + #ifdef IS_MPI
1266 +      double tmpMass = totalMass;
1267 +      Vector3d tmpCom(com);  
1268 +      Vector3d tmpComVel(comVel);
1269 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1270 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1271 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1272 + #endif
1273 +      
1274 +      com /= totalMass;
1275 +      comVel /= totalMass;
1276 +   }        
1277 +  
1278 +   /*
1279 +   Return intertia tensor for entire system and angular momentum Vector.
1280 +
1281 +
1282 +       [  Ixx -Ixy  -Ixz ]
1283 +  J =| -Iyx  Iyy  -Iyz |
1284 +       [ -Izx -Iyz   Izz ]
1285 +    */
1286 +
1287 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1288 +      
1289 +
1290 +      double xx = 0.0;
1291 +      double yy = 0.0;
1292 +      double zz = 0.0;
1293 +      double xy = 0.0;
1294 +      double xz = 0.0;
1295 +      double yz = 0.0;
1296 +      Vector3d com(0.0);
1297 +      Vector3d comVel(0.0);
1298 +      
1299 +      getComAll(com, comVel);
1300 +      
1301 +      SimInfo::MoleculeIterator i;
1302 +      Molecule* mol;
1303 +      
1304 +      Vector3d thisq(0.0);
1305 +      Vector3d thisv(0.0);
1306 +
1307 +      double thisMass = 0.0;
1308 +    
1309 +      
1310 +      
1311 +  
1312 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1313 +        
1314 +         thisq = mol->getCom()-com;
1315 +         thisv = mol->getComVel()-comVel;
1316 +         thisMass = mol->getMass();
1317 +         // Compute moment of intertia coefficients.
1318 +         xx += thisq[0]*thisq[0]*thisMass;
1319 +         yy += thisq[1]*thisq[1]*thisMass;
1320 +         zz += thisq[2]*thisq[2]*thisMass;
1321 +        
1322 +         // compute products of intertia
1323 +         xy += thisq[0]*thisq[1]*thisMass;
1324 +         xz += thisq[0]*thisq[2]*thisMass;
1325 +         yz += thisq[1]*thisq[2]*thisMass;
1326 +            
1327 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1328 +            
1329 +      }  
1330 +      
1331 +      
1332 +      inertiaTensor(0,0) = yy + zz;
1333 +      inertiaTensor(0,1) = -xy;
1334 +      inertiaTensor(0,2) = -xz;
1335 +      inertiaTensor(1,0) = -xy;
1336 +      inertiaTensor(1,1) = xx + zz;
1337 +      inertiaTensor(1,2) = -yz;
1338 +      inertiaTensor(2,0) = -xz;
1339 +      inertiaTensor(2,1) = -yz;
1340 +      inertiaTensor(2,2) = xx + yy;
1341 +      
1342 + #ifdef IS_MPI
1343 +      Mat3x3d tmpI(inertiaTensor);
1344 +      Vector3d tmpAngMom;
1345 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1346 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1347 + #endif
1348 +              
1349 +      return;
1350 +   }
1351 +
1352 +   //Returns the angular momentum of the system
1353 +   Vector3d SimInfo::getAngularMomentum(){
1354 +      
1355 +      Vector3d com(0.0);
1356 +      Vector3d comVel(0.0);
1357 +      Vector3d angularMomentum(0.0);
1358 +      
1359 +      getComAll(com,comVel);
1360 +      
1361 +      SimInfo::MoleculeIterator i;
1362 +      Molecule* mol;
1363 +      
1364 +      Vector3d thisr(0.0);
1365 +      Vector3d thisp(0.0);
1366 +      
1367 +      double thisMass;
1368 +      
1369 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1370 +        thisMass = mol->getMass();
1371 +        thisr = mol->getCom()-com;
1372 +        thisp = (mol->getComVel()-comVel)*thisMass;
1373 +        
1374 +        angularMomentum += cross( thisr, thisp );
1375 +        
1376 +      }  
1377 +      
1378 + #ifdef IS_MPI
1379 +      Vector3d tmpAngMom;
1380 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1381 + #endif
1382 +      
1383 +      return angularMomentum;
1384 +   }
1385 +  
1386 +  
1387   }//end namespace oopse
1388  

Diff Legend

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