Changes in src/analysis_bonds.cpp [220cf37:fe238c]
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/analysis_bonds.cpp
r220cf37 rfe238c 9 9 #include "atom.hpp" 10 10 #include "bond.hpp" 11 #include "element.hpp" 12 #include "info.hpp" 11 13 #include "log.hpp" 12 14 #include "molecule.hpp" … … 37 39 } 38 40 if (((int)Mean % 2) != 0) 39 eLog() << Verbose(1) << "Something is wrong with the bond structure, the number of bonds is not even!" << endl;41 DoeLog(1) && (eLog()<< Verbose(1) << "Something is wrong with the bond structure, the number of bonds is not even!" << endl); 40 42 Mean /= (double)AtomCount; 41 43 }; … … 79 81 } 80 82 }; 83 84 /** Calculate the angle between \a *first and \a *origin and \a *second and \a *origin. 85 * \param *first first Vector 86 * \param *origin origin of angle taking 87 * \param *second second Vector 88 * \return angle between \a *first and \a *second, both relative to origin at \a *origin. 89 */ 90 double CalculateAngle(Vector *first, Vector *central, Vector *second) 91 { 92 Vector OHBond; 93 Vector OOBond; 94 95 OHBond.CopyVector(first); 96 OHBond.SubtractVector(central); 97 OOBond.CopyVector(second); 98 OOBond.SubtractVector(central); 99 const double angle = OHBond.Angle(&OOBond); 100 return angle; 101 }; 102 103 /** Checks whether the angle between \a *Oxygen and \a *Hydrogen and \a *Oxygen and \a *OtherOxygen is less than 30 degrees. 104 * Note that distance criterion is not checked. 105 * \param *Oxygen first oxygen atom, bonded to \a *Hydrogen 106 * \param *Hydrogen hydrogen bonded to \a *Oxygen 107 * \param *OtherOxygen other oxygen atom 108 * \return true - angle criteria fulfilled, false - criteria not fulfilled, angle greater than 30 degrees. 109 */ 110 bool CheckHydrogenBridgeBondAngle(atom *Oxygen, atom *Hydrogen, atom *OtherOxygen) 111 { 112 Info FunctionInfo(__func__); 113 114 // check angle 115 if (CalculateAngle(&Hydrogen->x, &Oxygen->x, &OtherOxygen->x) < M_PI*(30./180.)) { 116 return true; 117 } else { 118 return false; 119 } 120 }; 121 122 /** Counts the number of hydrogen bridge bonds. 123 * With \a *InterfaceElement an extra element can be specified that identifies some boundary. 124 * Then, counting is for the h-bridges that connect to interface only. 125 * \param *molecules molecules to count bonds 126 * \param *InterfaceElement or NULL 127 */ 128 int CountHydrogenBridgeBonds(MoleculeListClass *molecules, element * InterfaceElement = NULL) 129 { 130 atom *Walker = NULL; 131 atom *Runner = NULL; 132 int count = 0; 133 int OtherHydrogens = 0; 134 double Otherangle = 0.; 135 bool InterfaceFlag = false; 136 bool OtherHydrogenFlag = true; 137 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) { 138 Walker = (*MolWalker)->start; 139 while (Walker->next != (*MolWalker)->end) { 140 Walker = Walker->next; 141 for (MoleculeList::const_iterator MolRunner = molecules->ListOfMolecules.begin();MolRunner != molecules->ListOfMolecules.end(); MolRunner++) { 142 Runner = (*MolRunner)->start; 143 while (Runner->next != (*MolRunner)->end) { 144 Runner = Runner->next; 145 if ((Walker->type->Z == 8) && (Runner->type->Z == 8)) { 146 // check distance 147 const double distance = Runner->x.DistanceSquared(&Walker->x); 148 if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) { // distance >0 means different atoms 149 // on other atom(Runner) we check for bond to interface element and 150 // check that O-O line is not in between the shanks of the two connected hydrogens (Otherangle > 104.5) 151 OtherHydrogenFlag = true; 152 Otherangle = 0.; 153 OtherHydrogens = 0; 154 InterfaceFlag = (InterfaceElement == NULL); 155 for (BondList::const_iterator BondRunner = Runner->ListOfBonds.begin(); BondRunner != Runner->ListOfBonds.end(); BondRunner++) { 156 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Runner); 157 // if hydrogen, check angle to be greater(!) than 30 degrees 158 if (OtherAtom->type->Z == 1) { 159 const double angle = CalculateAngle(&OtherAtom->x, &Runner->x, &Walker->x); 160 OtherHydrogenFlag = OtherHydrogenFlag && (angle > M_PI*(30./180.) + MYEPSILON); 161 Otherangle += angle; 162 OtherHydrogens++; 163 } 164 InterfaceFlag = InterfaceFlag || (OtherAtom->type == InterfaceElement); 165 } 166 DoLog(1) && (Log() << Verbose(1) << "Otherangle is " << Otherangle << " for " << OtherHydrogens << " hydrogens." << endl); 167 switch (OtherHydrogens) { 168 case 0: 169 case 1: 170 break; 171 case 2: 172 OtherHydrogenFlag = OtherHydrogenFlag && (Otherangle > M_PI*(104.5/180.) + MYEPSILON); 173 break; 174 default: // 3 or more hydrogens ... 175 OtherHydrogenFlag = false; 176 break; 177 } 178 if (InterfaceFlag && OtherHydrogenFlag) { 179 // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule 180 for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) { 181 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker); 182 if (OtherAtom->type->Z == 1) { 183 // check angle 184 if (CheckHydrogenBridgeBondAngle(Walker, OtherAtom, Runner)) { 185 DoLog(1) && (Log() << Verbose(1) << Walker->Name << ", " << OtherAtom->Name << " and " << Runner->Name << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(&OtherAtom->x, &Walker->x, &Runner->x)*(180./M_PI) << "." << endl); 186 count++; 187 break; 188 } 189 } 190 } 191 } 192 } 193 } 194 } 195 } 196 } 197 } 198 return count; 199 } 200 201 /** Counts the number of bonds between two given elements. 202 * \param *molecules list of molecules with all atoms 203 * \param *first pointer to first element 204 * \param *second pointer to second element 205 * \return number of found bonds (\a *first-\a *second) 206 */ 207 int CountBondsOfTwo(MoleculeListClass * const molecules, const element * const first, const element * const second) 208 { 209 atom *Walker = NULL; 210 int count = 0; 211 212 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) { 213 Walker = (*MolWalker)->start; 214 while (Walker->next != (*MolWalker)->end) { 215 Walker = Walker->next; 216 if ((Walker->type == first) || (Walker->type == second)) { // first element matches 217 for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) { 218 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker); 219 if (((OtherAtom->type == first) || (OtherAtom->type == second)) && (Walker->nr < OtherAtom->nr)) { 220 count++; 221 DoLog(1) && (Log() << Verbose(1) << first->name << "-" << second->name << " bond found between " << *Walker << " and " << *OtherAtom << "." << endl); 222 } 223 } 224 } 225 } 226 } 227 return count; 228 }; 229 230 /** Counts the number of bonds between three given elements. 231 * Note that we do not look for arbitrary sequence of given bonds, but \a *second will be the central atom and we check 232 * whether it has bonds to both \a *first and \a *third. 233 * \param *molecules list of molecules with all atoms 234 * \param *first pointer to first element 235 * \param *second pointer to second element 236 * \param *third pointer to third element 237 * \return number of found bonds (\a *first-\a *second-\a *third, \a *third-\a *second-\a *first, respectively) 238 */ 239 int CountBondsOfThree(MoleculeListClass * const molecules, const element * const first, const element * const second, const element * const third) 240 { 241 int count = 0; 242 bool MatchFlag[2]; 243 bool result = false; 244 atom *Walker = NULL; 245 const element * ElementArray[2]; 246 ElementArray[0] = first; 247 ElementArray[1] = third; 248 249 for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) { 250 Walker = (*MolWalker)->start; 251 while (Walker->next != (*MolWalker)->end) { 252 Walker = Walker->next; 253 if (Walker->type == second) { // first element matches 254 for (int i=0;i<2;i++) 255 MatchFlag[i] = false; 256 for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) { 257 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker); 258 for (int i=0;i<2;i++) 259 if ((!MatchFlag[i]) && (OtherAtom->type == ElementArray[i])) { 260 MatchFlag[i] = true; 261 break; // each bonding atom can match at most one element we are looking for 262 } 263 } 264 result = true; 265 for (int i=0;i<2;i++) // gather results 266 result = result && MatchFlag[i]; 267 if (result) { // check results 268 count++; 269 DoLog(1) && (Log() << Verbose(1) << first->name << "-" << second->name << "-" << third->name << " bond found at " << *Walker << "." << endl); 270 } 271 } 272 } 273 } 274 return count; 275 };
Note:
See TracChangeset
for help on using the changeset viewer.