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 598 by chrisfen, Thu Sep 15 00:14:35 2005 UTC vs.
Revision 1024 by tim, Wed Aug 30 18:42:29 2006 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 "primitives/StuntDouble.hpp"
57   #include "UseTheForce/fCutoffPolicy.h"
58 < #include "UseTheForce/fCoulombicCorrection.h"
58 > #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
59 > #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
60 > #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
61   #include "UseTheForce/doForces_interface.h"
62 < #include "UseTheForce/notifyCutoffs_interface.h"
62 > #include "UseTheForce/DarkSide/electrostatic_interface.h"
63 > #include "UseTheForce/DarkSide/switcheroo_interface.h"
64   #include "utils/MemoryUtils.hpp"
65   #include "utils/simError.h"
66   #include "selection/SelectionManager.hpp"
67 + #include "io/ForceFieldOptions.hpp"
68 + #include "UseTheForce/ForceField.hpp"
69  
70   #ifdef IS_MPI
71   #include "UseTheForce/mpiComponentPlan.h"
# Line 66 | Line 73 | namespace oopse {
73   #endif
74  
75   namespace oopse {
76 +  std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) {
77 +    std::map<int, std::set<int> >::iterator i = container.find(index);
78 +    std::set<int> result;
79 +    if (i != container.end()) {
80 +        result = i->second;
81 +    }
82  
83 <  SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
84 <                   ForceField* ff, Globals* simParams) :
85 <    stamps_(stamps), forceField_(ff), simParams_(simParams),
86 <    ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
83 >    return result;
84 >  }
85 >  
86 >  SimInfo::SimInfo(ForceField* ff, Globals* simParams) :
87 >    forceField_(ff), simParams_(simParams),
88 >    ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
89      nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
90      nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
91      nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
92      nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
93 <    sman_(NULL), fortranInitialized_(false) {
93 >    sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) {
94  
80            
81      std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
95        MoleculeStamp* molStamp;
96        int nMolWithSameStamp;
97        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
98 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
98 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
99        CutoffGroupStamp* cgStamp;    
100        RigidBodyStamp* rbStamp;
101        int nRigidAtoms = 0;
102 <    
103 <      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
104 <        molStamp = i->first;
105 <        nMolWithSameStamp = i->second;
102 >      std::vector<Component*> components = simParams->getComponents();
103 >      
104 >      for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
105 >        molStamp = (*i)->getMoleculeStamp();
106 >        nMolWithSameStamp = (*i)->getNMol();
107          
108          addMoleculeStamp(molStamp, nMolWithSameStamp);
109  
110          //calculate atoms in molecules
111          nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
98
112  
113          //calculate atoms in cutoff groups
114          int nAtomsInGroups = 0;
115          int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
116          
117          for (int j=0; j < nCutoffGroupsInStamp; j++) {
118 <          cgStamp = molStamp->getCutoffGroup(j);
118 >          cgStamp = molStamp->getCutoffGroupStamp(j);
119            nAtomsInGroups += cgStamp->getNMembers();
120          }
121  
122          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
123 +
124          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
125  
126          //calculate atoms in rigid bodies
# Line 114 | Line 128 | namespace oopse {
128          int nRigidBodiesInStamp = molStamp->getNRigidBodies();
129          
130          for (int j=0; j < nRigidBodiesInStamp; j++) {
131 <          rbStamp = molStamp->getRigidBody(j);
131 >          rbStamp = molStamp->getRigidBodyStamp(j);
132            nAtomsInRigidBodies += rbStamp->getNMembers();
133          }
134  
# Line 123 | Line 137 | namespace oopse {
137          
138        }
139  
140 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
141 <      //therefore the total number of cutoff groups in the system is equal to
142 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
143 <      //file plus the number of cutoff groups defined in meta-data file
140 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
141 >      //group therefore the total number of cutoff groups in the system is
142 >      //equal to the total number of atoms minus number of atoms belong to
143 >      //cutoff group defined in meta-data file plus the number of cutoff
144 >      //groups defined in meta-data file
145        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
146  
147 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
148 <      //therefore the total number of  integrable objects in the system is equal to
149 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
150 <      //file plus the number of  rigid bodies defined in meta-data file
151 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
152 <
147 >      //every free atom (atom does not belong to rigid bodies) is an
148 >      //integrable object therefore the total number of integrable objects
149 >      //in the system is equal to the total number of atoms minus number of
150 >      //atoms belong to rigid body defined in meta-data file plus the number
151 >      //of rigid bodies defined in meta-data file
152 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
153 >                                                + nGlobalRigidBodies_;
154 >  
155        nGlobalMols_ = molStampIds_.size();
156  
157   #ifdef IS_MPI    
# Line 150 | Line 167 | namespace oopse {
167      }
168      molecules_.clear();
169        
153    delete stamps_;
170      delete sman_;
171      delete simParams_;
172      delete forceField_;
# Line 257 | Line 273 | namespace oopse {
273            }
274          }
275              
276 <      }//end for (integrableObject)
277 <    }// end for (mol)
276 >      }
277 >    }
278      
279      // n_constraints is local, so subtract them on each processor
280      ndf_local -= nConstraints_;
# Line 275 | Line 291 | namespace oopse {
291  
292    }
293  
294 +  int SimInfo::getFdf() {
295 + #ifdef IS_MPI
296 +    MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
297 + #else
298 +    fdf_ = fdf_local;
299 + #endif
300 +    return fdf_;
301 +  }
302 +    
303    void SimInfo::calcNdfRaw() {
304      int ndfRaw_local;
305  
# Line 337 | Line 362 | namespace oopse {
362      int b;
363      int c;
364      int d;
365 +
366 +    std::map<int, std::set<int> > atomGroups;
367 +
368 +    Molecule::RigidBodyIterator rbIter;
369 +    RigidBody* rb;
370 +    Molecule::IntegrableObjectIterator ii;
371 +    StuntDouble* integrableObject;
372      
373 +    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
374 +           integrableObject = mol->nextIntegrableObject(ii)) {
375 +
376 +      if (integrableObject->isRigidBody()) {
377 +          rb = static_cast<RigidBody*>(integrableObject);
378 +          std::vector<Atom*> atoms = rb->getAtoms();
379 +          std::set<int> rigidAtoms;
380 +          for (int i = 0; i < atoms.size(); ++i) {
381 +            rigidAtoms.insert(atoms[i]->getGlobalIndex());
382 +          }
383 +          for (int i = 0; i < atoms.size(); ++i) {
384 +            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
385 +          }      
386 +      } else {
387 +        std::set<int> oneAtomSet;
388 +        oneAtomSet.insert(integrableObject->getGlobalIndex());
389 +        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
390 +      }
391 +    }  
392 +
393 +    
394 +    
395      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
396        a = bond->getAtomA()->getGlobalIndex();
397        b = bond->getAtomB()->getGlobalIndex();        
# Line 348 | Line 402 | namespace oopse {
402        a = bend->getAtomA()->getGlobalIndex();
403        b = bend->getAtomB()->getGlobalIndex();        
404        c = bend->getAtomC()->getGlobalIndex();
405 +      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
406 +      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
407 +      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
408  
409 <      exclude_.addPair(a, b);
410 <      exclude_.addPair(a, c);
411 <      exclude_.addPair(b, c);        
409 >      exclude_.addPairs(rigidSetA, rigidSetB);
410 >      exclude_.addPairs(rigidSetA, rigidSetC);
411 >      exclude_.addPairs(rigidSetB, rigidSetC);
412 >      
413 >      //exclude_.addPair(a, b);
414 >      //exclude_.addPair(a, c);
415 >      //exclude_.addPair(b, c);        
416      }
417  
418      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
# Line 359 | Line 420 | namespace oopse {
420        b = torsion->getAtomB()->getGlobalIndex();        
421        c = torsion->getAtomC()->getGlobalIndex();        
422        d = torsion->getAtomD()->getGlobalIndex();        
423 +      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
424 +      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
425 +      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
426 +      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
427  
428 +      exclude_.addPairs(rigidSetA, rigidSetB);
429 +      exclude_.addPairs(rigidSetA, rigidSetC);
430 +      exclude_.addPairs(rigidSetA, rigidSetD);
431 +      exclude_.addPairs(rigidSetB, rigidSetC);
432 +      exclude_.addPairs(rigidSetB, rigidSetD);
433 +      exclude_.addPairs(rigidSetC, rigidSetD);
434 +
435 +      /*
436 +      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
437 +      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
438 +      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
439 +      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
440 +      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
441 +      exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
442 +        
443 +      
444        exclude_.addPair(a, b);
445        exclude_.addPair(a, c);
446        exclude_.addPair(a, d);
447        exclude_.addPair(b, c);
448        exclude_.addPair(b, d);
449        exclude_.addPair(c, d);        
450 +      */
451      }
452  
371    Molecule::RigidBodyIterator rbIter;
372    RigidBody* rb;
453      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
454        std::vector<Atom*> atoms = rb->getAtoms();
455        for (int i = 0; i < atoms.size() -1 ; ++i) {
# Line 394 | Line 474 | namespace oopse {
474      int b;
475      int c;
476      int d;
477 +
478 +    std::map<int, std::set<int> > atomGroups;
479 +
480 +    Molecule::RigidBodyIterator rbIter;
481 +    RigidBody* rb;
482 +    Molecule::IntegrableObjectIterator ii;
483 +    StuntDouble* integrableObject;
484      
485 +    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
486 +           integrableObject = mol->nextIntegrableObject(ii)) {
487 +
488 +      if (integrableObject->isRigidBody()) {
489 +          rb = static_cast<RigidBody*>(integrableObject);
490 +          std::vector<Atom*> atoms = rb->getAtoms();
491 +          std::set<int> rigidAtoms;
492 +          for (int i = 0; i < atoms.size(); ++i) {
493 +            rigidAtoms.insert(atoms[i]->getGlobalIndex());
494 +          }
495 +          for (int i = 0; i < atoms.size(); ++i) {
496 +            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
497 +          }      
498 +      } else {
499 +        std::set<int> oneAtomSet;
500 +        oneAtomSet.insert(integrableObject->getGlobalIndex());
501 +        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
502 +      }
503 +    }  
504 +
505 +    
506      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
507        a = bond->getAtomA()->getGlobalIndex();
508        b = bond->getAtomB()->getGlobalIndex();        
# Line 406 | Line 514 | namespace oopse {
514        b = bend->getAtomB()->getGlobalIndex();        
515        c = bend->getAtomC()->getGlobalIndex();
516  
517 <      exclude_.removePair(a, b);
518 <      exclude_.removePair(a, c);
519 <      exclude_.removePair(b, c);        
517 >      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
518 >      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
519 >      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
520 >
521 >      exclude_.removePairs(rigidSetA, rigidSetB);
522 >      exclude_.removePairs(rigidSetA, rigidSetC);
523 >      exclude_.removePairs(rigidSetB, rigidSetC);
524 >      
525 >      //exclude_.removePair(a, b);
526 >      //exclude_.removePair(a, c);
527 >      //exclude_.removePair(b, c);        
528      }
529  
530      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
# Line 417 | Line 533 | namespace oopse {
533        c = torsion->getAtomC()->getGlobalIndex();        
534        d = torsion->getAtomD()->getGlobalIndex();        
535  
536 +      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
537 +      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
538 +      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
539 +      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
540 +
541 +      exclude_.removePairs(rigidSetA, rigidSetB);
542 +      exclude_.removePairs(rigidSetA, rigidSetC);
543 +      exclude_.removePairs(rigidSetA, rigidSetD);
544 +      exclude_.removePairs(rigidSetB, rigidSetC);
545 +      exclude_.removePairs(rigidSetB, rigidSetD);
546 +      exclude_.removePairs(rigidSetC, rigidSetD);
547 +
548 +      /*
549 +      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
550 +      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
551 +      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
552 +      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
553 +      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
554 +      exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
555 +
556 +      
557        exclude_.removePair(a, b);
558        exclude_.removePair(a, c);
559        exclude_.removePair(a, d);
560        exclude_.removePair(b, c);
561        exclude_.removePair(b, d);
562        exclude_.removePair(c, d);        
563 +      */
564      }
565  
428    Molecule::RigidBodyIterator rbIter;
429    RigidBody* rb;
566      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
567        std::vector<Atom*> atoms = rb->getAtoms();
568        for (int i = 0; i < atoms.size() -1 ; ++i) {
# Line 465 | Line 601 | namespace oopse {
601      /** @deprecate */    
602      int isError = 0;
603      
604 <    setupCoulombicCorrection( isError );
604 >    setupElectrostaticSummationMethod( isError );
605 >    setupSwitchingFunction();
606 >    setupAccumulateBoxDipole();
607  
608      if(isError){
609        sprintf( painCave.errMsg,
# Line 510 | Line 648 | namespace oopse {
648      int useLennardJones = 0;
649      int useElectrostatic = 0;
650      int useEAM = 0;
651 +    int useSC = 0;
652      int useCharge = 0;
653      int useDirectional = 0;
654      int useDipole = 0;
# Line 521 | Line 660 | namespace oopse {
660      int useDirectionalAtom = 0;    
661      int useElectrostatics = 0;
662      //usePBC and useRF are from simParams
663 <    int usePBC = simParams_->getPBC();
664 <    int useRF = simParams_->getUseRF();
663 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
664 >    int useRF;
665 >    int useSF;
666 >    int useSP;
667 >    int useBoxDipole;
668 >    std::string myMethod;
669  
670 +    // set the useRF logical
671 +    useRF = 0;
672 +    useSF = 0;
673 +
674 +
675 +    if (simParams_->haveElectrostaticSummationMethod()) {
676 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
677 +      toUpper(myMethod);
678 +      if (myMethod == "REACTION_FIELD"){
679 +        useRF=1;
680 +      } else if (myMethod == "SHIFTED_FORCE"){
681 +        useSF = 1;
682 +      } else if (myMethod == "SHIFTED_POTENTIAL"){
683 +        useSP = 1;
684 +      }
685 +    }
686 +    
687 +    if (simParams_->haveAccumulateBoxDipole())
688 +      if (simParams_->getAccumulateBoxDipole())
689 +        useBoxDipole = 1;
690 +
691      //loop over all of the atom types
692      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
693        useLennardJones |= (*i)->isLennardJones();
694        useElectrostatic |= (*i)->isElectrostatic();
695        useEAM |= (*i)->isEAM();
696 +      useSC |= (*i)->isSC();
697        useCharge |= (*i)->isCharge();
698        useDirectional |= (*i)->isDirectional();
699        useDipole |= (*i)->isDipole();
# Line 579 | Line 744 | namespace oopse {
744      temp = useEAM;
745      MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
746  
747 +    temp = useSC;
748 +    MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
749 +    
750      temp = useShape;
751      MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
752  
# Line 588 | Line 756 | namespace oopse {
756      temp = useRF;
757      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
758  
759 <    temp = useUW;
760 <    MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
759 >    temp = useSF;
760 >    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
761  
762 <    temp = useDW;
763 <    MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
764 <    
762 >    temp = useSP;
763 >    MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
764 >
765 >    temp = useBoxDipole;
766 >    MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
767 >
768   #endif
769  
770      fInfo_.SIM_uses_PBC = usePBC;    
# Line 606 | Line 777 | namespace oopse {
777      fInfo_.SIM_uses_StickyPower = useStickyPower;
778      fInfo_.SIM_uses_GayBerne = useGayBerne;
779      fInfo_.SIM_uses_EAM = useEAM;
780 +    fInfo_.SIM_uses_SC = useSC;
781      fInfo_.SIM_uses_Shapes = useShape;
782      fInfo_.SIM_uses_FLARB = useFLARB;
783      fInfo_.SIM_uses_RF = useRF;
784 +    fInfo_.SIM_uses_SF = useSF;
785 +    fInfo_.SIM_uses_SP = useSP;
786 +    fInfo_.SIM_uses_BoxDipole = useBoxDipole;
787  
788 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
789 <
788 >    if( myMethod == "REACTION_FIELD") {
789 >      
790        if (simParams_->haveDielectric()) {
791          fInfo_.dielect = simParams_->getDielectric();
792        } else {
# Line 621 | Line 796 | namespace oopse {
796                  "\tsetting a dielectric constant!\n");
797          painCave.isFatal = 1;
798          simError();
799 <      }
625 <        
626 <    } else {
627 <      fInfo_.dielect = 0.0;
799 >      }      
800      }
801  
802    }
# Line 643 | Line 815 | namespace oopse {
815      }
816  
817      //calculate mass ratio of cutoff group
818 <    std::vector<double> mfact;
818 >    std::vector<RealType> mfact;
819      SimInfo::MoleculeIterator mi;
820      Molecule* mol;
821      Molecule::CutoffGroupIterator ci;
822      CutoffGroup* cg;
823      Molecule::AtomIterator ai;
824      Atom* atom;
825 <    double totalMass;
825 >    RealType totalMass;
826  
827      //to avoid memory reallocation, reserve enough space for mfact
828      mfact.reserve(getNCutoffGroups());
# Line 660 | Line 832 | namespace oopse {
832  
833          totalMass = cg->getMass();
834          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
835 <          mfact.push_back(atom->getMass()/totalMass);
835 >          // Check for massless groups - set mfact to 1 if true
836 >          if (totalMass != 0)
837 >            mfact.push_back(atom->getMass()/totalMass);
838 >          else
839 >            mfact.push_back( 1.0 );
840          }
841  
842        }      
# Line 769 | Line 945 | namespace oopse {
945  
946   #endif
947  
948 <  double SimInfo::calcMaxCutoffRadius() {
948 >  void SimInfo::setupCutoff() {          
949 >    
950 >    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
951  
952 +    // Check the cutoff policy
953 +    int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default
954  
955 <    std::set<AtomType*> atomTypes;
956 <    std::set<AtomType*>::iterator i;
957 <    std::vector<double> cutoffRadius;
958 <
959 <    //get the unique atom types
780 <    atomTypes = getUniqueAtomTypes();
781 <
782 <    //query the max cutoff radius among these atom types
783 <    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
784 <      cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
955 >    std::string myPolicy;
956 >    if (forceFieldOptions_.haveCutoffPolicy()){
957 >      myPolicy = forceFieldOptions_.getCutoffPolicy();
958 >    }else if (simParams_->haveCutoffPolicy()) {
959 >      myPolicy = simParams_->getCutoffPolicy();
960      }
961  
962 <    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
963 < #ifdef IS_MPI
789 <    //pick the max cutoff radius among the processors
790 < #endif
791 <
792 <    return maxCutoffRadius;
793 <  }
794 <
795 <  void SimInfo::getCutoff(double& rcut, double& rsw) {
796 <    
797 <    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
798 <        
799 <      if (!simParams_->haveRcut()){
800 <        sprintf(painCave.errMsg,
801 <                "SimCreator Warning: No value was set for the cutoffRadius.\n"
802 <                "\tOOPSE will use a default value of 15.0 angstroms"
803 <                "\tfor the cutoffRadius.\n");
804 <        painCave.isFatal = 0;
805 <        simError();
806 <        rcut = 15.0;
807 <      } else{
808 <        rcut = simParams_->getRcut();
809 <      }
810 <
811 <      if (!simParams_->haveRsw()){
812 <        sprintf(painCave.errMsg,
813 <                "SimCreator Warning: No value was set for switchingRadius.\n"
814 <                "\tOOPSE will use a default value of\n"
815 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
816 <        painCave.isFatal = 0;
817 <        simError();
818 <        rsw = 0.95 * rcut;
819 <      } else{
820 <        rsw = simParams_->getRsw();
821 <      }
822 <
823 <    } else {
824 <      // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
825 <      //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
826 <        
827 <      if (simParams_->haveRcut()) {
828 <        rcut = simParams_->getRcut();
829 <      } else {
830 <        //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
831 <        rcut = calcMaxCutoffRadius();
832 <      }
833 <
834 <      if (simParams_->haveRsw()) {
835 <        rsw  = simParams_->getRsw();
836 <      } else {
837 <        rsw = rcut;
838 <      }
839 <    
840 <    }
841 <  }
842 <
843 <  void SimInfo::setupCutoff() {    
844 <    getCutoff(rcut_, rsw_);    
845 <    double rnblist = rcut_ + 1; // skin of neighbor list
846 <
847 <    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
848 <    
849 <    int cp =  TRADITIONAL_CUTOFF_POLICY;
850 <    if (simParams_->haveCutoffPolicy()) {
851 <      std::string myPolicy = simParams_->getCutoffPolicy();
962 >    if (!myPolicy.empty()){
963 >      toUpper(myPolicy);
964        if (myPolicy == "MIX") {
965          cp = MIX_CUTOFF_POLICY;
966        } else {
# Line 866 | Line 978 | namespace oopse {
978            }    
979          }          
980        }
981 +    }          
982 +    notifyFortranCutoffPolicy(&cp);
983 +
984 +    // Check the Skin Thickness for neighborlists
985 +    RealType skin;
986 +    if (simParams_->haveSkinThickness()) {
987 +      skin = simParams_->getSkinThickness();
988 +      notifyFortranSkinThickness(&skin);
989 +    }            
990 +        
991 +    // Check if the cutoff was set explicitly:
992 +    if (simParams_->haveCutoffRadius()) {
993 +      rcut_ = simParams_->getCutoffRadius();
994 +      if (simParams_->haveSwitchingRadius()) {
995 +        rsw_  = simParams_->getSwitchingRadius();
996 +      } else {
997 +        if (fInfo_.SIM_uses_Charges |
998 +            fInfo_.SIM_uses_Dipoles |
999 +            fInfo_.SIM_uses_RF) {
1000 +          
1001 +          rsw_ = 0.85 * rcut_;
1002 +          sprintf(painCave.errMsg,
1003 +                  "SimCreator Warning: No value was set for the switchingRadius.\n"
1004 +                  "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n"
1005 +                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1006 +        painCave.isFatal = 0;
1007 +        simError();
1008 +        } else {
1009 +          rsw_ = rcut_;
1010 +          sprintf(painCave.errMsg,
1011 +                  "SimCreator Warning: No value was set for the switchingRadius.\n"
1012 +                  "\tOOPSE will use the same value as the cutoffRadius.\n"
1013 +                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1014 +          painCave.isFatal = 0;
1015 +          simError();
1016 +        }
1017 +      }
1018 +      
1019 +      notifyFortranCutoffs(&rcut_, &rsw_);
1020 +      
1021 +    } else {
1022 +      
1023 +      // For electrostatic atoms, we'll assume a large safe value:
1024 +      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
1025 +        sprintf(painCave.errMsg,
1026 +                "SimCreator Warning: No value was set for the cutoffRadius.\n"
1027 +                "\tOOPSE will use a default value of 15.0 angstroms"
1028 +                "\tfor the cutoffRadius.\n");
1029 +        painCave.isFatal = 0;
1030 +        simError();
1031 +        rcut_ = 15.0;
1032 +      
1033 +        if (simParams_->haveElectrostaticSummationMethod()) {
1034 +          std::string myMethod = simParams_->getElectrostaticSummationMethod();
1035 +          toUpper(myMethod);
1036 +          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
1037 +            if (simParams_->haveSwitchingRadius()){
1038 +              sprintf(painCave.errMsg,
1039 +                      "SimInfo Warning: A value was set for the switchingRadius\n"
1040 +                      "\teven though the electrostaticSummationMethod was\n"
1041 +                      "\tset to %s\n", myMethod.c_str());
1042 +              painCave.isFatal = 1;
1043 +              simError();            
1044 +            }
1045 +          }
1046 +        }
1047 +      
1048 +        if (simParams_->haveSwitchingRadius()){
1049 +          rsw_ = simParams_->getSwitchingRadius();
1050 +        } else {        
1051 +          sprintf(painCave.errMsg,
1052 +                  "SimCreator Warning: No value was set for switchingRadius.\n"
1053 +                  "\tOOPSE will use a default value of\n"
1054 +                  "\t0.85 * cutoffRadius for the switchingRadius\n");
1055 +          painCave.isFatal = 0;
1056 +          simError();
1057 +          rsw_ = 0.85 * rcut_;
1058 +        }
1059 +        notifyFortranCutoffs(&rcut_, &rsw_);
1060 +      } else {
1061 +        // We didn't set rcut explicitly, and we don't have electrostatic atoms, so
1062 +        // We'll punt and let fortran figure out the cutoffs later.
1063 +        
1064 +        notifyFortranYouAreOnYourOwn();
1065 +
1066 +      }
1067      }
870    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
1068    }
1069  
1070 <  void SimInfo::setupCoulombicCorrection( int isError ) {    
1070 >  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
1071      
1072      int errorOut;
1073 <    int cc =  NONE;
1074 <    double alphaVal;
1073 >    int esm =  NONE;
1074 >    int sm = UNDAMPED;
1075 >    RealType alphaVal;
1076 >    RealType dielectric;
1077  
1078      errorOut = isError;
1079 +    alphaVal = simParams_->getDampingAlpha();
1080 +    dielectric = simParams_->getDielectric();
1081  
1082 <    if (simParams_->haveCoulombicCorrection()) {
1083 <      std::string myCorrection = simParams_->getCoulombicCorrection();
1084 <      if (myCorrection == "NONE") {
1085 <        cc = NONE;
1082 >    if (simParams_->haveElectrostaticSummationMethod()) {
1083 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
1084 >      toUpper(myMethod);
1085 >      if (myMethod == "NONE") {
1086 >        esm = NONE;
1087        } else {
1088 <        if (myCorrection == "UNDAMPED_WOLF") {
1089 <          cc = UNDAMPED_WOLF;
1088 >        if (myMethod == "SWITCHING_FUNCTION") {
1089 >          esm = SWITCHING_FUNCTION;
1090          } else {
1091 <          if (myCorrection == "WOLF") {            
1092 <            cc = WOLF;
1093 <            if (!simParams_->haveDampingAlpha()) {
1094 <              //throw error
1095 <              sprintf( painCave.errMsg,
894 <                       "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Wolf Coulombic Correction.", simParams_->getDampingAlpha());
895 <              painCave.isFatal = 0;
896 <              simError();
897 <            }
898 <            alphaVal = simParams_->getDampingAlpha();
899 <          } else {
900 <            if (myCorrection == "REACTION_FIELD") {
901 <              cc = REACTION_FIELD;
1091 >          if (myMethod == "SHIFTED_POTENTIAL") {
1092 >            esm = SHIFTED_POTENTIAL;
1093 >          } else {
1094 >            if (myMethod == "SHIFTED_FORCE") {            
1095 >              esm = SHIFTED_FORCE;
1096              } else {
1097 <              // throw error        
1098 <              sprintf( painCave.errMsg,
1099 <                       "SimInfo error: Unknown coulombicCorrection. (Input file specified %s .)\n\tcoulombicCorrection must be one of: \"none\", \"undamped_wolf\", \"wolf\", or \"reaction_field\".", myCorrection.c_str() );
1100 <              painCave.isFatal = 1;
1101 <              simError();
1102 <            }    
1103 <          }          
1097 >              if (myMethod == "REACTION_FIELD") {            
1098 >                esm = REACTION_FIELD;
1099 >              } else {
1100 >                // throw error        
1101 >                sprintf( painCave.errMsg,
1102 >                         "SimInfo error: Unknown electrostaticSummationMethod.\n"
1103 >                         "\t(Input file specified %s .)\n"
1104 >                         "\telectrostaticSummationMethod must be one of: \"none\",\n"
1105 >                         "\t\"shifted_potential\", \"shifted_force\", or \n"
1106 >                         "\t\"reaction_field\".\n", myMethod.c_str() );
1107 >                painCave.isFatal = 1;
1108 >                simError();
1109 >              }    
1110 >            }          
1111 >          }
1112          }
1113        }
1114      }
1115 <    initFortranFF( &fInfo_.SIM_uses_RF, &cc, &alphaVal, &errorOut );
1115 >    
1116 >    if (simParams_->haveElectrostaticScreeningMethod()) {
1117 >      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
1118 >      toUpper(myScreen);
1119 >      if (myScreen == "UNDAMPED") {
1120 >        sm = UNDAMPED;
1121 >      } else {
1122 >        if (myScreen == "DAMPED") {
1123 >          sm = DAMPED;
1124 >          if (!simParams_->haveDampingAlpha()) {
1125 >            //throw error
1126 >            sprintf( painCave.errMsg,
1127 >                     "SimInfo warning: dampingAlpha was not specified in the input file.\n"
1128 >                     "\tA default value of %f (1/ang) will be used.\n", alphaVal);
1129 >            painCave.isFatal = 0;
1130 >            simError();
1131 >          }
1132 >        } else {
1133 >          // throw error        
1134 >          sprintf( painCave.errMsg,
1135 >                   "SimInfo error: Unknown electrostaticScreeningMethod.\n"
1136 >                   "\t(Input file specified %s .)\n"
1137 >                   "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
1138 >                   "or \"damped\".\n", myScreen.c_str() );
1139 >          painCave.isFatal = 1;
1140 >          simError();
1141 >        }
1142 >      }
1143 >    }
1144 >    
1145 >    // let's pass some summation method variables to fortran
1146 >    setElectrostaticSummationMethod( &esm );
1147 >    setFortranElectrostaticMethod( &esm );
1148 >    setScreeningMethod( &sm );
1149 >    setDampingAlpha( &alphaVal );
1150 >    setReactionFieldDielectric( &dielectric );
1151 >    initFortranFF( &errorOut );
1152 >  }
1153 >
1154 >  void SimInfo::setupSwitchingFunction() {    
1155 >    int ft = CUBIC;
1156 >
1157 >    if (simParams_->haveSwitchingFunctionType()) {
1158 >      std::string funcType = simParams_->getSwitchingFunctionType();
1159 >      toUpper(funcType);
1160 >      if (funcType == "CUBIC") {
1161 >        ft = CUBIC;
1162 >      } else {
1163 >        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1164 >          ft = FIFTH_ORDER_POLY;
1165 >        } else {
1166 >          // throw error        
1167 >          sprintf( painCave.errMsg,
1168 >                   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1169 >          painCave.isFatal = 1;
1170 >          simError();
1171 >        }          
1172 >      }
1173 >    }
1174 >
1175 >    // send switching function notification to switcheroo
1176 >    setFunctionType(&ft);
1177 >
1178 >  }
1179 >
1180 >  void SimInfo::setupAccumulateBoxDipole() {    
1181 >
1182 >    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
1183 >    if ( simParams_->haveAccumulateBoxDipole() )
1184 >      if ( simParams_->getAccumulateBoxDipole() ) {
1185 >        setAccumulateBoxDipole();
1186 >        calcBoxDipole_ = true;
1187 >      }
1188 >
1189    }
1190  
1191    void SimInfo::addProperty(GenericData* genData) {
# Line 969 | Line 1244 | namespace oopse {
1244      Molecule* mol;
1245  
1246      Vector3d comVel(0.0);
1247 <    double totalMass = 0.0;
1247 >    RealType totalMass = 0.0;
1248      
1249  
1250      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1251 <      double mass = mol->getMass();
1251 >      RealType mass = mol->getMass();
1252        totalMass += mass;
1253        comVel += mass * mol->getComVel();
1254      }  
1255  
1256   #ifdef IS_MPI
1257 <    double tmpMass = totalMass;
1257 >    RealType tmpMass = totalMass;
1258      Vector3d tmpComVel(comVel);    
1259 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1260 <    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1259 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1260 >    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1261   #endif
1262  
1263      comVel /= totalMass;
# Line 995 | Line 1270 | namespace oopse {
1270      Molecule* mol;
1271  
1272      Vector3d com(0.0);
1273 <    double totalMass = 0.0;
1273 >    RealType totalMass = 0.0;
1274      
1275      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1276 <      double mass = mol->getMass();
1276 >      RealType mass = mol->getMass();
1277        totalMass += mass;
1278        com += mass * mol->getCom();
1279      }  
1280  
1281   #ifdef IS_MPI
1282 <    double tmpMass = totalMass;
1282 >    RealType tmpMass = totalMass;
1283      Vector3d tmpCom(com);    
1284 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1285 <    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1284 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1285 >    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1286   #endif
1287  
1288      com /= totalMass;
# Line 1031 | Line 1306 | namespace oopse {
1306        Molecule* mol;
1307        
1308      
1309 <      double totalMass = 0.0;
1309 >      RealType totalMass = 0.0;
1310      
1311  
1312        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1313 <         double mass = mol->getMass();
1313 >         RealType mass = mol->getMass();
1314           totalMass += mass;
1315           com += mass * mol->getCom();
1316           comVel += mass * mol->getComVel();          
1317        }  
1318        
1319   #ifdef IS_MPI
1320 <      double tmpMass = totalMass;
1320 >      RealType tmpMass = totalMass;
1321        Vector3d tmpCom(com);  
1322        Vector3d tmpComVel(comVel);
1323 <      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1324 <      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1325 <      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1323 >      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1324 >      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1325 >      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1326   #endif
1327        
1328        com /= totalMass;
# Line 1066 | Line 1341 | namespace oopse {
1341     void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1342        
1343  
1344 <      double xx = 0.0;
1345 <      double yy = 0.0;
1346 <      double zz = 0.0;
1347 <      double xy = 0.0;
1348 <      double xz = 0.0;
1349 <      double yz = 0.0;
1344 >      RealType xx = 0.0;
1345 >      RealType yy = 0.0;
1346 >      RealType zz = 0.0;
1347 >      RealType xy = 0.0;
1348 >      RealType xz = 0.0;
1349 >      RealType yz = 0.0;
1350        Vector3d com(0.0);
1351        Vector3d comVel(0.0);
1352        
# Line 1083 | Line 1358 | namespace oopse {
1358        Vector3d thisq(0.0);
1359        Vector3d thisv(0.0);
1360  
1361 <      double thisMass = 0.0;
1361 >      RealType thisMass = 0.0;
1362      
1363        
1364        
# Line 1121 | Line 1396 | namespace oopse {
1396   #ifdef IS_MPI
1397        Mat3x3d tmpI(inertiaTensor);
1398        Vector3d tmpAngMom;
1399 <      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1400 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1399 >      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1400 >      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1401   #endif
1402                
1403        return;
# Line 1143 | Line 1418 | namespace oopse {
1418        Vector3d thisr(0.0);
1419        Vector3d thisp(0.0);
1420        
1421 <      double thisMass;
1421 >      RealType thisMass;
1422        
1423        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1424          thisMass = mol->getMass();
# Line 1156 | Line 1431 | namespace oopse {
1431        
1432   #ifdef IS_MPI
1433        Vector3d tmpAngMom;
1434 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1434 >      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1435   #endif
1436        
1437        return angularMomentum;
1438     }
1439    
1440 <  
1440 >  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1441 >    return IOIndexToIntegrableObject.at(index);
1442 >  }
1443 >  
1444 >  void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
1445 >    IOIndexToIntegrableObject= v;
1446 >  }
1447 >
1448 > /*
1449 >   void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
1450 >      assert( v.size() == nAtoms_ + nRigidBodies_);
1451 >      sdByGlobalIndex_ = v;
1452 >    }
1453 >
1454 >    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
1455 >      //assert(index < nAtoms_ + nRigidBodies_);
1456 >      return sdByGlobalIndex_.at(index);
1457 >    }  
1458 > */  
1459   }//end namespace oopse
1460  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines