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 <cstdlib>
21 #include "wf_writer.h"
22 using namespace std;
print_wavefunction(ostream & inputfile)23 void Slat_wf_writer::print_wavefunction(ostream & inputfile ) {
24 int nmo=0;
25 if(calctype == "RHF" || calctype == "ROHF") {
26 nmo=max(nup, ndown);
27 }
28 else if(calctype == "UHF") {
29 nmo=spin_dwn_start+ndown;
30 }
31 else if(calctype=="GVB") {
32 assert(detwt.size() > 0);
33 nmo=0;
34 for(unsigned int i=0; i< occ_up.size(); i++) {
35 for(unsigned int j=0; j < occ_up[i].size(); j++)
36 if(nmo < occ_up[i][j]) nmo=occ_up[i][j];
37 }
38
39 for(unsigned int i=0; i< occ_down.size(); i++) {
40 for(unsigned int j=0; j < occ_down[i].size(); j++)
41 if(nmo < occ_down[i][j]) nmo=occ_down[i][j];
42 }
43 }
44 else {
45 cout << "strange calctype " << calctype << endl;
46 exit(1);
47 }
48
49 int nelectrons=nup+ndown;
50 if(magnification < 0)
51 magnification=max(nelectrons/20.0, 1.0);
52
53 inputfile << " SLATER \n"
54 << " " << orbtype << " {\n"
55 << " " << mo_matrix_type << endl
56 << " MAGNIFY " << magnification << endl
57 << " NMO " << nmo << endl
58 << " ORBFILE " << orbname << endl;
59 if(basisname!="") {
60 inputfile << " INCLUDE " << basisname << endl;
61 }
62
63 if(write_centers) {
64 inputfile << " CENTERS { READ " << centername << " } " << endl;
65 }
66 else if(use_global_centers) {
67 inputfile << " CENTERS { USEGLOBAL } " << endl;
68 }
69 else {
70 inputfile << " CENTERS { USEATOMS } " << endl;
71 }
72 inputfile << " }\n\n";
73
74 if(detwt.size()==0) {
75 //see below about the two-determinant UHF wave function
76 //if(calctype=="UHF" && ndown==nup)
77 // inputfile << " DETWT { .707 .707 } \n";
78 //else
79 inputfile << " DETWT { 1.0 } \n";
80
81
82 inputfile << " STATES {\n ";
83
84
85 inputfile << " #Spin up orbitals\n ";
86 for(int i=1; i< nup+1; i++) {
87 inputfile << i << " ";
88 if(i%20==0) inputfile << "\n ";
89 }
90 int downstart=1;
91 if(calctype=="UHF") {
92 downstart=spin_dwn_start+1;
93 }
94 inputfile << "\n ";
95
96 inputfile << " #Spin down orbitals \n ";
97 for(int i=downstart; i < downstart+ndown; i++) {
98 inputfile << i << " ";
99 if((i-downstart+1)%20==0) inputfile << "\n ";
100 }
101
102 //This doesn't necessarily help, and can sometimes hurt, so
103 //we don't put it in.
104 //if(calctype=="UHF" && ndown==nup) {
105 // inputfile <<"\n #for UHF, we put in two determinants \n";
106 // inputfile << " #Spin up orbitals\n";
107 // for(int i=downstart; i < downstart+ndown; i++) {
108 // inputfile << i << " ";
109 // if((i-downstart+1)%20==0) inputfile << "\n ";
110 // }
111 // inputfile << endl;
112 // inputfile << " #Spin down orbitals\n ";
113 // for(int i=1; i< nup+1; i++) {
114 // inputfile << i << " ";
115 // if(i%20==0) inputfile << "\n ";
116 // }
117 //}
118
119 inputfile << " }\n";
120 }
121 else {
122 int ndet=detwt.size();
123 assert(occ_up.size()==detwt.size());
124 assert(occ_down.size()==detwt.size());
125 inputfile << "DETWT { ";
126 for(int d=0; d< ndet; d++) {
127 inputfile << detwt[d] << " ";
128 }
129 inputfile << " } " << endl;
130
131 inputfile << " STATES { \n";
132 for(int d=0; d< ndet; d++) {
133 assert(occ_up[d].size()==nup);
134 inputfile << "# Up spin \n";
135 for(int i=0; i < nup; i++) {
136 inputfile << occ_up[d][i] << " ";
137 }
138 inputfile << endl;
139
140 assert(occ_down[d].size()==ndown);
141 inputfile << "# Down spin\n";
142 for(int i=0; i< ndown; i++) {
143 inputfile << occ_down[d][i] << " ";
144 }
145 inputfile << endl;
146 }
147 inputfile << "} " << endl;
148 }
149
150 }
151
152
153 //######################################################################
154
set_atoms(vector<Atom> & atoms)155 void Jastrow_wf_writer::set_atoms(vector <Atom> & atoms) {
156 int natoms=atoms.size();
157 for(int at=0; at < natoms; at++) {
158 bool unique=true;
159 int nunique=atomnames.size();
160 for(int i=0; i< nunique; i++) {
161 if(atoms[at].name==atomnames[i]) {
162 unique=false;
163 break;
164 }
165 }
166 if(unique) {
167 atomnames.push_back(atoms[at].name);
168 }
169 }
170 }
171
172
add_basis(Basis_writer & b)173 void Jastrow_wf_writer::add_basis(Basis_writer & b) {
174 basis.push_back(&b);
175 }
176
print_wavefunction(ostream & os)177 void Jastrow_wf_writer::print_wavefunction(ostream & os) {
178 os << "JASTROW \n\n";
179 for(vector <Basis_writer *>:: iterator i=basis.begin();
180 i!=basis.end(); i++) {
181 (*i)->print_basis(os);
182 }
183
184 for(vector <string>::iterator i=atomnames.begin();
185 i!=atomnames.end(); i++) {
186 os << "EICORRELATION { " << *i << " 0.0 0.0 0.0 }\n\n";
187 os << "EEICORRELATION { " << *i << endl
188 << " 1 1 0 0.\n"
189 << " 1 0 1 0.\n"
190 << " 1 1 1 0.\n"
191 << " 2 2 0 0.\n"
192 << " 2 0 1 0.\n"
193 << " 2 0 2 0.\n"
194 << " 2 2 2 0.\n"
195 << " 3 3 0 0.\n"
196 << " 3 0 2 0.\n"
197 << " 3 3 2 0.\n"
198 << " 1 2 2 0.\n"
199 << " 2 3 2 0.\n"
200 << "}\n\n";
201 }
202
203 }
204
205
206 //#######################################################################
set_atoms(vector<Atom> & atoms)207 void Jastrow2_wf_writer::set_atoms(vector <Atom> & atoms) {
208 int natoms=atoms.size();
209 for(int at=0; at < natoms; at++) {
210 bool unique=true;
211 int nunique=atomnames.size();
212 for(int i=0; i< nunique; i++) {
213 if(atoms[at].name==atomnames[i]) {
214 unique=false;
215 break;
216 }
217 }
218 if(unique) {
219 atomnames.push_back(atoms[at].name);
220 }
221 }
222 }
223
224
add_ee_basis(Basis_writer & b,int group)225 void Jastrow2_wf_writer::add_ee_basis(Basis_writer & b, int group) {
226 assert(group < ngroups);
227 ee_basis[group].push_back(&b);
228 }
229
add_ei_basis(Basis_writer & b,int group)230 void Jastrow2_wf_writer::add_ei_basis(Basis_writer & b, int group) {
231 assert(group < ngroups);
232 ei_basis[group].push_back(&b);
233 }
234
235
236
print_wavefunction(ostream & os)237 void Jastrow2_wf_writer::print_wavefunction(ostream & os) {
238 os << "JASTROW2 \n\n";
239 assert(ee_basis.size()==ei_basis.size());
240 assert(ee_basis.size()==ngroups);
241 //cout << "ngroups " << ngroups << endl;
242 string indent= " ";
243
244 for(int g=0; g< ngroups; g++) {
245 os << "GROUP { \n";
246 os << " OPTIMIZEBASIS " << endl;
247 for(vector <Basis_writer *>:: iterator i=ei_basis[g].begin();
248 i!=ei_basis[g].end(); i++) {
249 os << "EIBASIS { " << endl;
250 (*i)->print_basis(os);
251 os << "}" << endl;
252 }
253
254 if(ei_basis[g].size() > 0) {
255 os << "ONEBODY { " << endl;
256
257 for(vector <string>::iterator i=atomnames.begin();
258 i!=atomnames.end(); i++) {
259 int nfunc=0;
260 for(unsigned int bas=0;
261 bas< ei_basis[g].size(); bas++) {
262 if(*i == ei_basis[g][bas]->label) {
263 nfunc+=ei_basis[g][bas]->nfunc();
264 }
265 }
266
267 if(nfunc > 0) {
268 os << " COEFFICIENTS { " << *i << " ";
269 for(int i=0; i< nfunc; i++)
270 os << " 0.0 ";
271 os << " } " << endl;
272 }
273 }
274 os << " } " << endl;
275 }
276
277 for(vector <Basis_writer*> :: iterator i=ee_basis[g].begin();
278 i!=ee_basis[g].end(); i++) {
279 os << "EEBASIS { " << endl;
280 (*i)->print_basis(os);
281 os << "}" << endl;
282 }
283
284 if(spindiff[g]) {
285 assert(ee_basis[g].size()==2);
286
287 os << "TWOBODY_SPIN { \n";
288 os << " FREEZE \n";
289 os << " LIKE_COEFFICIENTS { 0.25 0.0 } \n";
290 os << " UNLIKE_COEFFICIENTS { 0.0 0.5 } \n";
291 os << " }\n";
292 }
293 else {
294 int nfunc=0;
295 for(vector <Basis_writer *>:: iterator i=ee_basis[g].begin();
296 i!=ee_basis[g].end(); i++) {
297 nfunc+=(*i)->nfunc();
298 }
299 if(nfunc > 0) {
300 os << "TWOBODY { " << endl;
301 os << " COEFFICIENTS { " ;
302 for(int j=0; j< nfunc; j++)
303 os << " 0.0 ";
304 os << " }" << endl;
305 os << " } " << endl;
306 }
307 }
308 os << "} " << endl;
309
310 }
311
312 }
313
314
315 //----------------------------------------------------------------------
316
print_std_jastrow2(Jastrow2_wf_writer & jast2writer,ostream & os,double basis_cutoff)317 void print_std_jastrow2(Jastrow2_wf_writer & jast2writer, ostream & os,
318 double basis_cutoff) {
319 jast2writer.set_groups(2);
320 jast2writer.set_spin_diff(0);
321
322 vector <Cutoff_cusp_basis> eecusp; //Group 0 is the cusp
323 eecusp.resize(2);
324 eecusp[0].cutoff=eecusp[1].cutoff=basis_cutoff;
325 jast2writer.add_ee_basis(eecusp[0], 0);
326 jast2writer.add_ee_basis(eecusp[1], 0);
327
328 //ugly hack.. The jast writer saves pointers to the basis functions,
329 //but when we modify a vector, those pointers can change, so
330 //we have to add all the basis functions, then assign them.
331 //This is pretty unstable..there has to be a better way of doing
332 //it without losing polymorphism.
333
334 vector <Poly_pade_basis> eibas; //group 2 is the regular basis
335 for(unsigned int i=0; i< jast2writer.atomnames.size(); i++) {
336 Poly_pade_basis tbas;
337 tbas.label=jast2writer.atomnames[i];
338 tbas.cutoff=basis_cutoff;
339 tbas.beta0=.2;
340 tbas.nfunc_=3;
341 eibas.push_back(tbas);
342
343 }
344 for(unsigned int i=0; i< eibas.size(); i++) {
345 jast2writer.add_ei_basis(eibas[i], 1);
346 }
347 Poly_pade_basis eebas;
348 eebas.cutoff=basis_cutoff;
349 eebas.beta0=.5;
350 eebas.nfunc_=3;
351 jast2writer.add_ee_basis(eebas,1);
352 jast2writer.print_wavefunction(os);
353 }
354
355
356 //----------------------------------------------------------------------
357 //this is done in a much more straightforward way than the Jastrow2_writer
358 //I don't think that we need all the polymorphism, etc, that we had before
print_3b_jastrow2(ostream & os,std::vector<std::string> & unique_atoms,double cutoff)359 void print_3b_jastrow2(ostream & os, std::vector<std::string> & unique_atoms, double cutoff) {
360 os << "JASTROW2" << endl;
361 os << "GROUP { \n"
362 << " OPTIMIZEBASIS\n"
363 << " EEBASIS { EE CUTOFF_CUSP GAMMA 24.0 CUSP 1.0 CUTOFF " << cutoff << " } \n"
364 << " EEBASIS { EE CUTOFF_CUSP GAMMA 24.0 CUSP 1.0 CUTOFF " << cutoff << " } \n"
365 << " TWOBODY_SPIN { FREEZE \n"
366 << " LIKE_COEFFICIENTS { 0.25 0.0 } \n"
367 << " UNLIKE_COEFFICIENTS { 0.0 0.5 } \n }\n"
368 << "}\n";
369 os << "GROUP { \n OPTIMIZEBASIS\n"
370 <<" EEBASIS { EE POLYPADE BETA0 0.5 NFUNC 3 RCUT " << cutoff << " }\n";
371 for(vector<string>::iterator i=unique_atoms.begin(); i!= unique_atoms.end();
372 i++) {
373 os << " EIBASIS { " << *i << " POLYPADE BETA0 0.2 NFUNC 4 RCUT " << cutoff << " } \n";
374 }
375 os << " ONEBODY { \n";
376 for(vector<string>::iterator i=unique_atoms.begin(); i!= unique_atoms.end();
377 i++) {
378 os << " COEFFICIENTS { " << *i << " 0 0 0 0 } \n";
379 }
380 os << " }\n";
381 os << " TWOBODY { \n";
382 os << " COEFFICIENTS { 0 0 0 } \n";
383 os << " }\n";
384 os << " THREEBODY { \n";
385 for(vector<string>::iterator i=unique_atoms.begin(); i!= unique_atoms.end();
386 i++) {
387 os << " COEFFICIENTS { " << *i << " 0 0 0 0 0 0 0 0 0 0 0 0 } \n";
388 }
389 os << " }\n";
390 os << "}\n";
391
392 }
393 //----------------------------------------------------------------------
394