source: src/analysis_bonds.cpp@ b32dbb

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since b32dbb was bfd839, checked in by Frederik Heber <heber@…>, 15 years ago

Changed and added counting of hydrogen bonds to menu.

  • CountHydrogenBonds() now receives two interface elements. This was needed to overcome double-couting in CSH-experiments.
  • Unit test CountHBonds has been changed accordingly, but not test for the new feature has been added.
  • case 'H' in builder.cpp is now used for counting hydrogen bonds not for giving help screen
  • option has been added to MeasureAtoms() menu

Signed-off-by: Frederik Heber <heber@…>

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