1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * HF/DFT from Hel
7 *
8 * Copyright © 2015 The Regents of the University of California
9 * All Rights Reserved
10 *
11 * Written by Susi Lehtola, Lawrence Berkeley National Laboratory
12 *
13 * This program is free software; you can redistribute it and/or
14 * modify it under the terms of the GNU General Public License
15 * as published by the Free Software Foundation; either version 2
16 * of the License, or (at your option) any later version.
17 */
18
19 #include "../settings.h"
20 #include "../checkpoint.h"
21
getlegend(size_t i,size_t n)22 std::string getlegend(size_t i, size_t n) {
23 if(n==1)
24 return "orbital";
25 else if(n==2) {
26 if(i)
27 return "beta orbital";
28 else
29 return "alpha orbital";
30 } else
31 throw std::logic_error("Should have not ended up here.\n");
32 }
33
34 Settings settings;
35
main_guarded(int argc,char ** argv)36 int main_guarded(int argc, char **argv) {
37 if(argc!=3) {
38 printf("Usage: %s chkpt1 chkpt2",argv[0]);
39 return 1;
40 }
41
42 Checkpoint lhchk(argv[1],false);
43 Checkpoint rhchk(argv[2],false);
44
45 BasisSet lhbas;
46 lhchk.read(lhbas);
47
48 BasisSet rhbas;
49 rhchk.read(rhbas);
50
51 // Get basis set overlaps
52 arma::mat S12=lhbas.overlap(rhbas);
53
54 // Restricted calcs?
55 int lhrestr, rhrestr;
56 lhchk.read("Restricted",lhrestr);
57 rhchk.read("Restricted",rhrestr);
58
59 // Orbital coefficient matrices
60 std::vector<arma::cx_mat> lmat, rmat;
61 if(lhrestr) {
62 lmat.resize(1);
63 if(lhchk.exist("CW.re"))
64 lhchk.cread("CW",lmat[0]);
65 else {
66 arma::mat C;
67 lhchk.read("C",C);
68 lmat[0]=C*COMPLEX1;
69 }
70
71 int Nel;
72 lhchk.read("Nel-a",Nel);
73 if(lmat[0].n_cols> (size_t) Nel)
74 lmat[0].shed_cols(Nel,lmat[0].n_cols-1);
75
76 } else {
77 lmat.resize(2);
78 if(lhchk.exist("CWa.re")) {
79 lhchk.cread("CWa",lmat[0]);
80 lhchk.cread("CWb",lmat[1]);
81 } else {
82 arma::mat Ca, Cb;
83 lhchk.read("Ca",Ca);
84 lhchk.read("Cb",Cb);
85 lmat[0]=Ca*COMPLEX1;
86 lmat[1]=Cb*COMPLEX1;
87 }
88
89 int Nel;
90 lhchk.read("Nel-a",Nel);
91 if(lmat[0].n_cols> (size_t) Nel)
92 lmat[0].shed_cols(Nel,lmat[0].n_cols-1);
93 lhchk.read("Nel-b",Nel);
94 if(lmat[1].n_cols> (size_t) Nel)
95 lmat[1].shed_cols(Nel,lmat[1].n_cols-1);
96 }
97
98 if(rhrestr) {
99 rmat.resize(1);
100 if(rhchk.exist("CW.re"))
101 rhchk.cread("CW",rmat[0]);
102 else {
103 arma::mat C;
104 rhchk.read("C",C);
105 rmat[0]=C*COMPLEX1;
106 }
107
108 int Nel;
109 rhchk.read("Nel-a",Nel);
110 if(rmat[0].n_cols> (size_t) Nel)
111 rmat[0].shed_cols(Nel,rmat[0].n_cols-1);
112 } else {
113 rmat.resize(2);
114 if(rhchk.exist("CWa.re")) {
115 rhchk.cread("CWa",rmat[0]);
116 rhchk.cread("CWb",rmat[1]);
117 } else {
118 arma::mat Ca, Cb;
119 rhchk.read("Ca",Ca);
120 rhchk.read("Cb",Cb);
121 rmat[0]=Ca*COMPLEX1;
122 rmat[1]=Cb*COMPLEX1;
123 }
124
125 int Nel;
126 rhchk.read("Nel-a",Nel);
127 if(rmat[0].n_cols> (size_t) Nel)
128 rmat[0].shed_cols(Nel,rmat[0].n_cols-1);
129 rhchk.read("Nel-b",Nel);
130 if(rmat[1].n_cols> (size_t) Nel)
131 rmat[1].shed_cols(Nel,rmat[1].n_cols-1);
132 }
133
134 for(size_t il=0;il<lmat.size();il++)
135 for(size_t ir=0;ir<rmat.size();ir++) {
136
137 if(lmat.size()==2 && rmat.size()==2 && ((il && !ir) || (ir && !il)))
138 continue;
139
140 // Orbital projections
141 arma::mat proj;
142 {
143 // With or without complex conjugate?
144 arma::cx_mat Smot(arma::trans(lmat[il])*S12*rmat[ir]);
145 arma::mat prot(arma::real(Smot%arma::conj(Smot)));
146
147 arma::cx_mat Smost(arma::strans(lmat[il])*S12*rmat[ir]);
148 arma::mat prost(arma::real(Smost%arma::conj(Smost)));
149
150 // Sums
151 double sumt(arma::sum(arma::sum(prot)));
152 double sumst(arma::sum(arma::sum(prost)));
153
154 if(sumt>sumst) {
155 printf("Determinants do not differ by complex conjugation\n");
156 proj=prot;
157 } else {
158 printf("Determinants differ by complex conjugation\n");
159 proj=prost;
160 }
161 }
162
163 printf("%s - %s projection\n",getlegend(il,lmat.size()).c_str(),getlegend(ir,rmat.size()).c_str());
164 //proj.print();
165
166 if(lmat[il].n_cols == rmat[ir].n_cols) {
167 printf("Trace of density projection is % 10.6f, deviation from norm is %e\n",arma::sum(arma::sum(proj)),rmat[ir].n_cols-arma::sum(arma::sum(proj)));
168
169 // Determine ordering
170 arma::uvec lorder(arma::linspace<arma::uvec>(0,lmat[il].n_cols-1,lmat[il].n_cols));
171 arma::uvec rorder(arma::linspace<arma::uvec>(0,rmat[ir].n_cols-1,rmat[ir].n_cols));
172 for(size_t i=0;i<lorder.n_elem;i++) {
173 // Sorted matrix
174 arma::mat sproj(proj(lorder,rorder));
175 // Take submatrix
176 sproj=sproj.submat(i,i,sproj.n_rows-1,sproj.n_cols-1);
177
178 // Find maximum element
179 arma::uword lind, rind;
180 sproj.max(lind,rind);
181 lind+=i;
182 rind+=i;
183
184 // Swap indices
185 std::swap(lorder(i),lorder(lind));
186 std::swap(rorder(i),rorder(rind));
187 }
188
189 lorder.t().print("Left ordering");
190 rorder.t().print("Right ordering");
191
192 // Sort projection matrix
193 arma::mat sproj(proj(lorder,rorder));
194 sproj.print("Projection matrix");
195
196 arma::vec osort(arma::sort(arma::diagvec(sproj),"descend"));
197 printf("Sum of diagonal orbital projections is % 10.6f, deviation from norm is %e\n",arma::sum(osort),rmat[ir].n_cols-arma::sum(osort));
198
199 osort-=arma::ones<arma::vec>(osort.n_elem);
200 osort.print("Minimal orbital differences");
201 } else {
202 proj.print();
203 }
204
205 printf("\n");
206 }
207
208 return 0;
209 }
210
main(int argc,char ** argv)211 int main(int argc, char **argv) {
212 try {
213 return main_guarded(argc, argv);
214 } catch (const std::exception &e) {
215 std::cerr << "error: " << e.what() << std::endl;
216 return 1;
217 }
218 }
219