source: src/analysis_bonds.cpp@ 3839e5

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 3839e5 was 36166d, checked in by Tillmann Crueger <crueger@…>, 14 years ago

Removed left over parts from old memory-tracker

  • Property mode set to 100644
File size: 11.3 KB
RevLine 
[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 */
25void 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]54void 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 */
89double 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 */
107bool 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]126int 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 */
204int 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 */
235int 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};
Note: See TracBrowser for help on using the repository browser.