Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F91686486
MoleculeSet.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Nov 13, 11:26
Size
9 KB
Mime Type
text/x-c
Expires
Fri, Nov 15, 11:26 (2 d)
Engine
blob
Format
Raw Data
Handle
22305775
Attached To
rLAMMPS lammps
MoleculeSet.cpp
View Options
// ATC_Method headers
#include "MoleculeSet.h"
#include "ATC_Method.h"
#include "LammpsInterface.h"
#include "ATC_Error.h"
#include <queue>
#include <utility>
#include <sstream>
using std::multimap;
using std::map;
using std::pair;
using std::set;
using std::stringstream;
namespace ATC {
//--------------------------------------------------------
//--------------------------------------------------------
// Class MoleculeSet
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
MoleculeSet::MoleculeSet(ATC_Method * atc, int groupBit) :
atc_(atc),
groupBit_(groupBit),
lammps_(atc->lammps_interface())
{
// do nothing
}
//--------------------------------------------------------
// Destructor
//--------------------------------------------------------
MoleculeSet::~MoleculeSet()
{
clear();
}
//--------------------------------------------------------
// clear
//--------------------------------------------------------
void MoleculeSet::clear()
{
moleculeToAtoms_.clear();
localMoleculeToAtoms_.clear();
}
//--------------------------------------------------------
// initialize
//--------------------------------------------------------
void MoleculeSet::initialize(map<int,double> * globalAtomsPerMolecule)
{
// determine the total number of molecules in this group
// essentially ripped Compute::molecules_in_group from lammps
// find lo/hi molecule ID for any atom in group
int i;
int *molecule = lammps_->atom_to_molecule();
const int *mask = lammps_->atom_mask();
int nlocal = lammps_->nlocal();
int lo = lammps_->natoms();
int hi = -1;
int flag = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupBit_) {
if (molecule[i] == 0) flag = 1;
lo = MIN(lo,molecule[i]);
hi = MAX(hi,molecule[i]);
}
}
int flagall;
lammps_->int_allsum(&flag,&flagall);
if (flagall) throw ATC_Error("Atom with molecule ID = 0 included in atc molecule group");
int globalLo, globalHi;
lammps_->int_allmin(&lo,&globalLo);
lammps_->int_allmax(&hi,&globalHi);
if (globalLo == lammps_->natoms()) throw ATC_Error("MoleculeSet:initialize - no molecules correspond to the group");
// molmap = vector of length nlen
// set to 1 for IDs that appear in group across all procs, else 0
int nlen = globalHi-globalLo+1;
int * localCount = new int[nlen];
for (i = 0; i < nlen; i++) localCount[i] = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupBit_)
localCount[molecule[i]-globalLo]++;
int * globalCount = new int[nlen];
lammps_->int_allsum(localCount,globalCount,nlen);
// nmolecules = # of non-zero IDs in molmap
// molmap[i] = index of molecule, skipping molecules not in group with -1
nMoleculesTotal_ = 0;
for (i = 0; i < nlen; i++)
if (globalCount[i]) nMoleculesTotal_++;
if (globalAtomsPerMolecule) {
for (i = 0; i < nlen; i++)
if (globalCount[i]) globalAtomsPerMolecule->insert(pair<int,double>(i+globalLo,double(globalCount[i])));
}
// deallocate storage
delete [] localCount;
delete [] globalCount;
}
//--------------------------------------------------------
// atoms_by_global_molecule
//--------------------------------------------------------
set<int> MoleculeSet::atoms_by_global_molecule(int id) const
{
if (need_reset()) reset();
typedef multimap<int, set<int> >::const_iterator MMIT;
typedef set<int>::const_iterator SIT;
pair<MMIT,MMIT> mol = moleculeToAtoms_.equal_range(id);
set<int> realAtoms;
const set<int> * myAtoms;
for (MMIT molIt = mol.first; molIt != mol.second; molIt++) {
myAtoms = &(molIt->second);
for (SIT atomIt = myAtoms->begin(); atomIt != myAtoms->end(); atomIt++)
realAtoms.insert(*atomIt);
}
return realAtoms;
}
//--------------------------------------------------------
// atoms_by_local_molecule
//--------------------------------------------------------
const set<int> & MoleculeSet::atoms_by_local_molecule(int id) const
{
if (need_reset()) reset();
return localMoleculeToAtoms_[id]->second;
}
//--------------------------------------------------------
// set_local_molecules_to_atoms
//--------------------------------------------------------
void MoleculeSet::set_local_molecules_to_atoms() const
{
localMoleculeToAtoms_.clear();
localMoleculeToAtoms_.reserve(moleculeToAtoms_.size());
multimap<int, set<int> >::const_iterator molecule;
for (molecule = moleculeToAtoms_.begin(); molecule != moleculeToAtoms_.end(); molecule++)
localMoleculeToAtoms_.push_back(molecule);
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class SmallMoleculeSet
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
SmallMoleculeSet::SmallMoleculeSet(ATC_Method * atc, int groupBit,
PerAtomQuantity<int> * bondList, PerAtomQuantity<int> * numBond) :
MoleculeSet(atc,groupBit),
bondList_(bondList)
{
}
//--------------------------------------------------------
// Destructor
//--------------------------------------------------------
SmallMoleculeSet::~SmallMoleculeSet()
{
}
//--------------------------------------------------------
// clear
//--------------------------------------------------------
void SmallMoleculeSet::clear()
{
MoleculeSet::clear();
}
//--------------------------------------------------------
// initialize
//--------------------------------------------------------
void SmallMoleculeSet::initialize(std::map<int, double> * globalAtomsPerMolecule)
{
// make sure newton_bond is off, otherwise use large molecule set
if (lammps_->newton_bond())
throw ATC_Error("Cannot use newton_bond with small molecules");
MoleculeSet::initialize(globalAtomsPerMolecule);
}
//--------------------------------------------------------
// reset
//--------------------------------------------------------
void SmallMoleculeSet::reset() const
{
stringstream message;
using std::queue;
moleculeToAtoms_.clear();
lammps_->forward_comm_fix();
int * numBond = lammps_->num_bond();
int ** bondAtom = lammps_->bond_atom();
// add in real atoms for molecules
int *molecule = lammps_->atom_to_molecule();
const int *mask = lammps_->atom_mask();
int nlocal = lammps_->nlocal();
_atomFound_.resize(atc_->nproc_ghost());
_atomFound_ = false;
int nmol = 0;
for (int i = 0; i < nlocal; i++) {
queue<int> myQueue;
if ((mask[i] & groupBit_) && !_atomFound_(i)) {
set<int> myAtoms;
myAtoms.insert(i);
_atomFound_(i) = true;
for (int j = 0; j < numBond[i]; j++) {
int localIdx = lammps_->local_to_global_map(bondAtom[i][j]);
if (!_atomFound_(localIdx)) {
myQueue.push(localIdx);
_atomFound_(localIdx) = true;
}
}
while (!myQueue.empty()) {
int myIdx = myQueue.front();
myQueue.pop();
myAtoms.insert(myIdx);
for (int j = 0; j < numBond[myIdx]; j++) {
int localIdx = lammps_->local_to_global_map(bondAtom[myIdx][j]);
if (!_atomFound_(localIdx)) {
myQueue.push(localIdx);
_atomFound_(localIdx) = true;
}
}
}
nmol++;
moleculeToAtoms_.insert(pair<int,set<int> >(molecule[i],myAtoms));
}
}
// set local molecule order
MoleculeSet::set_local_molecules_to_atoms();
needReset_ = false;
}
//--------------------------------------------------------
// atoms_by_global_molecule
//--------------------------------------------------------
set<int> SmallMoleculeSet::atoms_by_global_molecule(int id) const
{
// take all atoms and prune out ghosts
set<int> realAtoms = MoleculeSet::atoms_by_global_molecule(id);
remove_proc_ghosts(realAtoms);
return realAtoms;
}
//--------------------------------------------------------
// local_fraction
//--------------------------------------------------------
double SmallMoleculeSet::local_fraction(int id) const
{
if (need_reset()) reset();
set<int> realAtoms = (localMoleculeToAtoms_[id])->second;
int totalAtoms = realAtoms.size();
remove_proc_ghosts(realAtoms);
return double(realAtoms.size())/double(totalAtoms);
}
//--------------------------------------------------------
// remove_proc_ghosts
//--------------------------------------------------------
void SmallMoleculeSet::remove_proc_ghosts(set<int> & atomSet) const
{
int nlocalIdx = lammps_->nlocal() - 1;
set<int>::const_iterator atomIt;
for (atomIt = atomSet.begin(); atomIt != atomSet.end(); atomIt++) {
if (*atomIt > nlocalIdx)
atomSet.erase(atomIt);
}
}
};
Event Timeline
Log In to Comment