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