[bcf653] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2010 University of Bonn. All rights reserved.
|
---|
| 5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
| 6 | */
|
---|
| 7 |
|
---|
[c4d4df] | 8 | /*
|
---|
| 9 | * analysis.cpp
|
---|
| 10 | *
|
---|
| 11 | * Created on: Oct 13, 2009
|
---|
| 12 | * Author: heber
|
---|
| 13 | */
|
---|
| 14 |
|
---|
[bf3817] | 15 | // include config.h
|
---|
| 16 | #ifdef HAVE_CONFIG_H
|
---|
| 17 | #include <config.h>
|
---|
| 18 | #endif
|
---|
| 19 |
|
---|
[ad011c] | 20 | #include "CodePatterns/MemDebug.hpp"
|
---|
[112b09] | 21 |
|
---|
[c4d4df] | 22 | #include <iostream>
|
---|
[36166d] | 23 | #include <iomanip>
|
---|
[c4d4df] | 24 |
|
---|
[be945c] | 25 | #include "atom.hpp"
|
---|
[129204] | 26 | #include "Bond/bond.hpp"
|
---|
[d74077] | 27 | #include "BoundaryTriangleSet.hpp"
|
---|
[be945c] | 28 | #include "Box.hpp"
|
---|
[c4d4df] | 29 | #include "element.hpp"
|
---|
[ad011c] | 30 | #include "CodePatterns/Info.hpp"
|
---|
| 31 | #include "CodePatterns/Log.hpp"
|
---|
[ea430a] | 32 | #include "Formula.hpp"
|
---|
[c4d4df] | 33 | #include "molecule.hpp"
|
---|
| 34 | #include "tesselation.hpp"
|
---|
| 35 | #include "tesselationhelpers.hpp"
|
---|
[8db598] | 36 | #include "triangleintersectionlist.hpp"
|
---|
[be945c] | 37 | #include "World.hpp"
|
---|
[57f243] | 38 | #include "LinearAlgebra/Vector.hpp"
|
---|
[cca9ef] | 39 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
[ad011c] | 40 | #include "CodePatterns/Verbose.hpp"
|
---|
[b34306] | 41 | #include "World.hpp"
|
---|
[84c494] | 42 | #include "Box.hpp"
|
---|
[c4d4df] | 43 |
|
---|
[be945c] | 44 | #include "analysis_correlation.hpp"
|
---|
| 45 |
|
---|
| 46 | /** Calculates the dipole vector of a given atomSet.
|
---|
| 47 | *
|
---|
| 48 | * Note that we use the following procedure as rule of thumb:
|
---|
| 49 | * -# go through every bond of the atom
|
---|
[d1912f] | 50 | * -# calculate the difference of electronegativities \f$\Delta\mathrm{EN}\f$
|
---|
| 51 | * -# if \f$\Delta\mathrm{EN} > 0.5\f$, we align the bond vector in direction of the more negative element
|
---|
[be945c] | 52 | * -# sum up all vectors
|
---|
| 53 | * -# finally, divide by the number of summed vectors
|
---|
| 54 | *
|
---|
| 55 | * @param atomsbegin begin iterator of atomSet
|
---|
| 56 | * @param atomsend end iterator of atomset
|
---|
| 57 | * @return dipole vector
|
---|
| 58 | */
|
---|
| 59 | Vector getDipole(molecule::const_iterator atomsbegin, molecule::const_iterator atomsend)
|
---|
| 60 | {
|
---|
| 61 | Vector DipoleVector;
|
---|
| 62 | size_t SumOfVectors = 0;
|
---|
| 63 | // go through all atoms
|
---|
| 64 | for (molecule::const_iterator atomiter = atomsbegin;
|
---|
| 65 | atomiter != atomsend;
|
---|
| 66 | ++atomiter) {
|
---|
| 67 | // go through all bonds
|
---|
[9d83b6] | 68 | const BondList& ListOfBonds = (*atomiter)->getListOfBonds();
|
---|
| 69 | for (BondList::const_iterator bonditer = ListOfBonds.begin();
|
---|
| 70 | bonditer != ListOfBonds.end();
|
---|
[be945c] | 71 | ++bonditer) {
|
---|
| 72 | const atom * Otheratom = (*bonditer)->GetOtherAtom(*atomiter);
|
---|
| 73 | if (Otheratom->getId() > (*atomiter)->getId()) {
|
---|
| 74 | const double DeltaEN = (*atomiter)->getType()->getElectronegativity()
|
---|
| 75 | -Otheratom->getType()->getElectronegativity();
|
---|
| 76 | Vector BondDipoleVector = (*atomiter)->getPosition() - Otheratom->getPosition();
|
---|
| 77 | // DeltaEN is always positive, gives correct orientation of vector
|
---|
| 78 | BondDipoleVector.Normalize();
|
---|
| 79 | BondDipoleVector *= DeltaEN;
|
---|
| 80 | DipoleVector += BondDipoleVector;
|
---|
| 81 | SumOfVectors++;
|
---|
| 82 | }
|
---|
| 83 | }
|
---|
| 84 | }
|
---|
| 85 | DipoleVector *= 1./(double)SumOfVectors;
|
---|
| 86 | DoLog(1) && (Log() << Verbose(1) << "Resulting dipole vector is " << DipoleVector << std::endl);
|
---|
| 87 |
|
---|
| 88 | return DipoleVector;
|
---|
| 89 | };
|
---|
| 90 |
|
---|
[ea430a] | 91 | /** Calculates the dipole angular correlation for given molecule type.
|
---|
| 92 | * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
|
---|
[be945c] | 93 | * Angles are given in degrees.
|
---|
[ea430a] | 94 | * \param *molecules vector of molecules
|
---|
| 95 | * \return Map of doubles with values the pair of the two atoms.
|
---|
| 96 | */
|
---|
[be945c] | 97 | DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<molecule *> &molecules)
|
---|
[ea430a] | 98 | {
|
---|
| 99 | Info FunctionInfo(__func__);
|
---|
[caa30b] | 100 | DipoleAngularCorrelationMap *outmap = new DipoleAngularCorrelationMap;
|
---|
[ea430a] | 101 | // double distance = 0.;
|
---|
| 102 | // Box &domain = World::getInstance().getDomain();
|
---|
| 103 | //
|
---|
[be945c] | 104 | if (molecules.empty()) {
|
---|
| 105 | DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
|
---|
| 106 | return outmap;
|
---|
| 107 | }
|
---|
| 108 |
|
---|
| 109 | outmap = new DipoleAngularCorrelationMap;
|
---|
| 110 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin();
|
---|
[92e5cb] | 111 | MolWalker != molecules.end(); ++MolWalker) {
|
---|
[be945c] | 112 | DoLog(2) && (Log()<< Verbose(2) << "Current molecule is "
|
---|
| 113 | << (*MolWalker)->getId() << "." << endl);
|
---|
| 114 | const Vector Dipole = getDipole((*MolWalker)->begin(), (*MolWalker)->end());
|
---|
[92e5cb] | 115 | std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker;
|
---|
| 116 | for (++MolOtherWalker;
|
---|
[be945c] | 117 | MolOtherWalker != molecules.end();
|
---|
[92e5cb] | 118 | ++MolOtherWalker) {
|
---|
[be945c] | 119 | DoLog(2) && (Log() << Verbose(2) << "Current other molecule is "
|
---|
| 120 | << (*MolOtherWalker)->getId() << "." << endl);
|
---|
| 121 | const Vector OtherDipole = getDipole((*MolOtherWalker)->begin(), (*MolOtherWalker)->end());
|
---|
| 122 | const double angle = Dipole.Angle(OtherDipole) * (180./M_PI);
|
---|
| 123 | DoLog(1) && (Log() << Verbose(1) << "Angle is " << angle << "." << endl);
|
---|
| 124 | outmap->insert ( make_pair (angle, make_pair ((*MolWalker), (*MolOtherWalker)) ) );
|
---|
| 125 | }
|
---|
| 126 | }
|
---|
[ea430a] | 127 | return outmap;
|
---|
| 128 | };
|
---|
| 129 |
|
---|
[c4d4df] | 130 |
|
---|
| 131 | /** Calculates the pair correlation between given elements.
|
---|
| 132 | * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
|
---|
[e65de8] | 133 | * \param *molecules vector of molecules
|
---|
[c78d44] | 134 | * \param &elements vector of elements to correlate
|
---|
[c4d4df] | 135 | * \return Map of doubles with values the pair of the two atoms.
|
---|
| 136 | */
|
---|
[e5c0a1] | 137 | PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements)
|
---|
[c4d4df] | 138 | {
|
---|
[3930eb] | 139 | Info FunctionInfo(__func__);
|
---|
[caa30b] | 140 | PairCorrelationMap *outmap = new PairCorrelationMap;
|
---|
[c4d4df] | 141 | double distance = 0.;
|
---|
[014475] | 142 | Box &domain = World::getInstance().getDomain();
|
---|
[c4d4df] | 143 |
|
---|
[e65de8] | 144 | if (molecules.empty()) {
|
---|
[58ed4a] | 145 | DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
|
---|
[c4d4df] | 146 | return outmap;
|
---|
| 147 | }
|
---|
[e65de8] | 148 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
|
---|
[009607e] | 149 | (*MolWalker)->doCountAtoms();
|
---|
[c78d44] | 150 |
|
---|
| 151 | // create all possible pairs of elements
|
---|
[e5c0a1] | 152 | set <pair<const element *,const element *> > PairsOfElements;
|
---|
[c78d44] | 153 | if (elements.size() >= 2) {
|
---|
[e5c0a1] | 154 | for (vector<const element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
|
---|
| 155 | for (vector<const element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
|
---|
[c78d44] | 156 | if (type1 != type2) {
|
---|
[e5c0a1] | 157 | PairsOfElements.insert( make_pair(*type1,*type2) );
|
---|
[2fe971] | 158 | DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << *(*type1) << " and " << *(*type2) << "." << endl);
|
---|
[c78d44] | 159 | }
|
---|
| 160 | } else if (elements.size() == 1) { // one to all are valid
|
---|
[e5c0a1] | 161 | const element *elemental = *elements.begin();
|
---|
| 162 | PairsOfElements.insert( pair<const element *,const element*>(elemental,0) );
|
---|
| 163 | PairsOfElements.insert( pair<const element *,const element*>(0,elemental) );
|
---|
[c78d44] | 164 | } else { // all elements valid
|
---|
| 165 | PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
|
---|
| 166 | }
|
---|
| 167 |
|
---|
[c4d4df] | 168 | outmap = new PairCorrelationMap;
|
---|
[e65de8] | 169 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
|
---|
| 170 | DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
|
---|
| 171 | for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
|
---|
| 172 | DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
|
---|
| 173 | for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
|
---|
| 174 | DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
|
---|
| 175 | for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
|
---|
| 176 | DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
|
---|
| 177 | if ((*iter)->getId() < (*runner)->getId()){
|
---|
[b5c53d] | 178 | for (set <pair<const element *, const element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
|
---|
[d74077] | 179 | if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
|
---|
| 180 | distance = domain.periodicDistance((*iter)->getPosition(),(*runner)->getPosition());
|
---|
[e65de8] | 181 | //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
|
---|
| 182 | outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
|
---|
[a5551b] | 183 | }
|
---|
[c4d4df] | 184 | }
|
---|
[a5551b] | 185 | }
|
---|
[c4d4df] | 186 | }
|
---|
| 187 | }
|
---|
[24725c] | 188 | }
|
---|
[c4d4df] | 189 | return outmap;
|
---|
| 190 | };
|
---|
| 191 |
|
---|
[7ea9e6] | 192 | /** Calculates the pair correlation between given elements.
|
---|
| 193 | * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
|
---|
| 194 | * \param *molecules list of molecules structure
|
---|
[c78d44] | 195 | * \param &elements vector of elements to correlate
|
---|
[7ea9e6] | 196 | * \param ranges[NDIM] interval boundaries for the periodic images to scan also
|
---|
| 197 | * \return Map of doubles with values the pair of the two atoms.
|
---|
| 198 | */
|
---|
[e5c0a1] | 199 | PairCorrelationMap *PeriodicPairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const int ranges[NDIM] )
|
---|
[7ea9e6] | 200 | {
|
---|
[3930eb] | 201 | Info FunctionInfo(__func__);
|
---|
[caa30b] | 202 | PairCorrelationMap *outmap = new PairCorrelationMap;
|
---|
[7ea9e6] | 203 | double distance = 0.;
|
---|
| 204 | int n[NDIM];
|
---|
| 205 | Vector checkX;
|
---|
| 206 | Vector periodicX;
|
---|
| 207 | int Othern[NDIM];
|
---|
| 208 | Vector checkOtherX;
|
---|
| 209 | Vector periodicOtherX;
|
---|
| 210 |
|
---|
[e65de8] | 211 | if (molecules.empty()) {
|
---|
[58ed4a] | 212 | DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
|
---|
[7ea9e6] | 213 | return outmap;
|
---|
| 214 | }
|
---|
[e65de8] | 215 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
|
---|
[009607e] | 216 | (*MolWalker)->doCountAtoms();
|
---|
[c78d44] | 217 |
|
---|
| 218 | // create all possible pairs of elements
|
---|
[e5c0a1] | 219 | set <pair<const element *,const element *> > PairsOfElements;
|
---|
[c78d44] | 220 | if (elements.size() >= 2) {
|
---|
[e5c0a1] | 221 | for (vector<const element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
|
---|
| 222 | for (vector<const element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
|
---|
[c78d44] | 223 | if (type1 != type2) {
|
---|
[e5c0a1] | 224 | PairsOfElements.insert( make_pair(*type1,*type2) );
|
---|
[2fe971] | 225 | DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << *(*type1) << " and " << *(*type2) << "." << endl);
|
---|
[c78d44] | 226 | }
|
---|
| 227 | } else if (elements.size() == 1) { // one to all are valid
|
---|
[e5c0a1] | 228 | const element *elemental = *elements.begin();
|
---|
| 229 | PairsOfElements.insert( pair<const element *,const element*>(elemental,0) );
|
---|
| 230 | PairsOfElements.insert( pair<const element *,const element*>(0,elemental) );
|
---|
[c78d44] | 231 | } else { // all elements valid
|
---|
| 232 | PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
|
---|
| 233 | }
|
---|
| 234 |
|
---|
[7ea9e6] | 235 | outmap = new PairCorrelationMap;
|
---|
[e65de8] | 236 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
|
---|
[cca9ef] | 237 | RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
|
---|
| 238 | RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
|
---|
[e65de8] | 239 | DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
|
---|
| 240 | for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
|
---|
| 241 | DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
|
---|
[d74077] | 242 | periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
|
---|
[e65de8] | 243 | // go through every range in xyz and get distance
|
---|
| 244 | for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
|
---|
| 245 | for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
|
---|
| 246 | for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
|
---|
| 247 | checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
|
---|
| 248 | for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
|
---|
| 249 | DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
|
---|
| 250 | for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
|
---|
| 251 | DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
|
---|
| 252 | if ((*iter)->getId() < (*runner)->getId()){
|
---|
[e5c0a1] | 253 | for (set <pair<const element *,const element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
|
---|
[d74077] | 254 | if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
|
---|
| 255 | periodicOtherX = FullInverseMatrix * ((*runner)->getPosition()); // x now in [0,1)^3
|
---|
[e65de8] | 256 | // go through every range in xyz and get distance
|
---|
| 257 | for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
|
---|
| 258 | for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
|
---|
| 259 | for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
|
---|
| 260 | checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
|
---|
| 261 | distance = checkX.distance(checkOtherX);
|
---|
| 262 | //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
|
---|
| 263 | outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
|
---|
| 264 | }
|
---|
| 265 | }
|
---|
[c78d44] | 266 | }
|
---|
[7ea9e6] | 267 | }
|
---|
[c78d44] | 268 | }
|
---|
[7ea9e6] | 269 | }
|
---|
| 270 | }
|
---|
[c78d44] | 271 | }
|
---|
[7ea9e6] | 272 |
|
---|
| 273 | return outmap;
|
---|
| 274 | };
|
---|
| 275 |
|
---|
[c4d4df] | 276 | /** Calculates the distance (pair) correlation between a given element and a point.
|
---|
[a5551b] | 277 | * \param *molecules list of molecules structure
|
---|
[c78d44] | 278 | * \param &elements vector of elements to correlate with point
|
---|
[c4d4df] | 279 | * \param *point vector to the correlation point
|
---|
| 280 | * \return Map of dobules with values as pairs of atom and the vector
|
---|
| 281 | */
|
---|
[e5c0a1] | 282 | CorrelationToPointMap *CorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point )
|
---|
[c4d4df] | 283 | {
|
---|
[3930eb] | 284 | Info FunctionInfo(__func__);
|
---|
[caa30b] | 285 | CorrelationToPointMap *outmap = new CorrelationToPointMap;
|
---|
[c4d4df] | 286 | double distance = 0.;
|
---|
[014475] | 287 | Box &domain = World::getInstance().getDomain();
|
---|
[c4d4df] | 288 |
|
---|
[e65de8] | 289 | if (molecules.empty()) {
|
---|
[a67d19] | 290 | DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
|
---|
[c4d4df] | 291 | return outmap;
|
---|
| 292 | }
|
---|
[e65de8] | 293 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
|
---|
[009607e] | 294 | (*MolWalker)->doCountAtoms();
|
---|
[c4d4df] | 295 | outmap = new CorrelationToPointMap;
|
---|
[e65de8] | 296 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
|
---|
| 297 | DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
|
---|
| 298 | for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
|
---|
| 299 | DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
|
---|
[e5c0a1] | 300 | for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
|
---|
[d74077] | 301 | if ((*type == NULL) || ((*iter)->getType() == *type)) {
|
---|
| 302 | distance = domain.periodicDistance((*iter)->getPosition(),*point);
|
---|
[e65de8] | 303 | DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
|
---|
| 304 | outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
|
---|
| 305 | }
|
---|
[c4d4df] | 306 | }
|
---|
[e65de8] | 307 | }
|
---|
[c4d4df] | 308 |
|
---|
| 309 | return outmap;
|
---|
| 310 | };
|
---|
| 311 |
|
---|
[7ea9e6] | 312 | /** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
|
---|
| 313 | * \param *molecules list of molecules structure
|
---|
[c78d44] | 314 | * \param &elements vector of elements to correlate to point
|
---|
[7ea9e6] | 315 | * \param *point vector to the correlation point
|
---|
| 316 | * \param ranges[NDIM] interval boundaries for the periodic images to scan also
|
---|
| 317 | * \return Map of dobules with values as pairs of atom and the vector
|
---|
| 318 | */
|
---|
[e5c0a1] | 319 | CorrelationToPointMap *PeriodicCorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point, const int ranges[NDIM] )
|
---|
[7ea9e6] | 320 | {
|
---|
[3930eb] | 321 | Info FunctionInfo(__func__);
|
---|
[caa30b] | 322 | CorrelationToPointMap *outmap = new CorrelationToPointMap;
|
---|
[7ea9e6] | 323 | double distance = 0.;
|
---|
| 324 | int n[NDIM];
|
---|
| 325 | Vector periodicX;
|
---|
| 326 | Vector checkX;
|
---|
| 327 |
|
---|
[e65de8] | 328 | if (molecules.empty()) {
|
---|
[a67d19] | 329 | DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
|
---|
[7ea9e6] | 330 | return outmap;
|
---|
| 331 | }
|
---|
[e65de8] | 332 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
|
---|
[009607e] | 333 | (*MolWalker)->doCountAtoms();
|
---|
[7ea9e6] | 334 | outmap = new CorrelationToPointMap;
|
---|
[e65de8] | 335 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
|
---|
[cca9ef] | 336 | RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
|
---|
| 337 | RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
|
---|
[e65de8] | 338 | DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
|
---|
| 339 | for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
|
---|
| 340 | DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
|
---|
[e5c0a1] | 341 | for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
|
---|
[d74077] | 342 | if ((*type == NULL) || ((*iter)->getType() == *type)) {
|
---|
| 343 | periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
|
---|
[e65de8] | 344 | // go through every range in xyz and get distance
|
---|
| 345 | for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
|
---|
| 346 | for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
|
---|
| 347 | for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
|
---|
| 348 | checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
|
---|
| 349 | distance = checkX.distance(*point);
|
---|
| 350 | DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
|
---|
| 351 | outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
|
---|
| 352 | }
|
---|
| 353 | }
|
---|
[7ea9e6] | 354 | }
|
---|
[e65de8] | 355 | }
|
---|
[7ea9e6] | 356 |
|
---|
| 357 | return outmap;
|
---|
| 358 | };
|
---|
| 359 |
|
---|
[c4d4df] | 360 | /** Calculates the distance (pair) correlation between a given element and a surface.
|
---|
[a5551b] | 361 | * \param *molecules list of molecules structure
|
---|
[c78d44] | 362 | * \param &elements vector of elements to correlate to surface
|
---|
[c4d4df] | 363 | * \param *Surface pointer to Tesselation class surface
|
---|
| 364 | * \param *LC LinkedCell structure to quickly find neighbouring atoms
|
---|
| 365 | * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
|
---|
| 366 | */
|
---|
[e5c0a1] | 367 | CorrelationToSurfaceMap *CorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell *LC )
|
---|
[c4d4df] | 368 | {
|
---|
[3930eb] | 369 | Info FunctionInfo(__func__);
|
---|
[caa30b] | 370 | CorrelationToSurfaceMap *outmap = new CorrelationToSurfaceMap;
|
---|
[99593f] | 371 | double distance = 0;
|
---|
[c4d4df] | 372 | class BoundaryTriangleSet *triangle = NULL;
|
---|
| 373 | Vector centroid;
|
---|
[7ea9e6] | 374 |
|
---|
[e65de8] | 375 | if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
|
---|
[58ed4a] | 376 | DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
|
---|
[7ea9e6] | 377 | return outmap;
|
---|
| 378 | }
|
---|
[e65de8] | 379 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
|
---|
[009607e] | 380 | (*MolWalker)->doCountAtoms();
|
---|
[7ea9e6] | 381 | outmap = new CorrelationToSurfaceMap;
|
---|
[e65de8] | 382 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
|
---|
| 383 | DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << (*MolWalker)->name << "." << endl);
|
---|
| 384 | if ((*MolWalker)->empty())
|
---|
| 385 | DoLog(2) && (2) && (Log() << Verbose(2) << "\t is empty." << endl);
|
---|
| 386 | for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
|
---|
| 387 | DoLog(3) && (Log() << Verbose(3) << "\tCurrent atom is " << *(*iter) << "." << endl);
|
---|
[e5c0a1] | 388 | for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
|
---|
[d74077] | 389 | if ((*type == NULL) || ((*iter)->getType() == *type)) {
|
---|
| 390 | TriangleIntersectionList Intersections((*iter)->getPosition(),Surface,LC);
|
---|
[e65de8] | 391 | distance = Intersections.GetSmallestDistance();
|
---|
| 392 | triangle = Intersections.GetClosestTriangle();
|
---|
| 393 | outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
|
---|
| 394 | }
|
---|
[7fd416] | 395 | }
|
---|
[e65de8] | 396 | }
|
---|
[7ea9e6] | 397 |
|
---|
| 398 | return outmap;
|
---|
| 399 | };
|
---|
| 400 |
|
---|
| 401 | /** Calculates the distance (pair) correlation between a given element, all its periodic images and and a surface.
|
---|
| 402 | * Note that we also put all periodic images found in the cells given by [ -ranges[i], ranges[i] ] and i=0,...,NDIM-1.
|
---|
| 403 | * I.e. We multiply the atom::node with the inverse of the domain matrix, i.e. transform it to \f$[0,0^3\f$, then add per
|
---|
| 404 | * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
|
---|
| 405 | * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
|
---|
| 406 | * \param *molecules list of molecules structure
|
---|
[c78d44] | 407 | * \param &elements vector of elements to correlate to surface
|
---|
[7ea9e6] | 408 | * \param *Surface pointer to Tesselation class surface
|
---|
| 409 | * \param *LC LinkedCell structure to quickly find neighbouring atoms
|
---|
| 410 | * \param ranges[NDIM] interval boundaries for the periodic images to scan also
|
---|
| 411 | * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
|
---|
| 412 | */
|
---|
[e5c0a1] | 413 | CorrelationToSurfaceMap *PeriodicCorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
|
---|
[7ea9e6] | 414 | {
|
---|
[3930eb] | 415 | Info FunctionInfo(__func__);
|
---|
[caa30b] | 416 | CorrelationToSurfaceMap *outmap = new CorrelationToSurfaceMap;
|
---|
[7ea9e6] | 417 | double distance = 0;
|
---|
| 418 | class BoundaryTriangleSet *triangle = NULL;
|
---|
| 419 | Vector centroid;
|
---|
[99593f] | 420 | int n[NDIM];
|
---|
| 421 | Vector periodicX;
|
---|
| 422 | Vector checkX;
|
---|
[c4d4df] | 423 |
|
---|
[e65de8] | 424 | if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
|
---|
[a67d19] | 425 | DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
|
---|
[c4d4df] | 426 | return outmap;
|
---|
| 427 | }
|
---|
[e65de8] | 428 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
|
---|
[009607e] | 429 | (*MolWalker)->doCountAtoms();
|
---|
[c4d4df] | 430 | outmap = new CorrelationToSurfaceMap;
|
---|
[244a84] | 431 | double ShortestDistance = 0.;
|
---|
| 432 | BoundaryTriangleSet *ShortestTriangle = NULL;
|
---|
[e65de8] | 433 | for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
|
---|
[cca9ef] | 434 | RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
|
---|
| 435 | RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
|
---|
[e65de8] | 436 | DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
|
---|
| 437 | for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
|
---|
| 438 | DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
|
---|
[e5c0a1] | 439 | for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
|
---|
[d74077] | 440 | if ((*type == NULL) || ((*iter)->getType() == *type)) {
|
---|
| 441 | periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
|
---|
[e65de8] | 442 | // go through every range in xyz and get distance
|
---|
| 443 | ShortestDistance = -1.;
|
---|
| 444 | for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
|
---|
| 445 | for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
|
---|
| 446 | for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
|
---|
| 447 | checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
|
---|
[d74077] | 448 | TriangleIntersectionList Intersections(checkX,Surface,LC);
|
---|
[e65de8] | 449 | distance = Intersections.GetSmallestDistance();
|
---|
| 450 | triangle = Intersections.GetClosestTriangle();
|
---|
| 451 | if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
|
---|
| 452 | ShortestDistance = distance;
|
---|
| 453 | ShortestTriangle = triangle;
|
---|
[99593f] | 454 | }
|
---|
[e65de8] | 455 | }
|
---|
| 456 | // insert
|
---|
| 457 | outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
|
---|
| 458 | //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
|
---|
| 459 | }
|
---|
[c4d4df] | 460 | }
|
---|
[e65de8] | 461 | }
|
---|
[c4d4df] | 462 |
|
---|
| 463 | return outmap;
|
---|
| 464 | };
|
---|
| 465 |
|
---|
[bd61b41] | 466 | /** Returns the index of the bin for a given value.
|
---|
[c4d4df] | 467 | * \param value value whose bin to look for
|
---|
| 468 | * \param BinWidth width of bin
|
---|
| 469 | * \param BinStart first bin
|
---|
| 470 | */
|
---|
[bd61b41] | 471 | int GetBin ( const double value, const double BinWidth, const double BinStart )
|
---|
[c4d4df] | 472 | {
|
---|
[92e5cb] | 473 | //Info FunctionInfo(__func__);
|
---|
[bd61b41] | 474 | int bin =(int) (floor((value - BinStart)/BinWidth));
|
---|
| 475 | return (bin);
|
---|
[c4d4df] | 476 | };
|
---|
| 477 |
|
---|
| 478 |
|
---|
[92e5cb] | 479 | /** Adds header part that is unique to BinPairMap.
|
---|
| 480 | *
|
---|
| 481 | * @param file stream to print to
|
---|
[c4d4df] | 482 | */
|
---|
[92e5cb] | 483 | void OutputCorrelation_Header( ofstream * const file )
|
---|
[c4d4df] | 484 | {
|
---|
[92e5cb] | 485 | *file << "\tCount";
|
---|
[c4d4df] | 486 | };
|
---|
[b1f254] | 487 |
|
---|
[92e5cb] | 488 | /** Prints values stored in BinPairMap iterator.
|
---|
| 489 | *
|
---|
| 490 | * @param file stream to print to
|
---|
| 491 | * @param runner iterator pointing at values to print
|
---|
[be945c] | 492 | */
|
---|
[92e5cb] | 493 | void OutputCorrelation_Value( ofstream * const file, BinPairMap::const_iterator &runner )
|
---|
[be945c] | 494 | {
|
---|
[92e5cb] | 495 | *file << runner->second;
|
---|
[be945c] | 496 | };
|
---|
| 497 |
|
---|
[92e5cb] | 498 |
|
---|
| 499 | /** Adds header part that is unique to DipoleAngularCorrelationMap.
|
---|
| 500 | *
|
---|
| 501 | * @param file stream to print to
|
---|
[b1f254] | 502 | */
|
---|
[92e5cb] | 503 | void OutputDipoleAngularCorrelation_Header( ofstream * const file )
|
---|
[b1f254] | 504 | {
|
---|
[92e5cb] | 505 | *file << "\tAtom1\tAtom2";
|
---|
[b1f254] | 506 | };
|
---|
| 507 |
|
---|
[92e5cb] | 508 | /** Prints values stored in DipoleAngularCorrelationMap iterator.
|
---|
| 509 | *
|
---|
| 510 | * @param file stream to print to
|
---|
| 511 | * @param runner iterator pointing at values to print
|
---|
[b1f254] | 512 | */
|
---|
[92e5cb] | 513 | void OutputDipoleAngularCorrelation_Value( ofstream * const file, DipoleAngularCorrelationMap::const_iterator &runner )
|
---|
[b1f254] | 514 | {
|
---|
[92e5cb] | 515 | *file << runner->second.first->getId() << "\t" << runner->second.second->getId();
|
---|
[b1f254] | 516 | };
|
---|
| 517 |
|
---|
[92e5cb] | 518 |
|
---|
| 519 | /** Adds header part that is unique to PairCorrelationMap.
|
---|
| 520 | *
|
---|
| 521 | * @param file stream to print to
|
---|
[b1f254] | 522 | */
|
---|
[92e5cb] | 523 | void OutputPairCorrelation_Header( ofstream * const file )
|
---|
[b1f254] | 524 | {
|
---|
[92e5cb] | 525 | *file << "\tAtom1\tAtom2";
|
---|
| 526 | };
|
---|
| 527 |
|
---|
| 528 | /** Prints values stored in PairCorrelationMap iterator.
|
---|
| 529 | *
|
---|
| 530 | * @param file stream to print to
|
---|
| 531 | * @param runner iterator pointing at values to print
|
---|
| 532 | */
|
---|
| 533 | void OutputPairCorrelation_Value( ofstream * const file, PairCorrelationMap::const_iterator &runner )
|
---|
| 534 | {
|
---|
| 535 | *file << *(runner->second.first) << "\t" << *(runner->second.second);
|
---|
| 536 | };
|
---|
| 537 |
|
---|
| 538 |
|
---|
| 539 | /** Adds header part that is unique to CorrelationToPointMap.
|
---|
| 540 | *
|
---|
| 541 | * @param file stream to print to
|
---|
| 542 | */
|
---|
| 543 | void OutputCorrelationToPoint_Header( ofstream * const file )
|
---|
| 544 | {
|
---|
| 545 | *file << "\tAtom::x[i]-point.x[i]";
|
---|
| 546 | };
|
---|
| 547 |
|
---|
| 548 | /** Prints values stored in CorrelationToPointMap iterator.
|
---|
| 549 | *
|
---|
| 550 | * @param file stream to print to
|
---|
| 551 | * @param runner iterator pointing at values to print
|
---|
| 552 | */
|
---|
| 553 | void OutputCorrelationToPoint_Value( ofstream * const file, CorrelationToPointMap::const_iterator &runner )
|
---|
| 554 | {
|
---|
| 555 | for (int i=0;i<NDIM;i++)
|
---|
| 556 | *file << "\t" << setprecision(8) << (runner->second.first->at(i) - runner->second.second->at(i));
|
---|
[b1f254] | 557 | };
|
---|
| 558 |
|
---|
[92e5cb] | 559 |
|
---|
| 560 | /** Adds header part that is unique to CorrelationToSurfaceMap.
|
---|
| 561 | *
|
---|
| 562 | * @param file stream to print to
|
---|
| 563 | */
|
---|
| 564 | void OutputCorrelationToSurface_Header( ofstream * const file )
|
---|
| 565 | {
|
---|
| 566 | *file << "\tTriangle";
|
---|
| 567 | };
|
---|
| 568 |
|
---|
| 569 | /** Prints values stored in CorrelationToSurfaceMap iterator.
|
---|
| 570 | *
|
---|
| 571 | * @param file stream to print to
|
---|
| 572 | * @param runner iterator pointing at values to print
|
---|
| 573 | */
|
---|
| 574 | void OutputCorrelationToSurface_Value( ofstream * const file, CorrelationToSurfaceMap::const_iterator &runner )
|
---|
| 575 | {
|
---|
| 576 | *file << *(runner->second.first) << "\t" << *(runner->second.second);
|
---|
| 577 | };
|
---|