# | 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" |
62 | > | #include "UseTheForce/DarkSide/switcheroo_interface.h" |
63 | #include "utils/MemoryUtils.hpp" | |
64 | #include "utils/simError.h" | |
65 | #include "selection/SelectionManager.hpp" | |
# | Line 67 | Line 70 | namespace oopse { | |
70 | #endif | |
71 | ||
72 | namespace oopse { | |
73 | + | std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { |
74 | + | std::map<int, std::set<int> >::iterator i = container.find(index); |
75 | + | std::set<int> result; |
76 | + | if (i != container.end()) { |
77 | + | result = i->second; |
78 | + | } |
79 | ||
80 | < | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, |
81 | < | ForceField* ff, Globals* simParams) : |
82 | < | stamps_(stamps), forceField_(ff), simParams_(simParams), |
80 | > | return result; |
81 | > | } |
82 | > | |
83 | > | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : |
84 | > | forceField_(ff), simParams_(simParams), |
85 | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | |
86 | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | |
87 | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | |
# | Line 78 | Line 89 | namespace oopse { | |
89 | nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), | |
90 | sman_(NULL), fortranInitialized_(false) { | |
91 | ||
81 | – | |
82 | – | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; |
92 | MoleculeStamp* molStamp; | |
93 | int nMolWithSameStamp; | |
94 | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | |
# | Line 87 | Line 96 | namespace oopse { | |
96 | CutoffGroupStamp* cgStamp; | |
97 | RigidBodyStamp* rbStamp; | |
98 | int nRigidAtoms = 0; | |
99 | < | |
100 | < | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { |
101 | < | molStamp = i->first; |
102 | < | nMolWithSameStamp = i->second; |
99 | > | std::vector<Component*> components = simParams->getComponents(); |
100 | > | |
101 | > | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
102 | > | molStamp = (*i)->getMoleculeStamp(); |
103 | > | nMolWithSameStamp = (*i)->getNMol(); |
104 | ||
105 | addMoleculeStamp(molStamp, nMolWithSameStamp); | |
106 | ||
107 | //calculate atoms in molecules | |
108 | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | |
109 | ||
100 | – | |
110 | //calculate atoms in cutoff groups | |
111 | int nAtomsInGroups = 0; | |
112 | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | |
113 | ||
114 | for (int j=0; j < nCutoffGroupsInStamp; j++) { | |
115 | < | cgStamp = molStamp->getCutoffGroup(j); |
115 | > | cgStamp = molStamp->getCutoffGroupStamp(j); |
116 | nAtomsInGroups += cgStamp->getNMembers(); | |
117 | } | |
118 | ||
# | Line 116 | Line 125 | namespace oopse { | |
125 | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | |
126 | ||
127 | for (int j=0; j < nRigidBodiesInStamp; j++) { | |
128 | < | rbStamp = molStamp->getRigidBody(j); |
128 | > | rbStamp = molStamp->getRigidBodyStamp(j); |
129 | nAtomsInRigidBodies += rbStamp->getNMembers(); | |
130 | } | |
131 | ||
# | Line 155 | Line 164 | namespace oopse { | |
164 | } | |
165 | molecules_.clear(); | |
166 | ||
158 | – | delete stamps_; |
167 | delete sman_; | |
168 | delete simParams_; | |
169 | delete forceField_; | |
# | Line 262 | Line 270 | namespace oopse { | |
270 | } | |
271 | } | |
272 | ||
273 | < | }//end for (integrableObject) |
274 | < | }// end for (mol) |
273 | > | } |
274 | > | } |
275 | ||
276 | // n_constraints is local, so subtract them on each processor | |
277 | ndf_local -= nConstraints_; | |
# | Line 342 | Line 350 | namespace oopse { | |
350 | int b; | |
351 | int c; | |
352 | int d; | |
353 | + | |
354 | + | std::map<int, std::set<int> > atomGroups; |
355 | + | |
356 | + | Molecule::RigidBodyIterator rbIter; |
357 | + | RigidBody* rb; |
358 | + | Molecule::IntegrableObjectIterator ii; |
359 | + | StuntDouble* integrableObject; |
360 | ||
361 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
362 | + | integrableObject = mol->nextIntegrableObject(ii)) { |
363 | + | |
364 | + | if (integrableObject->isRigidBody()) { |
365 | + | rb = static_cast<RigidBody*>(integrableObject); |
366 | + | std::vector<Atom*> atoms = rb->getAtoms(); |
367 | + | std::set<int> rigidAtoms; |
368 | + | for (int i = 0; i < atoms.size(); ++i) { |
369 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); |
370 | + | } |
371 | + | for (int i = 0; i < atoms.size(); ++i) { |
372 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
373 | + | } |
374 | + | } else { |
375 | + | std::set<int> oneAtomSet; |
376 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
377 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
378 | + | } |
379 | + | } |
380 | + | |
381 | + | |
382 | + | |
383 | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | |
384 | a = bond->getAtomA()->getGlobalIndex(); | |
385 | b = bond->getAtomB()->getGlobalIndex(); | |
# | Line 353 | Line 390 | namespace oopse { | |
390 | a = bend->getAtomA()->getGlobalIndex(); | |
391 | b = bend->getAtomB()->getGlobalIndex(); | |
392 | c = bend->getAtomC()->getGlobalIndex(); | |
393 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
394 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
395 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
396 | ||
397 | < | exclude_.addPair(a, b); |
398 | < | exclude_.addPair(a, c); |
399 | < | exclude_.addPair(b, c); |
397 | > | exclude_.addPairs(rigidSetA, rigidSetB); |
398 | > | exclude_.addPairs(rigidSetA, rigidSetC); |
399 | > | exclude_.addPairs(rigidSetB, rigidSetC); |
400 | > | |
401 | > | //exclude_.addPair(a, b); |
402 | > | //exclude_.addPair(a, c); |
403 | > | //exclude_.addPair(b, c); |
404 | } | |
405 | ||
406 | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | |
# | Line 364 | Line 408 | namespace oopse { | |
408 | b = torsion->getAtomB()->getGlobalIndex(); | |
409 | c = torsion->getAtomC()->getGlobalIndex(); | |
410 | d = torsion->getAtomD()->getGlobalIndex(); | |
411 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
412 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
413 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
414 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); |
415 | ||
416 | + | exclude_.addPairs(rigidSetA, rigidSetB); |
417 | + | exclude_.addPairs(rigidSetA, rigidSetC); |
418 | + | exclude_.addPairs(rigidSetA, rigidSetD); |
419 | + | exclude_.addPairs(rigidSetB, rigidSetC); |
420 | + | exclude_.addPairs(rigidSetB, rigidSetD); |
421 | + | exclude_.addPairs(rigidSetC, rigidSetD); |
422 | + | |
423 | + | /* |
424 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); |
425 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); |
426 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); |
427 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); |
428 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); |
429 | + | exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); |
430 | + | |
431 | + | |
432 | exclude_.addPair(a, b); | |
433 | exclude_.addPair(a, c); | |
434 | exclude_.addPair(a, d); | |
435 | exclude_.addPair(b, c); | |
436 | exclude_.addPair(b, d); | |
437 | exclude_.addPair(c, d); | |
438 | + | */ |
439 | } | |
440 | ||
376 | – | Molecule::RigidBodyIterator rbIter; |
377 | – | RigidBody* rb; |
441 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
442 | std::vector<Atom*> atoms = rb->getAtoms(); | |
443 | for (int i = 0; i < atoms.size() -1 ; ++i) { | |
# | Line 399 | Line 462 | namespace oopse { | |
462 | int b; | |
463 | int c; | |
464 | int d; | |
465 | + | |
466 | + | std::map<int, std::set<int> > atomGroups; |
467 | + | |
468 | + | Molecule::RigidBodyIterator rbIter; |
469 | + | RigidBody* rb; |
470 | + | Molecule::IntegrableObjectIterator ii; |
471 | + | StuntDouble* integrableObject; |
472 | ||
473 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
474 | + | integrableObject = mol->nextIntegrableObject(ii)) { |
475 | + | |
476 | + | if (integrableObject->isRigidBody()) { |
477 | + | rb = static_cast<RigidBody*>(integrableObject); |
478 | + | std::vector<Atom*> atoms = rb->getAtoms(); |
479 | + | std::set<int> rigidAtoms; |
480 | + | for (int i = 0; i < atoms.size(); ++i) { |
481 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); |
482 | + | } |
483 | + | for (int i = 0; i < atoms.size(); ++i) { |
484 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
485 | + | } |
486 | + | } else { |
487 | + | std::set<int> oneAtomSet; |
488 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
489 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
490 | + | } |
491 | + | } |
492 | + | |
493 | + | |
494 | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | |
495 | a = bond->getAtomA()->getGlobalIndex(); | |
496 | b = bond->getAtomB()->getGlobalIndex(); | |
# | Line 411 | Line 502 | namespace oopse { | |
502 | b = bend->getAtomB()->getGlobalIndex(); | |
503 | c = bend->getAtomC()->getGlobalIndex(); | |
504 | ||
505 | < | exclude_.removePair(a, b); |
506 | < | exclude_.removePair(a, c); |
507 | < | exclude_.removePair(b, c); |
505 | > | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
506 | > | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
507 | > | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
508 | > | |
509 | > | exclude_.removePairs(rigidSetA, rigidSetB); |
510 | > | exclude_.removePairs(rigidSetA, rigidSetC); |
511 | > | exclude_.removePairs(rigidSetB, rigidSetC); |
512 | > | |
513 | > | //exclude_.removePair(a, b); |
514 | > | //exclude_.removePair(a, c); |
515 | > | //exclude_.removePair(b, c); |
516 | } | |
517 | ||
518 | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | |
# | Line 422 | Line 521 | namespace oopse { | |
521 | c = torsion->getAtomC()->getGlobalIndex(); | |
522 | d = torsion->getAtomD()->getGlobalIndex(); | |
523 | ||
524 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
525 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
526 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
527 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); |
528 | + | |
529 | + | exclude_.removePairs(rigidSetA, rigidSetB); |
530 | + | exclude_.removePairs(rigidSetA, rigidSetC); |
531 | + | exclude_.removePairs(rigidSetA, rigidSetD); |
532 | + | exclude_.removePairs(rigidSetB, rigidSetC); |
533 | + | exclude_.removePairs(rigidSetB, rigidSetD); |
534 | + | exclude_.removePairs(rigidSetC, rigidSetD); |
535 | + | |
536 | + | /* |
537 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); |
538 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); |
539 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); |
540 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); |
541 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); |
542 | + | exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); |
543 | + | |
544 | + | |
545 | exclude_.removePair(a, b); | |
546 | exclude_.removePair(a, c); | |
547 | exclude_.removePair(a, d); | |
548 | exclude_.removePair(b, c); | |
549 | exclude_.removePair(b, d); | |
550 | exclude_.removePair(c, d); | |
551 | + | */ |
552 | } | |
553 | ||
433 | – | Molecule::RigidBodyIterator rbIter; |
434 | – | RigidBody* rb; |
554 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
555 | std::vector<Atom*> atoms = rb->getAtoms(); | |
556 | for (int i = 0; i < atoms.size() -1 ; ++i) { | |
# | Line 471 | Line 590 | namespace oopse { | |
590 | int isError = 0; | |
591 | ||
592 | setupElectrostaticSummationMethod( isError ); | |
593 | + | setupSwitchingFunction(); |
594 | ||
595 | if(isError){ | |
596 | sprintf( painCave.errMsg, | |
# | Line 515 | Line 635 | namespace oopse { | |
635 | int useLennardJones = 0; | |
636 | int useElectrostatic = 0; | |
637 | int useEAM = 0; | |
638 | + | int useSC = 0; |
639 | int useCharge = 0; | |
640 | int useDirectional = 0; | |
641 | int useDipole = 0; | |
# | Line 528 | Line 649 | namespace oopse { | |
649 | //usePBC and useRF are from simParams | |
650 | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | |
651 | int useRF; | |
652 | + | int useSF; |
653 | std::string myMethod; | |
654 | ||
655 | // set the useRF logical | |
656 | useRF = 0; | |
657 | + | useSF = 0; |
658 | ||
659 | ||
660 | if (simParams_->haveElectrostaticSummationMethod()) { | |
# | Line 539 | Line 662 | namespace oopse { | |
662 | toUpper(myMethod); | |
663 | if (myMethod == "REACTION_FIELD") { | |
664 | useRF=1; | |
665 | + | } else { |
666 | + | if (myMethod == "SHIFTED_FORCE") { |
667 | + | useSF = 1; |
668 | + | } |
669 | } | |
670 | } | |
671 | ||
# | Line 547 | Line 674 | namespace oopse { | |
674 | useLennardJones |= (*i)->isLennardJones(); | |
675 | useElectrostatic |= (*i)->isElectrostatic(); | |
676 | useEAM |= (*i)->isEAM(); | |
677 | + | useSC |= (*i)->isSC(); |
678 | useCharge |= (*i)->isCharge(); | |
679 | useDirectional |= (*i)->isDirectional(); | |
680 | useDipole |= (*i)->isDipole(); | |
# | Line 597 | Line 725 | namespace oopse { | |
725 | temp = useEAM; | |
726 | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
727 | ||
728 | + | temp = useSC; |
729 | + | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
730 | + | |
731 | temp = useShape; | |
732 | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
733 | ||
# | Line 605 | Line 736 | namespace oopse { | |
736 | ||
737 | temp = useRF; | |
738 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
739 | + | |
740 | + | temp = useSF; |
741 | + | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
742 | ||
743 | #endif | |
744 | ||
# | Line 618 | Line 752 | namespace oopse { | |
752 | fInfo_.SIM_uses_StickyPower = useStickyPower; | |
753 | fInfo_.SIM_uses_GayBerne = useGayBerne; | |
754 | fInfo_.SIM_uses_EAM = useEAM; | |
755 | + | fInfo_.SIM_uses_SC = useSC; |
756 | fInfo_.SIM_uses_Shapes = useShape; | |
757 | fInfo_.SIM_uses_FLARB = useFLARB; | |
758 | fInfo_.SIM_uses_RF = useRF; | |
759 | + | fInfo_.SIM_uses_SF = useSF; |
760 | ||
761 | if( myMethod == "REACTION_FIELD") { | |
762 | ||
# | Line 635 | Line 771 | namespace oopse { | |
771 | simError(); | |
772 | } | |
773 | } | |
774 | + | |
775 | } | |
776 | ||
777 | void SimInfo::setupFortranSim() { | |
# | Line 778 | Line 915 | namespace oopse { | |
915 | ||
916 | ||
917 | } | |
781 | – | |
782 | – | #endif |
783 | – | |
784 | – | double SimInfo::calcMaxCutoffRadius() { |
785 | – | |
786 | – | |
787 | – | std::set<AtomType*> atomTypes; |
788 | – | std::set<AtomType*>::iterator i; |
789 | – | std::vector<double> cutoffRadius; |
790 | – | |
791 | – | //get the unique atom types |
792 | – | atomTypes = getUniqueAtomTypes(); |
918 | ||
794 | – | //query the max cutoff radius among these atom types |
795 | – | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
796 | – | cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); |
797 | – | } |
798 | – | |
799 | – | double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); |
800 | – | #ifdef IS_MPI |
801 | – | //pick the max cutoff radius among the processors |
919 | #endif | |
920 | ||
921 | < | return maxCutoffRadius; |
805 | < | } |
806 | < | |
807 | < | void SimInfo::getCutoff(double& rcut, double& rsw) { |
921 | > | void SimInfo::setupCutoff() { |
922 | ||
923 | < | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
810 | < | |
811 | < | if (!simParams_->haveCutoffRadius()){ |
812 | < | sprintf(painCave.errMsg, |
813 | < | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
814 | < | "\tOOPSE will use a default value of 15.0 angstroms" |
815 | < | "\tfor the cutoffRadius.\n"); |
816 | < | painCave.isFatal = 0; |
817 | < | simError(); |
818 | < | rcut = 15.0; |
819 | < | } else{ |
820 | < | rcut = simParams_->getCutoffRadius(); |
821 | < | } |
822 | < | |
823 | < | if (!simParams_->haveSwitchingRadius()){ |
824 | < | sprintf(painCave.errMsg, |
825 | < | "SimCreator Warning: No value was set for switchingRadius.\n" |
826 | < | "\tOOPSE will use a default value of\n" |
827 | < | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
828 | < | painCave.isFatal = 0; |
829 | < | simError(); |
830 | < | rsw = 0.95 * rcut; |
831 | < | } else{ |
832 | < | rsw = simParams_->getSwitchingRadius(); |
833 | < | } |
834 | < | |
835 | < | } else { |
836 | < | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in |
837 | < | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used |
838 | < | |
839 | < | if (simParams_->haveCutoffRadius()) { |
840 | < | rcut = simParams_->getCutoffRadius(); |
841 | < | } else { |
842 | < | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system |
843 | < | rcut = calcMaxCutoffRadius(); |
844 | < | } |
845 | < | |
846 | < | if (simParams_->haveSwitchingRadius()) { |
847 | < | rsw = simParams_->getSwitchingRadius(); |
848 | < | } else { |
849 | < | rsw = rcut; |
850 | < | } |
851 | < | |
852 | < | } |
853 | < | } |
854 | < | |
855 | < | void SimInfo::setupCutoff() { |
856 | < | getCutoff(rcut_, rsw_); |
857 | < | double rnblist = rcut_ + 1; // skin of neighbor list |
858 | < | |
859 | < | //Pass these cutoff radius etc. to fortran. This function should be called once and only once |
860 | < | |
923 | > | // Check the cutoff policy |
924 | int cp = TRADITIONAL_CUTOFF_POLICY; | |
925 | if (simParams_->haveCutoffPolicy()) { | |
926 | std::string myPolicy = simParams_->getCutoffPolicy(); | |
# | Line 879 | Line 942 | namespace oopse { | |
942 | } | |
943 | } | |
944 | } | |
945 | < | } |
945 | > | } |
946 | > | notifyFortranCutoffPolicy(&cp); |
947 | ||
948 | < | |
948 | > | // Check the Skin Thickness for neighborlists |
949 | > | double skin; |
950 | if (simParams_->haveSkinThickness()) { | |
951 | < | double skinThickness = simParams_->getSkinThickness(); |
952 | < | } |
951 | > | skin = simParams_->getSkinThickness(); |
952 | > | notifyFortranSkinThickness(&skin); |
953 | > | } |
954 | > | |
955 | > | // Check if the cutoff was set explicitly: |
956 | > | if (simParams_->haveCutoffRadius()) { |
957 | > | rcut_ = simParams_->getCutoffRadius(); |
958 | > | if (simParams_->haveSwitchingRadius()) { |
959 | > | rsw_ = simParams_->getSwitchingRadius(); |
960 | > | } else { |
961 | > | rsw_ = rcut_; |
962 | > | } |
963 | > | notifyFortranCutoffs(&rcut_, &rsw_); |
964 | > | |
965 | > | } else { |
966 | > | |
967 | > | // For electrostatic atoms, we'll assume a large safe value: |
968 | > | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
969 | > | sprintf(painCave.errMsg, |
970 | > | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
971 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
972 | > | "\tfor the cutoffRadius.\n"); |
973 | > | painCave.isFatal = 0; |
974 | > | simError(); |
975 | > | rcut_ = 15.0; |
976 | > | |
977 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
978 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
979 | > | toUpper(myMethod); |
980 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
981 | > | if (simParams_->haveSwitchingRadius()){ |
982 | > | sprintf(painCave.errMsg, |
983 | > | "SimInfo Warning: A value was set for the switchingRadius\n" |
984 | > | "\teven though the electrostaticSummationMethod was\n" |
985 | > | "\tset to %s\n", myMethod.c_str()); |
986 | > | painCave.isFatal = 1; |
987 | > | simError(); |
988 | > | } |
989 | > | } |
990 | > | } |
991 | > | |
992 | > | if (simParams_->haveSwitchingRadius()){ |
993 | > | rsw_ = simParams_->getSwitchingRadius(); |
994 | > | } else { |
995 | > | sprintf(painCave.errMsg, |
996 | > | "SimCreator Warning: No value was set for switchingRadius.\n" |
997 | > | "\tOOPSE will use a default value of\n" |
998 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
999 | > | painCave.isFatal = 0; |
1000 | > | simError(); |
1001 | > | rsw_ = 0.85 * rcut_; |
1002 | > | } |
1003 | > | notifyFortranCutoffs(&rcut_, &rsw_); |
1004 | > | } else { |
1005 | > | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so |
1006 | > | // We'll punt and let fortran figure out the cutoffs later. |
1007 | > | |
1008 | > | notifyFortranYouAreOnYourOwn(); |
1009 | ||
1010 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); |
1011 | < | // also send cutoff notification to electrostatics |
891 | < | setElectrostaticCutoffRadius(&rcut_, &rsw_); |
1010 | > | } |
1011 | > | } |
1012 | } | |
1013 | ||
1014 | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | |
1015 | ||
1016 | int errorOut; | |
1017 | int esm = NONE; | |
1018 | + | int sm = UNDAMPED; |
1019 | double alphaVal; | |
1020 | double dielectric; | |
1021 | ||
# | Line 908 | Line 1029 | namespace oopse { | |
1029 | if (myMethod == "NONE") { | |
1030 | esm = NONE; | |
1031 | } else { | |
1032 | < | if (myMethod == "UNDAMPED_WOLF") { |
1033 | < | esm = UNDAMPED_WOLF; |
1032 | > | if (myMethod == "SWITCHING_FUNCTION") { |
1033 | > | esm = SWITCHING_FUNCTION; |
1034 | } else { | |
1035 | < | if (myMethod == "DAMPED_WOLF") { |
1036 | < | esm = DAMPED_WOLF; |
1037 | < | if (!simParams_->haveDampingAlpha()) { |
1038 | < | //throw error |
1039 | < | sprintf( painCave.errMsg, |
919 | < | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", alphaVal); |
920 | < | painCave.isFatal = 0; |
921 | < | simError(); |
922 | < | } |
923 | < | } else { |
924 | < | if (myMethod == "REACTION_FIELD") { |
925 | < | esm = REACTION_FIELD; |
1035 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
1036 | > | esm = SHIFTED_POTENTIAL; |
1037 | > | } else { |
1038 | > | if (myMethod == "SHIFTED_FORCE") { |
1039 | > | esm = SHIFTED_FORCE; |
1040 | } else { | |
1041 | < | // throw error |
1042 | < | sprintf( painCave.errMsg, |
1043 | < | "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() ); |
1044 | < | painCave.isFatal = 1; |
1045 | < | simError(); |
1046 | < | } |
1047 | < | } |
1041 | > | if (myMethod == "REACTION_FIELD") { |
1042 | > | esm = REACTION_FIELD; |
1043 | > | } else { |
1044 | > | // throw error |
1045 | > | sprintf( painCave.errMsg, |
1046 | > | "SimInfo error: Unknown electrostaticSummationMethod.\n" |
1047 | > | "\t(Input file specified %s .)\n" |
1048 | > | "\telectrostaticSummationMethod must be one of: \"none\",\n" |
1049 | > | "\t\"shifted_potential\", \"shifted_force\", or \n" |
1050 | > | "\t\"reaction_field\".\n", myMethod.c_str() ); |
1051 | > | painCave.isFatal = 1; |
1052 | > | simError(); |
1053 | > | } |
1054 | > | } |
1055 | > | } |
1056 | } | |
1057 | } | |
1058 | } | |
1059 | + | |
1060 | + | if (simParams_->haveElectrostaticScreeningMethod()) { |
1061 | + | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
1062 | + | toUpper(myScreen); |
1063 | + | if (myScreen == "UNDAMPED") { |
1064 | + | sm = UNDAMPED; |
1065 | + | } else { |
1066 | + | if (myScreen == "DAMPED") { |
1067 | + | sm = DAMPED; |
1068 | + | if (!simParams_->haveDampingAlpha()) { |
1069 | + | //throw error |
1070 | + | sprintf( painCave.errMsg, |
1071 | + | "SimInfo warning: dampingAlpha was not specified in the input file.\n" |
1072 | + | "\tA default value of %f (1/ang) will be used.\n", alphaVal); |
1073 | + | painCave.isFatal = 0; |
1074 | + | simError(); |
1075 | + | } |
1076 | + | } else { |
1077 | + | // throw error |
1078 | + | sprintf( painCave.errMsg, |
1079 | + | "SimInfo error: Unknown electrostaticScreeningMethod.\n" |
1080 | + | "\t(Input file specified %s .)\n" |
1081 | + | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" |
1082 | + | "or \"damped\".\n", myScreen.c_str() ); |
1083 | + | painCave.isFatal = 1; |
1084 | + | simError(); |
1085 | + | } |
1086 | + | } |
1087 | + | } |
1088 | + | |
1089 | // let's pass some summation method variables to fortran | |
1090 | setElectrostaticSummationMethod( &esm ); | |
1091 | < | setDampedWolfAlpha( &alphaVal ); |
1091 | > | notifyFortranElectrostaticMethod( &esm ); |
1092 | > | setScreeningMethod( &sm ); |
1093 | > | setDampingAlpha( &alphaVal ); |
1094 | setReactionFieldDielectric( &dielectric ); | |
1095 | < | initFortranFF( &esm, &errorOut ); |
1095 | > | initFortranFF( &errorOut ); |
1096 | } | |
1097 | ||
1098 | + | void SimInfo::setupSwitchingFunction() { |
1099 | + | int ft = CUBIC; |
1100 | + | |
1101 | + | if (simParams_->haveSwitchingFunctionType()) { |
1102 | + | std::string funcType = simParams_->getSwitchingFunctionType(); |
1103 | + | toUpper(funcType); |
1104 | + | if (funcType == "CUBIC") { |
1105 | + | ft = CUBIC; |
1106 | + | } else { |
1107 | + | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
1108 | + | ft = FIFTH_ORDER_POLY; |
1109 | + | } else { |
1110 | + | // throw error |
1111 | + | sprintf( painCave.errMsg, |
1112 | + | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
1113 | + | painCave.isFatal = 1; |
1114 | + | simError(); |
1115 | + | } |
1116 | + | } |
1117 | + | } |
1118 | + | |
1119 | + | // send switching function notification to switcheroo |
1120 | + | setFunctionType(&ft); |
1121 | + | |
1122 | + | } |
1123 | + | |
1124 | void SimInfo::addProperty(GenericData* genData) { | |
1125 | properties_.addProperty(genData); | |
1126 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |