1 /*
2
3 Copyright (C) 2007 Lucas K. Wagner
4
5 This program 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; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License along
16 with this program; if not, write to the Free Software Foundation, Inc.,
17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18
19 */
20 #include "Array.h"
21
22
23 /*
24 * usage: compare_mo <nfunc> <nmo> <file1> <file2>
25 * nfunc: total number of basis functions
26 * nmo: total number of orbitals to compare between the files
27 * file1 and file2 should be QWalk .orb files
28 * Will output the dot products between the orbitals and check
29 * whether any of the orbitals have changed significantly.
30 *
31 * This initially was meant for Jeff Grossman's continuous DMC
32 * algorithm, in which it's necessary to check for changes in orbital
33 * occupation, but can be used for other analysis purposes as well.
34 * */
35
read_mo(string & orbin,int nmo,int nfunc,Array2<double> & moCoeff)36 void read_mo(string & orbin, int nmo, int nfunc, Array2 <double> & moCoeff)
37 {
38 ifstream in(orbin.c_str());
39 if(!in) error("Couldn't open ", orbin);
40 string dummy;
41 in >> dummy;
42 while(dummy != "COEFFICIENTS")
43 in >> dummy;
44
45 moCoeff.Resize(nmo, nfunc);
46 for(int m=0; m< nmo; m++) {
47 for(int f=0; f< nfunc; f++) {
48 in >> moCoeff(m,f);
49 }
50 }
51 }
52
53 #include <iomanip>
54 /*!
55 Compare the sets of molecular orbitals, and
56 return true if they're the same, and false if not
57 Currently it doesn't use the lcao_overlap matrix, but
58 should be modified to do that; it's approximate otherwise.
59 */
compare_mo(Array2<double> & oldMOCoeff,Array2<double> & newMOCoeff,Array2<double> & lcao_overlap,Array1<int> compare_list)60 bool compare_mo(Array2 <double> & oldMOCoeff,
61 Array2 <double> & newMOCoeff,
62 Array2 <double> & lcao_overlap,
63 Array1 <int> compare_list) {
64
65 assert(oldMOCoeff.GetDim(1)==newMOCoeff.GetDim(1));
66
67 int nfunctions=oldMOCoeff.GetDim(1);
68
69 int ncompare=compare_list.GetDim(0);
70 if(newMOCoeff.GetDim(0) < ncompare)
71 error("the punch file doesn't have enough molecular orbitals to compare.");
72 if(oldMOCoeff.GetDim(0) < ncompare)
73 error("the old punch file doesn't have enough MO's to compare.");
74
75 vector <int> unresolved_mos;
76
77 //First check to see if the mo's are in the same place
78 //(most should be)
79 cout.precision(15);
80 for(int i=0; i< ncompare; i++) {
81 int mo=compare_list(i);
82 double dot=0, mag_old=0, mag_new=0;
83 //cout << "-----------------------mo " << mo << endl;
84 for(int f=0; f< nfunctions; f++) {
85 //cout << setw(10) << "function " << setw(4) << f
86 // << " new " << setw(15) << newMOCoeff(mo,f)
87 // << " old " << setw(15) << oldMOCoeff(mo,f)
88 // << " diff "
89 // << setw(15) << newMOCoeff(mo,f)-oldMOCoeff(mo,f)
90 // << endl;
91 dot+=newMOCoeff(mo,f)*oldMOCoeff(mo,f);
92 mag_old+=oldMOCoeff(mo,f)*oldMOCoeff(mo,f);
93 mag_new+=newMOCoeff(mo,f)*newMOCoeff(mo,f);
94 }
95 dot /= sqrt(mag_old*mag_new);
96 dot =fabs(dot);
97 cout << "mo " << mo << " dot " << dot << endl;
98 if(fabs(dot-1) > .01) {
99 unresolved_mos.push_back(mo);
100 }
101 }
102
103 int nunresolved=unresolved_mos.size();
104 for(int i=0; i< nunresolved; i++) {
105 cout << "not matched: " << unresolved_mos[i] << endl;
106 }
107 bool are_same=true;
108 //See if any just swapped..
109 for(int i=0; i< nunresolved; i++) {
110 int mo1=unresolved_mos[i];
111 bool resolved_swapping=false;
112
113 for(int j=0; j< nunresolved; j++) {
114 int mo2=unresolved_mos[j];
115
116 double dot=0, mag_old=0, mag_new=0;
117 for(int f=0; f< nfunctions; f++) {
118 dot+=newMOCoeff(mo1,f)*oldMOCoeff(mo2,f);
119 mag_old+=oldMOCoeff(mo2,f)*oldMOCoeff(mo2,f);
120 mag_new+=newMOCoeff(mo1,f)*newMOCoeff(mo1,f);
121 }
122 dot /= sqrt(mag_old*mag_new);
123 dot=fabs(dot);
124 if(fabs(dot-1) < .01) {
125 cout << "switched orbital: mo " << mo2 << " went to " << mo1
126 << " dot product " << dot
127 << endl;
128 resolved_swapping=true;
129 }
130 }
131
132 if(!resolved_swapping) {
133 cout << "Unresolvable change in mo " << mo1 << endl;
134 are_same=false;
135 }
136
137 }
138
139 return are_same;
140 }
141
142 //----------------------------------------------------------------------
143
main(int argc,char * argv[])144 int main(int argc, char * argv[]) {
145
146 if(argc <5) {
147 cout << "usage: compare_mo <nfunc> <nmo> <file1> <file2>" << endl;
148 exit(1);
149 }
150 cout.precision(15);
151
152 string file1=argv[3];
153 string file2=argv[4];
154 int nmo=atoi(argv[2]);
155 int nfunc=atoi(argv[1]);
156 Array2 <doublevar> mo1(nmo, nfunc), mo2(nmo, nfunc);
157 Array2 <doublevar> lcao_overlap(nfunc, nfunc);
158 lcao_overlap=0.0;
159 for(int i=0; i< nfunc; i++)
160 lcao_overlap(i,i)=1.0;
161
162 Array1<int> compare_list(nmo);
163 for(int m=0; m< nmo; m++) compare_list(m)=m;
164
165 read_mo(file1, nmo, nfunc, mo1);
166 read_mo(file2, nmo, nfunc, mo2);
167 compare_mo(mo1, mo2, lcao_overlap, compare_list);
168 }
169