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 "../checkpoint.h"
20 #include "../density_fitting.h"
21 #include "../erichol.h"
22 #include "../localization.h"
23 #include "../timer.h"
24 #include "../mathf.h"
25 #include "../linalg.h"
26 #include "../eriworker.h"
27 #include "../settings.h"
28 #include "../stringutil.h"
29
30 #include <cstdio>
31
32 #ifdef _OPENMP
33 #include <omp.h>
34 #endif
35
36 #ifdef SVNRELEASE
37 #include "../version.h"
38 #endif
39
40 Settings settings;
41
main_guarded(int argc,char ** argv)42 int main_guarded(int argc, char **argv) {
43 #ifdef _OPENMP
44 printf("ERKALE - Checkpoints from Hel, OpenMP version, running on %i cores.\n",omp_get_max_threads());
45 #else
46 printf("ERKALE - Checkpoints from Hel, serial version.\n");
47 #endif
48 print_copyright();
49 print_license();
50 #ifdef SVNRELEASE
51 printf("At svn revision %s.\n\n",SVNREVISION);
52 #endif
53 print_hostname();
54
55 if(argc!=2) {
56 printf("Usage: $ %s (runfile)\n",argv[0]);
57 return 0;
58 }
59
60 // Initialize libint
61 init_libint_base();
62
63 Timer t;
64
65 // Parse settings
66 settings.add_string("LoadChk","Checkpoint file that contains basis set etc","erkale.chk");
67 settings.add_string("SaveChk","Checkpoint to which orbitals are saved","chkpt.chk");
68 settings.add_string("Model","Model from which orbitals are","");
69 settings.add_bool("Binary","Use binary I/O?",true);
70 settings.add_string("Deleted","Add in deleted orbitals from file","");
71 settings.add_bool("NO","Are these natural orbitals?",false);
72
73 settings.parse(argv[1]);
74 settings.print();
75
76 // Load checkpoint
77 std::string loadchk(settings.get_string("LoadChk"));
78 std::string savechk(settings.get_string("SaveChk"));
79 std::string model(settings.get_string("Model"));
80 std::string del(settings.get_string("Deleted"));
81 bool binary(settings.get_bool("Binary"));
82 bool NO(settings.get_bool("NO"));
83
84 if(savechk.size()) {
85 // Copy checkpoint data
86 std::ostringstream cmd;
87 cmd << "cp " << loadchk << " " << savechk;
88 if(system(cmd.str().c_str()))
89 throw std::runtime_error("Error copying checkpoint file.\n");
90 }
91
92 // Save checkpoint
93 Checkpoint chkpt(savechk,true,false);
94
95 // File type
96 arma::file_type atype = binary ? arma::arma_binary : arma::raw_ascii;
97
98 // Load orbitals
99 arma::mat Ca, Cb;
100 arma::mat Pa, Pb;
101
102 if(NO) {
103 // MOs are in
104 Ca.load("Ca.dat",atype);
105 Cb.load("Cb.dat",atype);
106
107 // Load density matrices
108 Pa.load("Pa_" + model + "_MO.dat",atype);
109 Pb.load("Pb_" + model + "_MO.dat",atype);
110 // Convert to AO basis
111 Pa=Ca*Pa*arma::trans(Ca);
112 Pb=Cb*Pb*arma::trans(Cb);
113
114 // and the NO coefficients are in
115 arma::mat Wa, Wb;
116 Wa.load("Ca_"+model+"_MO.dat",atype);
117 Wb.load("Cb_"+model+"_MO.dat",atype);
118
119 // Rotate MOs to NO basis
120 Ca*=Wa;
121 Cb*=Wb;
122
123 } else {
124 Ca.load("Ca_"+model+".dat",atype);
125 Cb.load("Cb_"+model+".dat",atype);
126
127 // Load density matrices
128 Pa.load("Pa_" + model + "_MO.dat",atype);
129 Pb.load("Pb_" + model + "_MO.dat",atype);
130 // Convert to AO basis
131 Pa=Ca*Pa*arma::trans(Ca);
132 Pb=Cb*Pb*arma::trans(Cb);
133 }
134
135 // Was earlier checkpoint restricted?
136 int orestr;
137 chkpt.read("Restricted",orestr);
138 if(orestr) {
139 chkpt.remove("C");
140 } else {
141 chkpt.remove("Ca");
142 chkpt.remove("Cb");
143 }
144
145 // Save orbitals
146 int restr;
147 if(arma::norm(Ca-Cb,2)==0.0) {
148 restr=1;
149 chkpt.remove("Ha");
150 chkpt.remove("Hb");
151 chkpt.write("C",Ca);
152 chkpt.write("P",Pa+Pb);
153 } else {
154 restr=0;
155 chkpt.remove("H");
156 chkpt.write("Ca",Ca);
157 chkpt.write("Cb",Cb);
158 chkpt.write("P",Pa+Pb);
159 chkpt.write("Pa",Pa);
160 chkpt.write("Pb",Pb);
161 }
162 chkpt.write("Restricted",restr);
163
164 if(orestr && !restr) {
165 arma::vec E;
166 chkpt.read("E",E);
167 chkpt.remove("E");
168 chkpt.write("Ea",E);
169 chkpt.write("Eb",E);
170 } else if(!orestr && restr) {
171 arma::vec E;
172 chkpt.read("Ea",E);
173 chkpt.remove("Ea");
174 chkpt.remove("Eb");
175 chkpt.write("E",E);
176 }
177
178 if(NO) {
179 if(restr) {
180 // Load occupation numbers
181 arma::vec E;
182 E.load("NOON_" + model + "_alpha.dat",atype);
183 chkpt.write("E",E);
184 } else {
185 // Load occupation numbers
186 arma::vec Ea, Eb;
187 Ea.load("NOON_" + model + "_alpha.dat",atype);
188 Eb.load("NOON_" + model + "_beta.dat",atype);
189 chkpt.write("Ea",Ea);
190 chkpt.write("Eb",Eb);
191 }
192 }
193
194 if(del.size()) {
195 // Load orbitals from file
196 arma::mat Cd;
197 Cd.load(del,atype);
198
199 arma::vec zerocc(Cd.n_cols);
200 zerocc.zeros();
201
202 if(restr) {
203 arma::mat C;
204 chkpt.read("C",C);
205
206 arma::vec E;
207 chkpt.read("E",E);
208
209 arma::mat Cnew(arma::join_rows(C,Cd));
210 chkpt.write("C",Cnew);
211
212 if(NO) {
213 arma::vec Enew(E.n_elem+zerocc.n_elem);
214 Enew.subvec(0,E.n_elem-1)=E;
215 Enew.subvec(E.n_elem,Enew.n_elem-1)=zerocc;
216 chkpt.write("E",Enew);
217 }
218 } else {
219 chkpt.read("Ca",Ca);
220 chkpt.read("Cb",Cb);
221
222 arma::vec Ea, Eb;
223 chkpt.read("Ea",Ea);
224 chkpt.read("Eb",Eb);
225
226 arma::mat Canew(arma::join_rows(Ca,Cd));
227 chkpt.write("Ca",Canew);
228 arma::mat Cbnew(arma::join_rows(Cb,Cd));
229 chkpt.write("Cb",Cbnew);
230
231 if(NO) {
232 arma::vec Eanew(Ea.n_elem+zerocc.n_elem);
233 Eanew.subvec(0,Ea.n_elem-1)=Ea;
234 Eanew.subvec(Ea.n_elem,Eanew.n_elem-1)=zerocc;
235 chkpt.write("Ea",Eanew);
236
237 arma::vec Ebnew(Eb.n_elem+zerocc.n_elem);
238 Ebnew.subvec(0,Eb.n_elem-1)=Eb;
239 Ebnew.subvec(Eb.n_elem,Ebnew.n_elem-1)=zerocc;
240 chkpt.write("Eb",Ebnew);
241 }
242 }
243
244 printf("Added in %i deleted virtuals from %s.\n",(int) Cd.n_cols,del.c_str());
245 }
246
247 return 0;
248 }
249
main(int argc,char ** argv)250 int main(int argc, char **argv) {
251 try {
252 return main_guarded(argc, argv);
253 } catch (const std::exception &e) {
254 std::cerr << "error: " << e.what() << std::endl;
255 return 1;
256 }
257 }
258