1 /* Copyright (c) 2015 Gerald Knizia
2 *
3 * This file is part of the IboView program (see: http://www.iboview.org)
4 *
5 * IboView is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * IboView is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with bfint (LICENSE). If not, see http://www.gnu.org/licenses/
16 *
17 * Please see IboView documentation in README.txt for:
18 * -- A list of included external software and their licenses. The included
19 * external software's copyright is not touched by this agreement.
20 * -- Notes on re-distribution and contributions to/further development of
21 * the IboView software
22 */
23
24 #include "Iv.h"
25 #include <boost/format.hpp>
26 #include <iostream>
27 #include <math.h> // for sqrt.
28
29 #include "CtIo.h"
30 #include "CtMatrix.h"
31 #include "CxColor.h"
32
33 #include "IvDocument.h"
34 #include "IvScript.h"
35 #include "IvIrc.h"
36 // #include "IvOrbitalFile.h"
37
38 #include "CtAtomData.h"
39 #include "CtBasisLibrary.h"
40 #include "CxPodArray.h"
41 using namespace ct;
42
43 static const bool g_UseOccupationNumbers = false;
44
45 // make the difference in charge for IBOs along all frames.
MakeIboChangeCurve(TArray<float> & CurveData,uint iRow,FDocument * document)46 bool MakeIboChangeCurve(TArray<float> &CurveData, uint iRow, FDocument *document)
47 {
48 // WARNING: due to the way the document loading is done, this function is
49 // invoking at least O(N^2) cost (any maybe O(N^3)) when loading lots of frames,
50 // since it is called after every frame.
51 // Should be fixed... might get expensive?
52 FOrbital
53 *pRefOrb = dynamic_cast<FOrbital*>(document->GetRowCol(iRow, 0, false));
54 if (!pRefOrb)
55 // data set is not there or is not an orbital.
56 return false;
57 uint
58 nCols = document->GetNumFrames();
59 CurveData.clear();
60 CurveData.resize(nCols, 0.);
61 TArray<double>
62 RefChg = pRefOrb->MakeIaoCharges(g_UseOccupationNumbers);
63 // PrintChargeArray(std::cout, "ref charge", RefChg);
64 if (RefChg.empty())
65 return false; // IAO coeffs not made.
66 for (uint iCol = 0; iCol < nCols; ++ iCol) {
67 FOrbital
68 *pCurOrb = dynamic_cast<FOrbital*>(document->GetRowCol(iRow, iCol, false));
69 TArray<double>
70 CurChg = pCurOrb->MakeIaoCharges(g_UseOccupationNumbers);
71 // PrintChargeArray(std::cout, "cur charge", RefChg);
72 double f = 0;
73 assert(RefChg.size() == CurChg.size());
74 for (uint iAt = 0; iAt < RefChg.size(); ++ iAt)
75 f += sqr(RefChg[iAt] - CurChg[iAt]);
76 CurveData[iCol] = 1.0 * std::sqrt(f);
77 }
78 return true;
79 }
80
81 // root mean square deviation between two vectors (i.e., 2-norm of their difference)
Rmsd(double const * pA,double const * pB,size_t N)82 static double Rmsd(double const *pA, double const *pB, size_t N) {
83 double RmsdSq = 0;
84 for (size_t i = 0; i < N; ++ i)
85 RmsdSq += sqr(pA[i] - pB[i]);
86 return sqrt(RmsdSq);
87 }
88
89 // estimate the IRC arc positions for all loaded frames. Optionally do not do the mass weighting.
MakeIrcArcLengths(TArray<float> & ArcLengths,FDocument * document,uint Flags)90 void MakeIrcArcLengths(TArray<float> &ArcLengths, FDocument *document, uint Flags)
91 {
92 ArcLengths.clear();
93
94 // get a pointer to the initial geometry.
95 IFrame
96 *pFrame0 = document->GetFrame(0);
97 if (pFrame0 == 0)
98 return;
99 FAtomSet
100 *pAtoms0 = pFrame0->pGetAtoms();
101 if (pAtoms0 == 0)
102 return;
103
104 // make a list of the atomic masses
105 FAtomicMassType
106 MassType = ATMASS_StandardAtomicWeight;
107 if (Flags & ARCLENGTH_UseIsotopeMasses)
108 MassType = ATMASS_MostCommonIsotope;
109 uint
110 nAt = pAtoms0->size();
111 TArray<double>
112 Masses(nAt, 1.);
113 if ((Flags & ARCLENGTH_NoMassWeighting) == 0) {
114 for (uint iAt = 0; iAt < nAt; ++ iAt)
115 Masses[iAt] = GetAtomicMass((*pAtoms0)[iAt].AtomicNumber, MassType);
116 }
117
118 // assemble a matrix of the mass weighted coordinates for all frames.
119 // The numbers and types of atoms in all frames must be equal.
120 uint
121 nFrames = (unsigned)document->GetNumFrames();
122 FMemoryStack2
123 Mem(3 * nAt * nFrames * sizeof(double) + 2000000);
124 FStackMatrix
125 Mwcs(3 * nAt, nFrames, &Mem);
126 for (uint iFrame = 0; iFrame < nFrames; ++ iFrame) {
127 IFrame
128 *pFrame = document->GetFrame(iFrame);
129 assert(pFrame != 0);
130 FAtomSet
131 *pAtoms = pFrame->pGetAtoms();
132 if (pAtoms == 0 || pAtoms->size() != nAt)
133 IvNotify(NOTIFY_Error, "IrcArcLength: Number of atoms not consistent across different frames.");
134 for (uint iAt = 0; iAt < nAt; ++ iAt) {
135 if ((*pAtoms)[iAt].AtomicNumber != (*pAtoms0)[iAt].AtomicNumber)
136 IvNotify(NOTIFY_Error, "IrcArcLength: Type of atoms not consistent across different frames.");
137 for (uint ixyz = 0; ixyz < 3; ++ ixyz)
138 Mwcs(3*iAt + ixyz, iFrame) = sqrt(Masses[iAt]) * (*pAtoms)[iAt].vPos[ixyz];
139 }
140 }
141
142 // approximate the IRC arc lengths via the differences of the 2-norm of the mass-weighted
143 // coordinates (note: if we have gradients we could calculate them more accurately).
144 ArcLengths.resize(nFrames, 0.);
145 double ArcLength = 0.;
146 ArcLengths[0] = ArcLength;
147 for (uint iFrame = 1; iFrame < nFrames; ++ iFrame) {
148 double Delta = Rmsd(&Mwcs(0, iFrame), &Mwcs(0, iFrame-1), Mwcs.nRows);
149 ArcLength += Delta;
150 ArcLengths[iFrame] = float(ArcLength);
151 }
152
153 // that's it already.
154 }
155