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