[0331ee] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2013 Frederik Heber. All rights reserved.
|
---|
| 5 | *
|
---|
| 6 | *
|
---|
| 7 | * This file is part of MoleCuilder.
|
---|
| 8 | *
|
---|
| 9 | * MoleCuilder is free software: you can redistribute it and/or modify
|
---|
| 10 | * it under the terms of the GNU General Public License as published by
|
---|
| 11 | * the Free Software Foundation, either version 2 of the License, or
|
---|
| 12 | * (at your option) any later version.
|
---|
| 13 | *
|
---|
| 14 | * MoleCuilder is distributed in the hope that it will be useful,
|
---|
| 15 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 16 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 17 | * GNU General Public License for more details.
|
---|
| 18 | *
|
---|
| 19 | * You should have received a copy of the GNU General Public License
|
---|
| 20 | * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
|
---|
| 21 | */
|
---|
| 22 |
|
---|
| 23 | /*
|
---|
| 24 | * Interfragmenter.cpp
|
---|
| 25 | *
|
---|
| 26 | * Created on: Jul 5, 2013
|
---|
| 27 | * Author: heber
|
---|
| 28 | */
|
---|
| 29 |
|
---|
| 30 | // include config.h
|
---|
| 31 | #ifdef HAVE_CONFIG_H
|
---|
| 32 | #include <config.h>
|
---|
| 33 | #endif
|
---|
| 34 |
|
---|
| 35 | #include "CodePatterns/MemDebug.hpp"
|
---|
| 36 |
|
---|
| 37 | #include "Interfragmenter.hpp"
|
---|
| 38 |
|
---|
[d410e25] | 39 | #include <algorithm>
|
---|
| 40 |
|
---|
[0331ee] | 41 | #include "CodePatterns/Assert.hpp"
|
---|
| 42 | #include "CodePatterns/Log.hpp"
|
---|
| 43 |
|
---|
| 44 | #include "LinearAlgebra/Vector.hpp"
|
---|
| 45 |
|
---|
[cee9e8] | 46 | #include "Descriptors/AtomDescriptor.hpp"
|
---|
[0331ee] | 47 | #include "Element/element.hpp"
|
---|
| 48 | #include "Fragmentation/Graph.hpp"
|
---|
| 49 | #include "Fragmentation/KeySet.hpp"
|
---|
| 50 | #include "LinkedCell/LinkedCell_View.hpp"
|
---|
| 51 | #include "LinkedCell/types.hpp"
|
---|
| 52 | #include "World.hpp"
|
---|
| 53 |
|
---|
[0d9053] | 54 | static Vector getAtomIdSetCenter(
|
---|
| 55 | const AtomIdSet &_atoms)
|
---|
| 56 | {
|
---|
| 57 | const molecule * const _mol = (*_atoms.begin())->getMolecule();
|
---|
| 58 | const size_t atoms_size = _atoms.getAtomIds().size();
|
---|
| 59 | Vector center;
|
---|
| 60 | for (AtomIdSet::const_iterator iter = _atoms.begin();
|
---|
| 61 | iter != _atoms.end(); ++iter) {
|
---|
| 62 | center += (*iter)->getPosition();
|
---|
| 63 | ASSERT ( _mol == (*iter)->getMolecule(),
|
---|
| 64 | "getAtomIdSetCenter() - ids in same keyset belong to different molecule.");
|
---|
| 65 | }
|
---|
| 66 | center *= 1./(double)atoms_size;
|
---|
| 67 |
|
---|
| 68 | return center;
|
---|
| 69 | }
|
---|
| 70 |
|
---|
| 71 | Interfragmenter::candidates_t Interfragmenter::getNeighborsOutsideMolecule(
|
---|
| 72 | const AtomIdSet &_atoms,
|
---|
| 73 | double _Rcut,
|
---|
| 74 | const enum HydrogenTreatment _treatment) const
|
---|
| 75 | {
|
---|
| 76 | /// go through linked cell and get all neighboring atoms up to Rcut
|
---|
| 77 | const LinkedCell::LinkedCell_View view = World::getInstance().getLinkedCell(_Rcut);
|
---|
| 78 | const Vector center = getAtomIdSetCenter(_atoms);
|
---|
| 79 | const LinkedCell::LinkedList neighbors = view.getAllNeighbors(_Rcut, center);
|
---|
| 80 | LOG(4, "DEBUG: Obtained " << neighbors.size() << " neighbors in distance of "
|
---|
| 81 | << _Rcut << " around " << center << ".");
|
---|
| 82 |
|
---|
| 83 | /// remove all atoms that belong to same molecule as does one of the
|
---|
| 84 | /// fragment's atoms
|
---|
| 85 | const molecule * const _mol = (*_atoms.begin())->getMolecule();
|
---|
| 86 | candidates_t candidates;
|
---|
| 87 | candidates.reserve(neighbors.size());
|
---|
| 88 | for (LinkedCell::LinkedList::const_iterator iter = neighbors.begin();
|
---|
| 89 | iter != neighbors.end(); ++iter) {
|
---|
| 90 | const atom * const _atom = static_cast<const atom * const >(*iter);
|
---|
| 91 | ASSERT( _atom != NULL,
|
---|
| 92 | "Interfragmenter::getNeighborsOutsideMolecule() - a neighbor is not actually an atom?");
|
---|
| 93 | if ((_atom->getMolecule() != _mol)
|
---|
| 94 | && (_atom->getPosition().DistanceSquared(center) < _Rcut*_Rcut)
|
---|
| 95 | && ((_treatment == IncludeHydrogen) || (_atom->getType()->getAtomicNumber() != 1))) {
|
---|
| 96 | candidates.push_back(_atom);
|
---|
| 97 | }
|
---|
| 98 | }
|
---|
| 99 | LOG(3, "DEBUG: There remain " << candidates.size() << " candidates.");
|
---|
| 100 |
|
---|
| 101 | return candidates;
|
---|
| 102 | }
|
---|
| 103 |
|
---|
| 104 | void Interfragmenter::combineFragments(
|
---|
[b1c5f46] | 105 | const size_t _MaxOrder,
|
---|
[0d9053] | 106 | const candidates_t &_candidates,
|
---|
[bd6e5c] | 107 | const atomkeyset_t &_fragmentmap,
|
---|
[0d9053] | 108 | const KeySet &_keyset,
|
---|
| 109 | Graph &_InterFragments,
|
---|
| 110 | int &_counter)
|
---|
| 111 | {
|
---|
| 112 | for (candidates_t::const_iterator candidateiter = _candidates.begin();
|
---|
| 113 | candidateiter != _candidates.end(); ++candidateiter) {
|
---|
| 114 | const atom *_atom = *candidateiter;
|
---|
| 115 | LOG(3, "DEBUG: Current candidate is " << *_atom << ".");
|
---|
[bd6e5c] | 116 | atomkeyset_t::const_iterator finditer = _fragmentmap.find(_atom->getId());
|
---|
[0d9053] | 117 | ASSERT( finditer != _fragmentmap.end(),
|
---|
| 118 | "Interfragmenter::combineFragments() - could not find atom "
|
---|
| 119 | +toString(_atom->getId())+" in fragment specific lookup.");
|
---|
[bd6e5c] | 120 | // copy set to allow erase
|
---|
| 121 | keysets_t othersets(finditer->second);
|
---|
[b1c5f46] | 122 | ASSERT( !othersets.empty(),
|
---|
| 123 | "Interfragmenter::combineFragments() - keysets to "+toString(_atom->getId())+
|
---|
| 124 | "is empty.");
|
---|
[0d9053] | 125 | keysets_t::iterator otheriter = othersets.begin();
|
---|
| 126 | while (otheriter != othersets.end()) {
|
---|
[bd6e5c] | 127 | const KeySet &otherset = *otheriter;
|
---|
[0d9053] | 128 | LOG(3, "DEBUG: Current keyset is " << otherset << ".");
|
---|
| 129 | // only add them one way round and not the other
|
---|
| 130 | if (otherset < _keyset) {
|
---|
| 131 | ++otheriter;
|
---|
| 132 | continue;
|
---|
| 133 | }
|
---|
[b1c5f46] | 134 | // only add if combined they don't exceed the desired maxorder
|
---|
| 135 | if (otherset.size() + _keyset.size() > _MaxOrder) {
|
---|
| 136 | LOG(3, "INFO: Rejecting " << otherset << " as in sum their orders exceed " << _MaxOrder);
|
---|
| 137 | ++otheriter;
|
---|
| 138 | continue;
|
---|
| 139 | }
|
---|
[0d9053] | 140 | KeySet newset(otherset);
|
---|
| 141 | newset.insert(_keyset.begin(), _keyset.end());
|
---|
| 142 | LOG(3, "DEBUG: Inserting new combined set " << newset << ".");
|
---|
| 143 | _InterFragments.insert( std::make_pair(newset, std::make_pair(++_counter, 1.)));
|
---|
| 144 | // finally, remove the set such that no other combination exists
|
---|
| 145 | otheriter = othersets.erase(otheriter);
|
---|
| 146 | }
|
---|
| 147 | }
|
---|
| 148 | }
|
---|
| 149 |
|
---|
| 150 | void Interfragmenter::operator()(
|
---|
[d410e25] | 151 | Graph &TotalGraph,
|
---|
[0d9053] | 152 | const size_t MaxOrder,
|
---|
| 153 | const double Rcut,
|
---|
| 154 | const enum HydrogenTreatment treatment)
|
---|
| 155 | {
|
---|
[d713ce] | 156 | AtomFragmentsMap& atomfragments_container = AtomFragmentsMap::getInstance();
|
---|
| 157 | const atomkeyset_t &atomkeyset = atomfragments_container.getMap();
|
---|
[0d9053] | 158 |
|
---|
[0331ee] | 159 | Graph InterFragments;
|
---|
[d410e25] | 160 | int counter = atomkeyset.size();
|
---|
[0331ee] | 161 |
|
---|
| 162 | /// go through all fragments up to MaxOrder
|
---|
| 163 | LOG(1, "INFO: Creating inter-fragments.");
|
---|
[d410e25] | 164 | for (atomkeyset_t::const_iterator atomiter = atomkeyset.begin();
|
---|
| 165 | atomiter != atomkeyset.end(); ++atomiter) {
|
---|
| 166 | const atomId_t &atomid = atomiter->first;
|
---|
| 167 | LOG(2, "DEBUG: Current atomid is " << atomid);
|
---|
| 168 |
|
---|
| 169 | const AtomFragmentsMap::keysets_t &keysets = atomiter->second;
|
---|
| 170 | for (AtomFragmentsMap::keysets_t::const_iterator keyiter = keysets.begin();
|
---|
| 171 | keyiter != keysets.end(); ++keyiter) {
|
---|
| 172 | const KeySet &keyset = *keyiter;
|
---|
| 173 | const AtomIdSet atoms(keyset);
|
---|
| 174 | const size_t atoms_size = atoms.getAtomIds().size();
|
---|
| 175 | if ((atoms_size > MaxOrder) || (atoms_size == 0))
|
---|
| 176 | continue;
|
---|
| 177 |
|
---|
| 178 | // get neighboring atoms outside the current molecule
|
---|
| 179 | candidates_t candidates = getNeighborsOutsideMolecule(atoms, Rcut, treatment);
|
---|
| 180 |
|
---|
| 181 | // create a lookup specific to this fragment
|
---|
| 182 | std::vector<atomId_t> atomids(candidates.size());
|
---|
| 183 | std::transform(
|
---|
| 184 | candidates.begin(), candidates.end(),
|
---|
| 185 | atomids.begin(),
|
---|
| 186 | boost::bind(&atom::getId, _1));
|
---|
| 187 | atomkeyset_t fragmentmap = atomfragments_container.getMap(atomids, MaxOrder);
|
---|
| 188 |
|
---|
| 189 | /// combine each remaining fragment with current fragment to a new fragment
|
---|
| 190 | /// if keyset is less (to prevent addding same inter-fragment twice)
|
---|
| 191 | combineFragments(MaxOrder, candidates, fragmentmap, keyset, InterFragments, counter);
|
---|
| 192 | }
|
---|
[0331ee] | 193 | }
|
---|
| 194 |
|
---|
| 195 | /// eventually, add all new fragments to the Graph
|
---|
[d410e25] | 196 | counter = atomkeyset.size();
|
---|
[0331ee] | 197 | TotalGraph.InsertGraph(InterFragments, counter);
|
---|
| 198 | }
|
---|
[cee9e8] | 199 |
|
---|
| 200 | double Interfragmenter::findLargestCutoff(
|
---|
| 201 | const size_t _MaxOrder,
|
---|
| 202 | const double _upperbound,
|
---|
| 203 | const enum HydrogenTreatment _treatment) const
|
---|
| 204 | {
|
---|
| 205 | double Rcut = _upperbound*_upperbound;
|
---|
| 206 | std::pair<atomId_t, atomId_t> ClosestPair;
|
---|
| 207 | // place all atoms into LC grid with some upper bound
|
---|
| 208 | const LinkedCell::LinkedCell_View view = World::getInstance().getLinkedCell(_upperbound);
|
---|
| 209 |
|
---|
| 210 | // go through each atom and find closest atom not in the same keyset
|
---|
[d410e25] | 211 | AtomFragmentsMap& atomfragments_container = AtomFragmentsMap::getInstance();
|
---|
| 212 | const atomkeyset_t &atomkeyset = atomfragments_container.getMap();
|
---|
| 213 | for (atomkeyset_t::const_iterator atomiter = atomkeyset.begin();
|
---|
| 214 | atomiter != atomkeyset.end(); ++atomiter) {
|
---|
| 215 | const atomId_t &atomid = atomiter->first;
|
---|
| 216 | LOG(2, "DEBUG: Current atomid is " << atomid);
|
---|
| 217 |
|
---|
| 218 | const AtomFragmentsMap::keysets_t &keysets = atomiter->second;
|
---|
| 219 | for (AtomFragmentsMap::keysets_t::const_iterator keyiter = keysets.begin();
|
---|
| 220 | keyiter != keysets.end(); ++keyiter) {
|
---|
| 221 | const KeySet &keyset = *keyiter;
|
---|
| 222 | const AtomIdSet atoms(keyset);
|
---|
| 223 |
|
---|
| 224 | // get neighboring atoms outside the current molecule
|
---|
| 225 | const candidates_t candidates = getNeighborsOutsideMolecule(atoms, _upperbound, _treatment);
|
---|
| 226 | const Vector center = getAtomIdSetCenter(atoms);
|
---|
| 227 |
|
---|
| 228 | for (candidates_t::const_iterator candidateiter = candidates.begin();
|
---|
| 229 | candidateiter != candidates.end(); ++candidateiter) {
|
---|
| 230 | atom const * const Walker = *candidateiter;
|
---|
| 231 | // go through each atom in set and pick minimal distance
|
---|
| 232 | for (AtomIdSet::const_iterator setiter = atoms.begin(); setiter != atoms.end(); ++setiter) {
|
---|
| 233 | const double distanceSquared = Walker->getPosition().DistanceSquared(center);
|
---|
| 234 | // pick the smallest compared to current Rcut if smaller
|
---|
| 235 | if (distanceSquared < Rcut) {
|
---|
| 236 | Rcut = distanceSquared;
|
---|
| 237 | ClosestPair.first = (*setiter)->getId();
|
---|
| 238 | ClosestPair.second = Walker->getId();
|
---|
| 239 | LOG(2, "DEBUG: Found new pair " << ClosestPair << " with distance " << sqrt(Rcut));
|
---|
| 240 | }
|
---|
[cee9e8] | 241 | }
|
---|
| 242 | }
|
---|
| 243 | }
|
---|
| 244 | }
|
---|
[17e4fd] | 245 | const double largest_distance = sqrt(Rcut);
|
---|
| 246 | LOG(1, "INFO: Largest inter-fragment distance to cause no additional fragments: "
|
---|
| 247 | << largest_distance);
|
---|
[cee9e8] | 248 |
|
---|
[17e4fd] | 249 | return largest_distance;
|
---|
[cee9e8] | 250 | }
|
---|