[96c961] | 1 | /*
|
---|
| 2 | * analysis_bonds.cpp
|
---|
| 3 | *
|
---|
| 4 | * Created on: Nov 7, 2009
|
---|
| 5 | * Author: heber
|
---|
| 6 | */
|
---|
| 7 |
|
---|
[112b09] | 8 | #include "Helpers/MemDebug.hpp"
|
---|
| 9 |
|
---|
[220cf37] | 10 | #include "analysis_bonds.hpp"
|
---|
| 11 | #include "atom.hpp"
|
---|
| 12 | #include "bond.hpp"
|
---|
[388049] | 13 | #include "element.hpp"
|
---|
| 14 | #include "info.hpp"
|
---|
[36166d] | 15 | #include "verbose.hpp"
|
---|
[220cf37] | 16 | #include "log.hpp"
|
---|
| 17 | #include "molecule.hpp"
|
---|
| 18 |
|
---|
| 19 | /** Calculates the min, mean and maximum bond counts for the given molecule.
|
---|
| 20 | * \param *mol molecule with atoms and atom::ListOfBonds
|
---|
| 21 | * \param &Min minimum count on return
|
---|
| 22 | * \param &Mean mean count on return
|
---|
| 23 | * \param &Max maximum count on return
|
---|
| 24 | */
|
---|
| 25 | void GetMaxMinMeanBondCount(const molecule * const mol, double &Min, double &Mean, double &Max)
|
---|
| 26 | {
|
---|
| 27 | Min = 2e+6;
|
---|
| 28 | Max = -2e+5;
|
---|
| 29 | Mean = 0.;
|
---|
| 30 |
|
---|
| 31 | int AtomCount = 0;
|
---|
[9879f6] | 32 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
| 33 | const int count = (*iter)->ListOfBonds.size();
|
---|
[220cf37] | 34 | if (Max < count)
|
---|
| 35 | Max = count;
|
---|
| 36 | if (Min > count)
|
---|
| 37 | Min = count;
|
---|
| 38 | Mean += count;
|
---|
| 39 | AtomCount++;
|
---|
| 40 | }
|
---|
| 41 | if (((int)Mean % 2) != 0)
|
---|
[58ed4a] | 42 | DoeLog(1) && (eLog()<< Verbose(1) << "Something is wrong with the bond structure, the number of bonds is not even!" << endl);
|
---|
[220cf37] | 43 | Mean /= (double)AtomCount;
|
---|
| 44 | };
|
---|
| 45 |
|
---|
| 46 | /** Calculates the min and max bond distance of all atoms of two given elements.
|
---|
| 47 | * \param *mol molecule with atoms
|
---|
| 48 | * \param *type1 one element
|
---|
| 49 | * \param *type2 other element
|
---|
| 50 | * \param &Min minimum distance on return, 0 if no bond between the two elements
|
---|
| 51 | * \param &Mean mean distance (i.e. sum of distance for matching element pairs, divided by number) on return, 0 if no bond between the two elements
|
---|
| 52 | * \param &Max maximum distance on return, 0 if no bond between the two elements
|
---|
| 53 | */
|
---|
[4eb4fe] | 54 | void MinMeanMaxBondDistanceBetweenElements(const molecule *mol, const element *type1, const element *type2, double &Min, double &Mean, double &Max)
|
---|
[220cf37] | 55 | {
|
---|
| 56 | Min = 2e+6;
|
---|
| 57 | Mean = 0.;
|
---|
| 58 | Max = -2e+6;
|
---|
| 59 |
|
---|
| 60 | int AtomNo = 0;
|
---|
[9879f6] | 61 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
| 62 | if ((*iter)->type == type1)
|
---|
| 63 | for (BondList::const_iterator BondRunner = (*iter)->ListOfBonds.begin(); BondRunner != (*iter)->ListOfBonds.end(); BondRunner++)
|
---|
| 64 | if ((*BondRunner)->GetOtherAtom((*iter))->type == type2) {
|
---|
[220cf37] | 65 | const double distance = (*BondRunner)->GetDistanceSquared();
|
---|
| 66 | if (Min > distance)
|
---|
| 67 | Min = distance;
|
---|
| 68 | if (Max < distance)
|
---|
| 69 | Max = distance;
|
---|
| 70 | Mean += sqrt(distance);
|
---|
| 71 | AtomNo++;
|
---|
| 72 | }
|
---|
| 73 | }
|
---|
| 74 | if (Max < 0) {
|
---|
| 75 | Max = Min = 0.;
|
---|
| 76 | } else {
|
---|
| 77 | Max = sqrt(Max);
|
---|
| 78 | Min = sqrt(Min);
|
---|
| 79 | Mean = Mean/(double)AtomNo;
|
---|
| 80 | }
|
---|
| 81 | };
|
---|
[388049] | 82 |
|
---|
[fe238c] | 83 | /** Calculate the angle between \a *first and \a *origin and \a *second and \a *origin.
|
---|
| 84 | * \param *first first Vector
|
---|
| 85 | * \param *origin origin of angle taking
|
---|
| 86 | * \param *second second Vector
|
---|
| 87 | * \return angle between \a *first and \a *second, both relative to origin at \a *origin.
|
---|
| 88 | */
|
---|
| 89 | double CalculateAngle(Vector *first, Vector *central, Vector *second)
|
---|
| 90 | {
|
---|
| 91 | Vector OHBond;
|
---|
| 92 | Vector OOBond;
|
---|
| 93 |
|
---|
[8cbb97] | 94 | OHBond = (*first) - (*central);
|
---|
| 95 | OOBond = (*second) - (*central);
|
---|
| 96 | const double angle = OHBond.Angle(OOBond);
|
---|
[fe238c] | 97 | return angle;
|
---|
| 98 | };
|
---|
| 99 |
|
---|
| 100 | /** Checks whether the angle between \a *Oxygen and \a *Hydrogen and \a *Oxygen and \a *OtherOxygen is less than 30 degrees.
|
---|
| 101 | * Note that distance criterion is not checked.
|
---|
| 102 | * \param *Oxygen first oxygen atom, bonded to \a *Hydrogen
|
---|
| 103 | * \param *Hydrogen hydrogen bonded to \a *Oxygen
|
---|
| 104 | * \param *OtherOxygen other oxygen atom
|
---|
| 105 | * \return true - angle criteria fulfilled, false - criteria not fulfilled, angle greater than 30 degrees.
|
---|
| 106 | */
|
---|
| 107 | bool CheckHydrogenBridgeBondAngle(atom *Oxygen, atom *Hydrogen, atom *OtherOxygen)
|
---|
| 108 | {
|
---|
| 109 | Info FunctionInfo(__func__);
|
---|
| 110 |
|
---|
| 111 | // check angle
|
---|
| 112 | if (CalculateAngle(&Hydrogen->x, &Oxygen->x, &OtherOxygen->x) < M_PI*(30./180.)) {
|
---|
| 113 | return true;
|
---|
| 114 | } else {
|
---|
| 115 | return false;
|
---|
| 116 | }
|
---|
| 117 | };
|
---|
[388049] | 118 |
|
---|
| 119 | /** Counts the number of hydrogen bridge bonds.
|
---|
| 120 | * With \a *InterfaceElement an extra element can be specified that identifies some boundary.
|
---|
| 121 | * Then, counting is for the h-bridges that connect to interface only.
|
---|
| 122 | * \param *molecules molecules to count bonds
|
---|
| 123 | * \param *InterfaceElement or NULL
|
---|
[bfd839] | 124 | * \param *Interface2Element or NULL
|
---|
[388049] | 125 | */
|
---|
[bfd839] | 126 | int CountHydrogenBridgeBonds(MoleculeListClass *molecules, const element * InterfaceElement = NULL, const element * Interface2Element = NULL)
|
---|
[388049] | 127 | {
|
---|
| 128 | int count = 0;
|
---|
[fe238c] | 129 | int OtherHydrogens = 0;
|
---|
| 130 | double Otherangle = 0.;
|
---|
[388049] | 131 | bool InterfaceFlag = false;
|
---|
[bfd839] | 132 | bool Interface2Flag = false;
|
---|
[fe238c] | 133 | bool OtherHydrogenFlag = true;
|
---|
[a7b761b] | 134 | for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); ++MolWalker) {
|
---|
| 135 | molecule::iterator Walker = (*MolWalker)->begin();
|
---|
| 136 | for(;Walker!=(*MolWalker)->end();++Walker){
|
---|
| 137 | for (MoleculeList::const_iterator MolRunner = molecules->ListOfMolecules.begin();MolRunner != molecules->ListOfMolecules.end(); ++MolRunner) {
|
---|
| 138 | molecule::iterator Runner = (*MolRunner)->begin();
|
---|
| 139 | for(;Runner!=(*MolRunner)->end();++Runner){
|
---|
| 140 | if (((*Walker)->type->Z == 8) && ((*Runner)->type->Z == 8)) {
|
---|
[388049] | 141 | // check distance
|
---|
[a7b761b] | 142 | const double distance = (*Runner)->x.DistanceSquared((*Walker)->x);
|
---|
[388049] | 143 | if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) { // distance >0 means different atoms
|
---|
[fe238c] | 144 | // on other atom(Runner) we check for bond to interface element and
|
---|
| 145 | // check that O-O line is not in between the shanks of the two connected hydrogens (Otherangle > 104.5)
|
---|
| 146 | OtherHydrogenFlag = true;
|
---|
| 147 | Otherangle = 0.;
|
---|
| 148 | OtherHydrogens = 0;
|
---|
[388049] | 149 | InterfaceFlag = (InterfaceElement == NULL);
|
---|
[bfd839] | 150 | Interface2Flag = (Interface2Element == NULL);
|
---|
[a7b761b] | 151 | for (BondList::const_iterator BondRunner = (*Runner)->ListOfBonds.begin(); BondRunner != (*Runner)->ListOfBonds.end(); BondRunner++) {
|
---|
| 152 | atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Runner);
|
---|
[fe238c] | 153 | // if hydrogen, check angle to be greater(!) than 30 degrees
|
---|
| 154 | if (OtherAtom->type->Z == 1) {
|
---|
[a7b761b] | 155 | const double angle = CalculateAngle(&OtherAtom->x, &(*Runner)->x, &(*Walker)->x);
|
---|
[fe238c] | 156 | OtherHydrogenFlag = OtherHydrogenFlag && (angle > M_PI*(30./180.) + MYEPSILON);
|
---|
| 157 | Otherangle += angle;
|
---|
| 158 | OtherHydrogens++;
|
---|
| 159 | }
|
---|
[388049] | 160 | InterfaceFlag = InterfaceFlag || (OtherAtom->type == InterfaceElement);
|
---|
[bfd839] | 161 | Interface2Flag = Interface2Flag || (OtherAtom->type == Interface2Element);
|
---|
[388049] | 162 | }
|
---|
[fe238c] | 163 | DoLog(1) && (Log() << Verbose(1) << "Otherangle is " << Otherangle << " for " << OtherHydrogens << " hydrogens." << endl);
|
---|
| 164 | switch (OtherHydrogens) {
|
---|
| 165 | case 0:
|
---|
| 166 | case 1:
|
---|
| 167 | break;
|
---|
| 168 | case 2:
|
---|
| 169 | OtherHydrogenFlag = OtherHydrogenFlag && (Otherangle > M_PI*(104.5/180.) + MYEPSILON);
|
---|
| 170 | break;
|
---|
| 171 | default: // 3 or more hydrogens ...
|
---|
| 172 | OtherHydrogenFlag = false;
|
---|
| 173 | break;
|
---|
| 174 | }
|
---|
[bfd839] | 175 | if (InterfaceFlag && Interface2Flag && OtherHydrogenFlag) {
|
---|
[388049] | 176 | // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
|
---|
[a7b761b] | 177 | for (BondList::const_iterator BondRunner = (*Walker)->ListOfBonds.begin(); BondRunner != (*Walker)->ListOfBonds.end(); BondRunner++) {
|
---|
| 178 | atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Walker);
|
---|
[388049] | 179 | if (OtherAtom->type->Z == 1) {
|
---|
| 180 | // check angle
|
---|
[a7b761b] | 181 | if (CheckHydrogenBridgeBondAngle(*Walker, OtherAtom, *Runner)) {
|
---|
| 182 | DoLog(1) && (Log() << Verbose(1) << (*Walker)->getName() << ", " << OtherAtom->getName() << " and " << (*Runner)->getName() << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(&OtherAtom->x, &(*Walker)->x, &(*Runner)->x)*(180./M_PI) << "." << endl);
|
---|
[388049] | 183 | count++;
|
---|
| 184 | break;
|
---|
| 185 | }
|
---|
| 186 | }
|
---|
| 187 | }
|
---|
| 188 | }
|
---|
| 189 | }
|
---|
| 190 | }
|
---|
| 191 | }
|
---|
| 192 | }
|
---|
| 193 | }
|
---|
| 194 | }
|
---|
| 195 | return count;
|
---|
| 196 | }
|
---|
| 197 |
|
---|
| 198 | /** Counts the number of bonds between two given elements.
|
---|
| 199 | * \param *molecules list of molecules with all atoms
|
---|
| 200 | * \param *first pointer to first element
|
---|
| 201 | * \param *second pointer to second element
|
---|
| 202 | * \return number of found bonds (\a *first-\a *second)
|
---|
| 203 | */
|
---|
| 204 | int CountBondsOfTwo(MoleculeListClass * const molecules, const element * const first, const element * const second)
|
---|
| 205 | {
|
---|
| 206 | int count = 0;
|
---|
| 207 |
|
---|
| 208 | for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
|
---|
[a7b761b] | 209 | molecule::iterator Walker = (*MolWalker)->begin();
|
---|
| 210 | for(;Walker!=(*MolWalker)->end();++Walker){
|
---|
| 211 | atom * theAtom = *Walker;
|
---|
| 212 | if ((theAtom->type == first) || (theAtom->type == second)) { // first element matches
|
---|
| 213 | for (BondList::const_iterator BondRunner = theAtom->ListOfBonds.begin(); BondRunner != theAtom->ListOfBonds.end(); BondRunner++) {
|
---|
| 214 | atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
|
---|
| 215 | if (((OtherAtom->type == first) || (OtherAtom->type == second)) && (theAtom->nr < OtherAtom->nr)) {
|
---|
[388049] | 216 | count++;
|
---|
| 217 | DoLog(1) && (Log() << Verbose(1) << first->name << "-" << second->name << " bond found between " << *Walker << " and " << *OtherAtom << "." << endl);
|
---|
| 218 | }
|
---|
| 219 | }
|
---|
| 220 | }
|
---|
| 221 | }
|
---|
| 222 | }
|
---|
| 223 | return count;
|
---|
| 224 | };
|
---|
| 225 |
|
---|
| 226 | /** Counts the number of bonds between three given elements.
|
---|
| 227 | * Note that we do not look for arbitrary sequence of given bonds, but \a *second will be the central atom and we check
|
---|
| 228 | * whether it has bonds to both \a *first and \a *third.
|
---|
| 229 | * \param *molecules list of molecules with all atoms
|
---|
| 230 | * \param *first pointer to first element
|
---|
| 231 | * \param *second pointer to second element
|
---|
| 232 | * \param *third pointer to third element
|
---|
| 233 | * \return number of found bonds (\a *first-\a *second-\a *third, \a *third-\a *second-\a *first, respectively)
|
---|
| 234 | */
|
---|
| 235 | int CountBondsOfThree(MoleculeListClass * const molecules, const element * const first, const element * const second, const element * const third)
|
---|
| 236 | {
|
---|
| 237 | int count = 0;
|
---|
| 238 | bool MatchFlag[2];
|
---|
| 239 | bool result = false;
|
---|
| 240 | const element * ElementArray[2];
|
---|
| 241 | ElementArray[0] = first;
|
---|
| 242 | ElementArray[1] = third;
|
---|
| 243 |
|
---|
| 244 | for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
|
---|
[a7b761b] | 245 | molecule::iterator Walker = (*MolWalker)->begin();
|
---|
| 246 | for(;Walker!=(*MolWalker)->end();++Walker){
|
---|
| 247 | atom *theAtom = *Walker;
|
---|
| 248 | if (theAtom->type == second) { // first element matches
|
---|
[388049] | 249 | for (int i=0;i<2;i++)
|
---|
| 250 | MatchFlag[i] = false;
|
---|
[a7b761b] | 251 | for (BondList::const_iterator BondRunner = theAtom->ListOfBonds.begin(); BondRunner != theAtom->ListOfBonds.end(); BondRunner++) {
|
---|
| 252 | atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
|
---|
[388049] | 253 | for (int i=0;i<2;i++)
|
---|
| 254 | if ((!MatchFlag[i]) && (OtherAtom->type == ElementArray[i])) {
|
---|
| 255 | MatchFlag[i] = true;
|
---|
| 256 | break; // each bonding atom can match at most one element we are looking for
|
---|
| 257 | }
|
---|
| 258 | }
|
---|
| 259 | result = true;
|
---|
| 260 | for (int i=0;i<2;i++) // gather results
|
---|
| 261 | result = result && MatchFlag[i];
|
---|
| 262 | if (result) { // check results
|
---|
| 263 | count++;
|
---|
| 264 | DoLog(1) && (Log() << Verbose(1) << first->name << "-" << second->name << "-" << third->name << " bond found at " << *Walker << "." << endl);
|
---|
| 265 | }
|
---|
| 266 | }
|
---|
| 267 | }
|
---|
| 268 | }
|
---|
| 269 | return count;
|
---|
| 270 | };
|
---|