[730d7a] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
[0aa122] | 4 | * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
|
---|
[730d7a] | 5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
| 6 | */
|
---|
| 7 |
|
---|
| 8 | /*
|
---|
| 9 | * AdaptivityMap.cpp
|
---|
| 10 | *
|
---|
| 11 | * Created on: Oct 20, 2011
|
---|
| 12 | * Author: heber
|
---|
| 13 | */
|
---|
| 14 |
|
---|
| 15 | #ifdef HAVE_CONFIG_H
|
---|
| 16 | #include <config.h>
|
---|
| 17 | #endif
|
---|
| 18 |
|
---|
| 19 | #include "CodePatterns/MemDebug.hpp"
|
---|
| 20 |
|
---|
| 21 | #include "AdaptivityMap.hpp"
|
---|
| 22 |
|
---|
| 23 | #include <fstream>
|
---|
| 24 |
|
---|
[851be8] | 25 | #include "CodePatterns/Assert.hpp"
|
---|
[730d7a] | 26 | #include "CodePatterns/Log.hpp"
|
---|
| 27 |
|
---|
[6f0841] | 28 | #include "Atom/atom.hpp"
|
---|
[730d7a] | 29 | #include "Helpers/defs.hpp"
|
---|
| 30 | #include "Helpers/helpers.hpp"
|
---|
| 31 | #include "molecule.hpp"
|
---|
| 32 |
|
---|
| 33 | /** Constructor of class AdaptivityMap.
|
---|
| 34 | *
|
---|
| 35 | */
|
---|
| 36 | AdaptivityMap::AdaptivityMap()
|
---|
| 37 | {}
|
---|
| 38 |
|
---|
| 39 | /** Destructor of class AdaptivityMap.
|
---|
| 40 | *
|
---|
| 41 | */
|
---|
| 42 | AdaptivityMap::~AdaptivityMap()
|
---|
| 43 | {}
|
---|
| 44 |
|
---|
| 45 | /** Inserts a (\a No, \a value) pair into the list, overwriting present one.
|
---|
| 46 | * Note if values are equal, No will decided on which is first
|
---|
| 47 | * \param *out output stream for debugging
|
---|
| 48 | * \param &AdaptiveCriteriaList list to insert into
|
---|
| 49 | * \param &IndexedKeySetList list to find key set for a given index \a No
|
---|
| 50 | * \param FragOrder current bond order of fragment
|
---|
| 51 | * \param No index of keyset
|
---|
| 52 | * \param value energy value
|
---|
| 53 | */
|
---|
[851be8] | 54 | void AdaptivityMap::InsertIntoAdaptiveCriteriaList(int FragOrder, int No, double Value)
|
---|
[730d7a] | 55 | {
|
---|
[851be8] | 56 | ASSERT( AdaptiveCriteriaList != NULL,
|
---|
| 57 | "AdaptivityMap::InsertIntoAdaptiveCriteriaList() - AdaptiveCriteriaList is not allocated yet.");
|
---|
[730d7a] | 58 | const_iterator marker = find(No); // find keyset to Frag No.
|
---|
| 59 | if (marker != end()) { // if found
|
---|
| 60 | Value *= 1 + MYEPSILON*(*((*marker).second.begin())); // in case of equal energies this makes them not equal without changing anything actually
|
---|
| 61 | // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
|
---|
| 62 | std::pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList->insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
|
---|
| 63 | std::map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
|
---|
| 64 | if (!InsertedElement.second) { // this root is already present
|
---|
| 65 | if ((*PresentItem).second.second < FragOrder) // if order there is lower, update entry with higher-order term
|
---|
| 66 | //if ((*PresentItem).second.first < (*runner).first) // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
|
---|
| 67 | { // if value is smaller, update value and order
|
---|
| 68 | (*PresentItem).second.first = fabs(Value);
|
---|
| 69 | (*PresentItem).second.second = FragOrder;
|
---|
[47d041] | 70 | LOG(2, "Updated element (" << (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "]).");
|
---|
[730d7a] | 71 | } else {
|
---|
[47d041] | 72 | LOG(2, "Did not update element " << (*PresentItem).first << " as " << FragOrder << " is less than or equal to " << (*PresentItem).second.second << ".");
|
---|
[730d7a] | 73 | }
|
---|
| 74 | } else {
|
---|
[47d041] | 75 | LOG(2, "Inserted element (" << (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "]).");
|
---|
[730d7a] | 76 | }
|
---|
| 77 | } else {
|
---|
[47d041] | 78 | LOG(1, "No Fragment under No. " << No << "found.");
|
---|
[730d7a] | 79 | }
|
---|
| 80 | };
|
---|
| 81 |
|
---|
| 82 |
|
---|
| 83 | /** Scans the adaptive order file and insert (index, value) into map.
|
---|
| 84 | * \param &path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
|
---|
| 85 | */
|
---|
[851be8] | 86 | void AdaptivityMap::ScanAdaptiveFileIntoMap(std::string &path)
|
---|
[730d7a] | 87 | {
|
---|
| 88 | int No = 0, FragOrder = 0;
|
---|
| 89 | double Value = 0.;
|
---|
| 90 | char buffer[MAXSTRINGSIZE];
|
---|
| 91 | std::string filename = path + ENERGYPERFRAGMENT;
|
---|
| 92 | std::ifstream InputFile(filename.c_str());
|
---|
| 93 |
|
---|
| 94 | if (InputFile.fail()) {
|
---|
[47d041] | 95 | ELOG(1, "Cannot find file " << filename << ".");
|
---|
[851be8] | 96 | return;
|
---|
[730d7a] | 97 | }
|
---|
| 98 |
|
---|
| 99 | if (CountLinesinFile(InputFile) > 0) {
|
---|
| 100 | // each line represents a fragment root (Atom::Nr) id and its energy contribution
|
---|
| 101 | InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
|
---|
| 102 | InputFile.getline(buffer, MAXSTRINGSIZE);
|
---|
| 103 | while(!InputFile.eof()) {
|
---|
| 104 | InputFile.getline(buffer, MAXSTRINGSIZE);
|
---|
| 105 | if (strlen(buffer) > 2) {
|
---|
[47d041] | 106 | //LOG(2, "Scanning: " << buffer);
|
---|
[730d7a] | 107 | stringstream line(buffer);
|
---|
| 108 | line >> FragOrder;
|
---|
| 109 | line >> ws >> No;
|
---|
| 110 | line >> ws >> Value; // skip time entry
|
---|
| 111 | line >> ws >> Value;
|
---|
| 112 | No -= 1; // indices start at 1 in file, not 0
|
---|
[47d041] | 113 | //LOG(2, " - yields (" << No << "," << Value << ", " << FragOrder << ")");
|
---|
[730d7a] | 114 |
|
---|
| 115 | // clean the list of those entries that have been superceded by higher order terms already
|
---|
[851be8] | 116 | InsertIntoAdaptiveCriteriaList(FragOrder, No, Value);
|
---|
[730d7a] | 117 | }
|
---|
| 118 | }
|
---|
| 119 | // close and done
|
---|
| 120 | InputFile.close();
|
---|
| 121 | InputFile.clear();
|
---|
| 122 | }
|
---|
| 123 | };
|
---|
| 124 |
|
---|
| 125 | /** Maps adaptive criteria list back onto (Value, (Root Nr., Order))
|
---|
| 126 | * (i.e. sorted by value to pick the highest ones)
|
---|
| 127 | * \param *mol molecule with atoms
|
---|
| 128 | * \return remapped list
|
---|
| 129 | */
|
---|
[851be8] | 130 | void AdaptivityMap::ReMapAdaptiveCriteriaListToValue(molecule *mol)
|
---|
[730d7a] | 131 | {
|
---|
| 132 | atom *Walker = NULL;
|
---|
[851be8] | 133 | ASSERT( AdaptiveCriteriaList != NULL,
|
---|
| 134 | "AdaptivityMap::ReMapAdaptiveCriteriaListToValue() - AdaptiveCriteriaList is not allocated yet.");
|
---|
| 135 | FinalRootCandidates = new AdaptiveCriteriaValueMap;
|
---|
[47d041] | 136 | LOG(1, "Root candidate list is: ");
|
---|
[851be8] | 137 | for(AdaptiveCriteriaIndexMap::const_iterator runner = AdaptiveCriteriaList->begin(); runner != AdaptiveCriteriaList->end(); runner++) {
|
---|
[730d7a] | 138 | Walker = mol->FindAtom((*runner).first);
|
---|
| 139 | if (Walker != NULL) {
|
---|
| 140 | //if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
|
---|
| 141 | if (!Walker->MaxOrder) {
|
---|
[47d041] | 142 | LOG(2, "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])");
|
---|
[730d7a] | 143 | FinalRootCandidates->insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
|
---|
| 144 | } else {
|
---|
[47d041] | 145 | LOG(2, "Excluding (" << *Walker << ", " << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "]), as it has reached its maximum order.");
|
---|
[730d7a] | 146 | }
|
---|
| 147 | } else {
|
---|
[47d041] | 148 | ELOG(0, "Atom No. " << (*runner).second.first << " was not found in this molecule.");
|
---|
[730d7a] | 149 | performCriticalExit();
|
---|
| 150 | }
|
---|
| 151 | }
|
---|
| 152 | };
|
---|
| 153 |
|
---|
| 154 | /** Counts lines in file.
|
---|
| 155 | * Note we are scanning lines from current position, not from beginning.
|
---|
| 156 | * \param InputFile file to be scanned.
|
---|
| 157 | */
|
---|
| 158 | int AdaptivityMap::CountLinesinFile(std::ifstream &InputFile) const
|
---|
| 159 | {
|
---|
| 160 | char *buffer = new char[MAXSTRINGSIZE];
|
---|
| 161 | int lines=0;
|
---|
| 162 |
|
---|
| 163 | int PositionMarker = InputFile.tellg(); // not needed as Inputfile is copied, given by value, not by ref
|
---|
| 164 | // count the number of lines, i.e. the number of fragments
|
---|
| 165 | InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
|
---|
| 166 | InputFile.getline(buffer, MAXSTRINGSIZE);
|
---|
| 167 | while(!InputFile.eof()) {
|
---|
| 168 | InputFile.getline(buffer, MAXSTRINGSIZE);
|
---|
| 169 | lines++;
|
---|
| 170 | }
|
---|
| 171 | InputFile.seekg(PositionMarker, ios::beg);
|
---|
| 172 | delete[](buffer);
|
---|
| 173 | return lines;
|
---|
| 174 | };
|
---|
[851be8] | 175 |
|
---|
| 176 | /** Marks all candidate sites for update if below adaptive threshold.
|
---|
| 177 | * Picks a given number of highest values and set *AtomMask to true.
|
---|
| 178 | * \param *AtomMask defines true/false per global Atom::Nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
|
---|
| 179 | * \param Order desired order
|
---|
| 180 | * \param *mol molecule with atoms
|
---|
| 181 | * \return true - if update is necessary, false - not
|
---|
| 182 | */
|
---|
| 183 | bool AdaptivityMap::MarkUpdateCandidates(bool *AtomMask, int Order, molecule *mol) const
|
---|
| 184 | {
|
---|
| 185 | atom *Walker = NULL;
|
---|
| 186 | int No = -1;
|
---|
| 187 | bool status = false;
|
---|
| 188 | ASSERT( FinalRootCandidates != NULL,
|
---|
| 189 | "AdaptivityMap::MarkUpdateCandidates() - FinalRootCandidates is not allocated yet.");
|
---|
| 190 | for(AdaptiveCriteriaValueMap::const_iterator runner = FinalRootCandidates->upper_bound(pow(10.,Order)); runner != FinalRootCandidates->end(); runner++) {
|
---|
| 191 | No = (*runner).second.first;
|
---|
| 192 | Walker = mol->FindAtom(No);
|
---|
| 193 | //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->getNr()]) {
|
---|
[47d041] | 194 | LOG(2, "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true.");
|
---|
[851be8] | 195 | AtomMask[No] = true;
|
---|
| 196 | status = true;
|
---|
| 197 | //} else
|
---|
[47d041] | 198 | //LOG(2, "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->getNr()] << " does not allow further adaptive increase.");
|
---|
[851be8] | 199 | }
|
---|
| 200 | return status;
|
---|
| 201 | };
|
---|
| 202 |
|
---|
| 203 | /** Checks whether there are any adaptive items currently.
|
---|
| 204 | *
|
---|
| 205 | * @return true - there are items for adaptive refinement, false - there are none
|
---|
| 206 | */
|
---|
| 207 | bool AdaptivityMap::IsAdaptiveCriteriaListEmpty() const
|
---|
| 208 | {
|
---|
| 209 | ASSERT( AdaptiveCriteriaList != NULL,
|
---|
| 210 | "AdaptivityMap::ReMapAdaptiveCriteriaListToValue() - AdaptiveCriteriaList is not allocated yet.");
|
---|
| 211 | return AdaptiveCriteriaList->empty();
|
---|
| 212 | }
|
---|