| 56 | 
  | 
#include "UseTheForce/notifyCutoffs_interface.h" | 
| 57 | 
  | 
#include "utils/MemoryUtils.hpp" | 
| 58 | 
  | 
#include "utils/simError.h" | 
| 59 | 
+ | 
#include "selection/SelectionManager.hpp" | 
| 60 | 
  | 
 | 
| 61 | 
  | 
#ifdef IS_MPI | 
| 62 | 
  | 
#include "UseTheForce/mpiComponentPlan.h" | 
| 65 | 
  | 
 | 
| 66 | 
  | 
namespace oopse { | 
| 67 | 
  | 
 | 
| 68 | 
< | 
SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,  | 
| 68 | 
> | 
SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,  | 
| 69 | 
  | 
                                ForceField* ff, Globals* simParams) :  | 
| 70 | 
< | 
                                forceField_(ff), simParams_(simParams),  | 
| 70 | 
> | 
                                stamps_(stamps), forceField_(ff), simParams_(simParams),  | 
| 71 | 
  | 
                                ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 72 | 
  | 
                                nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),  | 
| 73 | 
  | 
                                nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 109 | 
  | 
 | 
| 110 | 
  | 
        //calculate atoms in rigid bodies | 
| 111 | 
  | 
        int nAtomsInRigidBodies = 0; | 
| 112 | 
< | 
        int nRigidBodiesInStamp = molStamp->getNCutoffGroups(); | 
| 112 | 
> | 
        int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | 
| 113 | 
  | 
         | 
| 114 | 
  | 
        for (int j=0; j < nRigidBodiesInStamp; j++) { | 
| 115 | 
  | 
            rbStamp = molStamp->getRigidBody(j); | 
| 138 | 
  | 
#ifdef IS_MPI     | 
| 139 | 
  | 
    molToProcMap_.resize(nGlobalMols_); | 
| 140 | 
  | 
#endif | 
| 141 | 
< | 
     | 
| 141 | 
> | 
 | 
| 142 | 
  | 
} | 
| 143 | 
  | 
 | 
| 144 | 
  | 
SimInfo::~SimInfo() { | 
| 145 | 
< | 
    //MemoryUtils::deleteVectorOfPointer(molecules_); | 
| 146 | 
< | 
 | 
| 147 | 
< | 
    MemoryUtils::deleteVectorOfPointer(moleculeStamps_); | 
| 148 | 
< | 
     | 
| 145 | 
> | 
    std::map<int, Molecule*>::iterator i; | 
| 146 | 
> | 
    for (i = molecules_.begin(); i != molecules_.end(); ++i) { | 
| 147 | 
> | 
        delete i->second; | 
| 148 | 
> | 
    } | 
| 149 | 
> | 
    molecules_.clear(); | 
| 150 | 
> | 
        | 
| 151 | 
> | 
    delete stamps_; | 
| 152 | 
  | 
    delete sman_; | 
| 153 | 
  | 
    delete simParams_; | 
| 154 | 
  | 
    delete forceField_; | 
| 151 | 
– | 
 | 
| 155 | 
  | 
} | 
| 156 | 
  | 
 | 
| 157 | 
  | 
int SimInfo::getNGlobalConstraints() { | 
| 366 | 
  | 
        exclude_.addPair(c, d);         | 
| 367 | 
  | 
    } | 
| 368 | 
  | 
 | 
| 369 | 
< | 
     | 
| 369 | 
> | 
    Molecule::RigidBodyIterator rbIter; | 
| 370 | 
> | 
    RigidBody* rb; | 
| 371 | 
> | 
    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 372 | 
> | 
        std::vector<Atom*> atoms = rb->getAtoms(); | 
| 373 | 
> | 
        for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 374 | 
> | 
            for (int j = i + 1; j < atoms.size(); ++j) { | 
| 375 | 
> | 
                a = atoms[i]->getGlobalIndex(); | 
| 376 | 
> | 
                b = atoms[j]->getGlobalIndex(); | 
| 377 | 
> | 
                exclude_.addPair(a, b); | 
| 378 | 
> | 
            } | 
| 379 | 
> | 
        } | 
| 380 | 
> | 
    }         | 
| 381 | 
> | 
 | 
| 382 | 
  | 
} | 
| 383 | 
  | 
 | 
| 384 | 
  | 
void SimInfo::removeExcludePairs(Molecule* mol) { | 
| 422 | 
  | 
        exclude_.removePair(b, d); | 
| 423 | 
  | 
        exclude_.removePair(c, d);         | 
| 424 | 
  | 
    } | 
| 425 | 
+ | 
 | 
| 426 | 
+ | 
    Molecule::RigidBodyIterator rbIter; | 
| 427 | 
+ | 
    RigidBody* rb; | 
| 428 | 
+ | 
    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 429 | 
+ | 
        std::vector<Atom*> atoms = rb->getAtoms(); | 
| 430 | 
+ | 
        for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 431 | 
+ | 
            for (int j = i + 1; j < atoms.size(); ++j) { | 
| 432 | 
+ | 
                a = atoms[i]->getGlobalIndex(); | 
| 433 | 
+ | 
                b = atoms[j]->getGlobalIndex(); | 
| 434 | 
+ | 
                exclude_.removePair(a, b); | 
| 435 | 
+ | 
            } | 
| 436 | 
+ | 
        } | 
| 437 | 
+ | 
    }         | 
| 438 | 
  | 
 | 
| 439 | 
  | 
} | 
| 440 | 
  | 
 | 
| 670 | 
  | 
    } | 
| 671 | 
  | 
     | 
| 672 | 
  | 
    //setup fortran simulation | 
| 645 | 
– | 
    //gloalExcludes and molMembershipArray should go away (They are never used) | 
| 646 | 
– | 
    //why the hell fortran need to know molecule? | 
| 647 | 
– | 
    //OOPSE = Object-Obfuscated Parallel Simulation Engine | 
| 673 | 
  | 
    int nGlobalExcludes = 0; | 
| 674 | 
  | 
    int* globalExcludes = NULL;  | 
| 675 | 
  | 
    int* excludeList = exclude_.getExcludeList(); | 
| 776 | 
  | 
    return maxCutoffRadius; | 
| 777 | 
  | 
} | 
| 778 | 
  | 
 | 
| 779 | 
< | 
void SimInfo::setupCutoff() { | 
| 755 | 
< | 
    double rcut_;  //cutoff radius | 
| 756 | 
< | 
    double rsw_; //switching radius | 
| 779 | 
> | 
void SimInfo::getCutoff(double& rcut, double& rsw) { | 
| 780 | 
  | 
     | 
| 781 | 
  | 
    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 782 | 
  | 
         | 
| 787 | 
  | 
                "\tfor the cutoffRadius.\n"); | 
| 788 | 
  | 
            painCave.isFatal = 0; | 
| 789 | 
  | 
            simError(); | 
| 790 | 
< | 
            rcut_ = 15.0; | 
| 790 | 
> | 
            rcut = 15.0; | 
| 791 | 
  | 
        } else{ | 
| 792 | 
< | 
            rcut_ = simParams_->getRcut(); | 
| 792 | 
> | 
            rcut = simParams_->getRcut(); | 
| 793 | 
  | 
        } | 
| 794 | 
  | 
 | 
| 795 | 
  | 
        if (!simParams_->haveRsw()){ | 
| 799 | 
  | 
                "\t0.95 * cutoffRadius for the switchingRadius\n"); | 
| 800 | 
  | 
            painCave.isFatal = 0; | 
| 801 | 
  | 
            simError(); | 
| 802 | 
< | 
            rsw_ = 0.95 * rcut_; | 
| 802 | 
> | 
            rsw = 0.95 * rcut; | 
| 803 | 
  | 
        } else{ | 
| 804 | 
< | 
            rsw_ = simParams_->getRsw(); | 
| 804 | 
> | 
            rsw = simParams_->getRsw(); | 
| 805 | 
  | 
        } | 
| 806 | 
  | 
 | 
| 807 | 
  | 
    } else { | 
| 809 | 
  | 
        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | 
| 810 | 
  | 
         | 
| 811 | 
  | 
        if (simParams_->haveRcut()) { | 
| 812 | 
< | 
            rcut_ = simParams_->getRcut(); | 
| 812 | 
> | 
            rcut = simParams_->getRcut(); | 
| 813 | 
  | 
        } else { | 
| 814 | 
  | 
            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | 
| 815 | 
< | 
            rcut_ = calcMaxCutoffRadius(); | 
| 815 | 
> | 
            rcut = calcMaxCutoffRadius(); | 
| 816 | 
  | 
        } | 
| 817 | 
  | 
 | 
| 818 | 
  | 
        if (simParams_->haveRsw()) { | 
| 819 | 
< | 
            rsw_  = simParams_->getRsw(); | 
| 819 | 
> | 
            rsw  = simParams_->getRsw(); | 
| 820 | 
  | 
        } else { | 
| 821 | 
< | 
            rsw_ = rcut_; | 
| 821 | 
> | 
            rsw = rcut; | 
| 822 | 
  | 
        } | 
| 823 | 
  | 
     | 
| 824 | 
  | 
    } | 
| 825 | 
< | 
         | 
| 825 | 
> | 
} | 
| 826 | 
> | 
 | 
| 827 | 
> | 
void SimInfo::setupCutoff() { | 
| 828 | 
> | 
    getCutoff(rcut_, rsw_);     | 
| 829 | 
  | 
    double rnblist = rcut_ + 1; // skin of neighbor list | 
| 830 | 
  | 
 | 
| 831 | 
  | 
    //Pass these cutoff radius etc. to fortran. This function should be called once and only once | 
| 857 | 
  | 
} | 
| 858 | 
  | 
 | 
| 859 | 
  | 
void SimInfo::setSnapshotManager(SnapshotManager* sman) { | 
| 860 | 
+ | 
    if (sman_ == sman) { | 
| 861 | 
+ | 
        return; | 
| 862 | 
+ | 
    }     | 
| 863 | 
+ | 
    delete sman_; | 
| 864 | 
  | 
    sman_ = sman; | 
| 865 | 
  | 
 | 
| 866 | 
  | 
    Molecule* mol; |