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