source: ThirdParty/mpqc_open/src/bin/mpqc/mpqcin.cc@ e9f307

Candidate_v1.7.0 stable
Last change on this file since e9f307 was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 17.4 KB
RevLine 
[0b990d]1
2#include <scconfig.h>
3#ifdef HAVE_SSTREAM
4# include <sstream>
5#else
6# include <strstream.h>
7#endif
8#include <stdlib.h>
9
10#include <util/misc/formio.h>
11
12using namespace std;
13using namespace sc;
14
15#undef yyFlexLexer
16#define yyFlexLexer MPQCInFlexLexer
17#include <FlexLexer.h>
18
19#include "mpqcin.h"
20#include "parse.h"
21
22int MPQCIn::checking_ = 0;
23
24MPQCIn::MPQCIn():
25 nirrep_(0),
26 mult_(1),
27 charge_(0),
28 basis_(0),
29 auxbasis_(0),
30 method_(0),
31 optimize_(0),
32 gradient_(0),
33 frequencies_(0),
34 opt_type_(T_INTERNAL),
35 redund_coor_(0),
36 restart_(0),
37 checkpoint_(1),
38 atom_charge_(0),
39 method_xc_(0),
40 method_grid_(0),
41 symmetry_(0),
42 memory_(0),
43 molecule_bohr_(0),
44 alpha_(0),
45 beta_(0),
46 docc_(0),
47 socc_(0),
48 frozen_docc_(0),
49 frozen_uocc_(0),
50 method_absmethod_(0),
51 method_ebc_(0),
52 method_gbc_(0)
53{
54 lexer_ = new MPQCInFlexLexer;
55}
56
57MPQCIn::~MPQCIn()
58{
59 delete lexer_;
60 if (basis_.val()) free(basis_.val());
61 if (auxbasis_.val()) free(auxbasis_.val());
62 if (method_.val()) free(method_.val());
63 if (method_xc_.val()) free(method_xc_.val());
64 if (method_grid_.val()) free(method_grid_.val());
65 if (symmetry_.val()) free(symmetry_.val());
66 if (memory_.val()) free(memory_.val());
67 if (alpha_.val()) free(alpha_.val());
68 if (beta_.val()) free(beta_.val());
69 if (docc_.val()) free(docc_.val());
70 if (socc_.val()) free(socc_.val());
71 if (frozen_docc_.val()) free(frozen_docc_.val());
72 if (frozen_uocc_.val()) free(frozen_uocc_.val());
73}
74
75void
76MPQCIn::error(const char* s)
77{
78 ExEnv::outn() << ExEnv::program_name()
79 << ": error: " << s
80 << endl;
81 abort();
82}
83
84void
85MPQCIn::error2(const char* s, const char *s2)
86{
87 ExEnv::outn() << ExEnv::program_name()
88 << ": error: " << s << "\"" << s2 << "\""
89 << endl;
90 abort();
91}
92
93void
94MPQCIn::yerror(const char* s)
95{
96 ExEnv::outn() << ExEnv::program_name()
97 << ": " << s
98 << " at line " << lexer_->lineno()+1
99 << endl;
100 abort();
101}
102
103void
104MPQCIn::yerror2(const char* s, const char *s2)
105{
106 ExEnv::outn() << ExEnv::program_name()
107 << ": " << s
108 << " \"" << s2 << "\" at line " << lexer_->lineno()+1
109 << endl;
110 abort();
111}
112
113int
114MPQCIn::ylex()
115{
116 return lexer_->yylex();
117}
118
119void
120MPQCIn::begin_molecule()
121{
122 if (mol_.nonnull()) {
123 ExEnv::outn() << ExEnv::program_name()
124 << ": error: second molecule given at line "
125 << lexer_->lineno()+1
126 << endl;
127 abort();
128 }
129 mol_ = new Molecule;
130}
131
132void
133MPQCIn::end_molecule()
134{
135 //double symtol = 1e-4;
136 //mol_->set_point_group(mol_->highest_point_group(symtol), symtol*10.0);
137}
138
139void
140MPQCIn::add_atom(char *sym, char *xs, char *ys, char *zs)
141{
142 int Z = mol_->atominfo()->string_to_Z(sym, 0);
143 if (Z == 0) yerror2("bad element", sym);
144 free(sym);
145
146 char *xse;
147 double x = strtod(xs,&xse);
148 if (xse == xs) yerror2("bad x coordinate", xs);
149 free(xs);
150
151 char *yse;
152 double y = strtod(ys,&yse);
153 if (yse == ys) yerror2("bad y coordinate", ys);
154 free(ys);
155
156 char *zse;
157 double z = strtod(zs,&zse);
158 if (zse == zs) yerror2("bad z coordinate", zs);
159 free(zs);
160
161 mol_->add_atom(Z, x, y, z, 0, 0, atom_charge_.set(), atom_charge_.val());
162 atom_charge_.reset(0);
163}
164
165void
166MPQCIn::set_charge(char *cs)
167{
168 char *cse;
169 int c = strtol(cs,&cse,10);
170 if (cse == cs) yerror2("bad charge", cs);
171 charge_ = c;
172 free(cs);
173}
174
175void
176MPQCIn::set_atom_charge(char *cs)
177{
178 char *cse;
179 int c = strtol(cs,&cse,10);
180 if (cse == cs) yerror2("bad atom charge", cs);
181 atom_charge_ = c;
182 free(cs);
183}
184
185void
186MPQCIn::set_method(char *m)
187{
188 method_ = m;
189}
190
191void
192MPQCIn::set_method_xc(char *m)
193{
194 method_xc_ = m;
195}
196
197void
198MPQCIn::set_method_grid(char *m)
199{
200 method_grid_ = m;
201}
202
203void
204MPQCIn::set_molecule_bohr(int i)
205{
206 molecule_bohr_ = i;
207}
208
209void
210MPQCIn::set_basis(char *b)
211{
212 basis_ = b;
213}
214
215void
216MPQCIn::set_auxbasis(char *b)
217{
218 auxbasis_ = b;
219}
220
221void
222MPQCIn::set_symmetry(char *s)
223{
224 symmetry_ = s;
225 if (strcmp(s,"auto")) {
226 Ref<PointGroup> p = new PointGroup(s);
227 nirrep_ = p->char_table().nirrep();
228 }
229}
230
231void
232MPQCIn::set_memory(char *s)
233{
234 memory_ = s;
235}
236
237void
238MPQCIn::set_multiplicity(char *ms)
239{
240 char *mse;
241 int m = strtol(ms,&mse,10);
242 if (mse == ms || m <= 0) yerror2("bad multiplicity", ms);
243 mult_ = m;
244 free(ms);
245}
246
247void
248MPQCIn::set_optimize(int i)
249{
250 optimize_ = i;
251}
252
253void
254MPQCIn::set_gradient(int i)
255{
256 gradient_ = i;
257}
258
259void
260MPQCIn::set_frequencies(int i)
261{
262 frequencies_ = i;
263}
264
265void
266MPQCIn::set_restart(int i)
267{
268 restart_ = i;
269}
270
271void
272MPQCIn::set_checkpoint(int i)
273{
274 checkpoint_ = i;
275}
276
277void
278MPQCIn::set_redund_coor(int i)
279{
280 redund_coor_ = i;
281}
282
283void
284MPQCIn::set_opt_type(int i)
285{
286 opt_type_ = i;
287}
288
289void
290MPQCIn::set_docc(std::vector<int> *a)
291{
292 docc_ = a;
293}
294
295void
296MPQCIn::set_socc(std::vector<int> *a)
297{
298 socc_ = a;
299}
300
301void
302MPQCIn::set_alpha(std::vector<int> *a)
303{
304 alpha_ = a;
305}
306
307void
308MPQCIn::set_beta(std::vector<int> *a)
309{
310 beta_ = a;
311}
312
313void
314MPQCIn::set_frozen_docc(std::vector<int> *a)
315{
316 frozen_docc_ = a;
317}
318
319void
320MPQCIn::set_frozen_uocc(std::vector<int> *a)
321{
322 frozen_uocc_ = a;
323}
324
325void
326MPQCIn::set_method_absmethod(const char *m)
327{
328 method_absmethod_ = m;
329}
330
331void
332MPQCIn::set_method_ebc(const char *m)
333{
334 method_ebc_ = m;
335}
336
337void
338MPQCIn::set_method_gbc(const char *m)
339{
340 method_gbc_ = m;
341}
342
343std::vector<int> *
344MPQCIn::make_nnivec(std::vector<int> *a, char *ms)
345{
346 if (ms == 0) return new std::vector<int>;
347
348 char *mse;
349 int m = strtol(ms,&mse,10);
350 if (mse == ms || m < 0) yerror2("bad positive integer", ms);
351 free(ms);
352
353 if (a == 0) a = new std::vector<int>;
354 a->push_back(m);
355 return a;
356}
357
358int
359MPQCIn::check_string(const char *s)
360{
361 checking_ = 1;
362#ifdef HAVE_SSTREAM
363 istringstream in(s);
364#else
365 istrstream in(s);
366#endif
367 lexer_->switch_streams(&in, &ExEnv::outn());
368 int token;
369 while ((token = ylex())) {
370 if (token == T_OO_INPUT_KEYWORD) return 0;
371 }
372 checking_ = 0;
373 return 1;
374}
375
376char *
377MPQCIn::parse_string(const char *s)
378{
379 // read in easy input
380#ifdef HAVE_SSTREAM
381 istringstream in(s);
382#else
383 istrstream in(s);
384#endif
385 lexer_->switch_streams(&in, &ExEnv::outn());
386 yparse();
387
388 // form the oo input
389#ifdef HAVE_SSTREAM
390 ostringstream ostrs;
391#else
392 ostrstream ostrs;
393#endif
394 SCFormIO::init_ostream(ostrs);
395 ostrs << decindent;
396 if (mol_.null()) error("no molecule given");
397 if (symmetry_.set() && strcmp(symmetry_.val(),"auto") != 0) {
398 mol_->symmetrize(new PointGroup(symmetry_.val()));
399 }
400 ostrs << indent << "molecule<Molecule>: (" << endl;
401 ostrs << incindent;
402 ostrs << indent << "symmetry = "
403 << (symmetry_.set()?symmetry_.val():"auto") << endl;
404 ostrs << indent << "unit = \""
405 << (molecule_bohr_.val()?"bohr":"angstrom")
406 << "\"" << endl;
407 mol_->print_parsedkeyval(ostrs, 0, 0, 0);
408 ostrs << decindent;
409 ostrs << indent << ")" << endl;
410 write_basis_object(ostrs, "basis", basis_.val());
411 ostrs << indent << "mpqc: (" << endl;
412 ostrs << incindent;
413 ostrs << indent << "do_gradient = " << gradient_.val() << endl;
414 ostrs << indent << "optimize = " << optimize_.val() << endl;
415 ostrs << indent << "restart = " << restart_.val() << endl;
416 ostrs << indent << "checkpoint = " << checkpoint_.val() << endl;
417 ostrs << indent << "savestate = " << checkpoint_.val() << endl;
418 bool need_cints = false;
419 write_energy_object(ostrs, "mole", method_.val(), 0, optimize_.val(),
420 need_cints);
421 if (need_cints)
422 ostrs << indent << "integrals<IntegralCints>: ()" << std::endl;
423 if (optimize_.val()) {
424 const char *coortype = "SymmMolecularCoor";
425 if (opt_type_.val() == T_CARTESIAN) coortype = "CartMolecularCoor";
426 else if (redund_coor_.val()) coortype = "RedundMolecularCoor";
427 ostrs << indent << "coor<" << coortype << ">: (" << endl;
428 ostrs << indent << " molecule = $:molecule" << endl;
429 if (opt_type_.val() == T_INTERNAL) {
430 ostrs << indent << " generator<IntCoorGen>: (" << endl;
431 ostrs << indent << " molecule = $:molecule" << endl;
432 ostrs << indent << " )" << endl;
433 }
434 ostrs << indent << ")" << endl;
435 ostrs << indent << "opt<QNewtonOpt>: (" << endl;
436 ostrs << indent << " function = $:mpqc:mole" << endl;
437 ostrs << indent << " update<BFGSUpdate>: ()" << endl;
438 ostrs << indent << " convergence<MolEnergyConvergence>: (" << endl;
439 ostrs << indent << " cartesian = yes" << endl;
440 ostrs << indent << " energy = $:mpqc:mole" << endl;
441 ostrs << indent << " )" << endl;
442 ostrs << indent << ")" << endl;
443 }
444
445 if (frequencies_.val()) {
446 ostrs << indent << "freq<MolecularFrequencies>: (" << endl;
447 ostrs << indent << " molecule = $:molecule" << endl;
448 ostrs << indent << ")" << endl;
449 }
450
451 ostrs << decindent;
452 ostrs << indent << ")" << endl;
453 ostrs << ends;
454
455#ifdef HAVE_SSTREAM
456 int n = 1 + strlen(ostrs.str().c_str());
457 char *in_char_array = strcpy(new char[n],ostrs.str().c_str());
458#else
459 char *in_char_array = ostrs.str();
460#endif
461 return in_char_array;
462}
463
464void
465MPQCIn::write_vector(ostream &ostrs,
466 const char *keyvalname,
467 const char *name, MPQCInDatum<std::vector<int> *>&vec,
468 int require_nirrep)
469{
470 if (vec.set()) {
471 ostrs << indent << keyvalname << " = ";
472 if (!require_nirrep && vec.val()->size() == 1) {
473 ostrs << (*vec.val())[0] << endl;
474 }
475 else if (nirrep_ && vec.val()->size() == nirrep_) {
476 ostrs << "[";
477 for (int i=0; i<nirrep_; i++) {
478 ostrs << " " << (*vec.val())[i];
479 }
480 ostrs << "]" << endl;
481 }
482 else {
483 if (!require_nirrep) error2("need 1 or n_irrep elements in ", name);
484 else {
485 error2("need n_irrep (must give symmetry) elements in ", name);
486 }
487 }
488 }
489}
490
491void
492MPQCIn::write_energy_object(ostream &ostrs,
493 const char *keyword,
494 const char *method,
495 const char *basis,
496 int coor, bool &need_cints)
497{
498 int nelectron = int(mol_->nuclear_charge()+1e-6) - charge_.val();
499 if (nelectron < 0) {
500 error("charge is impossibly large");
501 }
502 if (nelectron%2 == 0 && mult_.val()%2 == 0
503 ||nelectron%2 == 1 && mult_.val()%2 == 1) {
504 error("given multiplicity is not possible");
505 }
506
507 const char *method_object = 0;
508 const char *reference_method = 0;
509 const char *guess_method = method;
510 const char *auxbasis_key = 0;
511 int dft = 0;
512 int uscf = 0;
513 ostringstream o_extra;
514 SCFormIO::init_ostream(o_extra);
515 o_extra << incindent;
516 if (method) {
517 // Hartree-Fock methods
518 if (!strcmp(method, "HF")) {
519 if (mult_.val() == 1) method_object = "CLHF";
520 else { uscf = 1; method_object = "UHF"; }
521 }
522 else if (!strcmp(method, "RHF")) {
523 if (mult_.val() == 1) method_object = "CLHF";
524 else method_object = "HSOSHF";
525 }
526 else if (!strcmp(method, "UHF")) {
527 method_object = "UHF";
528 uscf = 1;
529 }
530 // Density Functional Methods
531 else if (!strcmp(method, "KS")) {
532 guess_method = "HF";
533 if (mult_.val() == 1) method_object = "CLKS";
534 else { uscf = 1; method_object = "UKS"; }
535 dft = 1;
536 }
537 else if (!strcmp(method, "RKS")) {
538 guess_method = "RHF";
539 if (mult_.val() == 1) method_object = "CLKS";
540 else method_object = "HSOSKS";
541 dft = 1;
542 }
543 else if (!strcmp(method, "UKS")) {
544 guess_method = "UHF";
545 method_object = "UKS";
546 dft = 1;
547 uscf = 1;
548 }
549 // Perturbation Theory
550 else if (!strcmp(method, "MP2")) {
551 guess_method = 0;
552 method_object = "MBPT2";
553 reference_method = "HF";
554 if (mult_.val() != 1) {
555 error("MP2 can only be used with multiplicity 1: try ZAPT2");
556 }
557 }
558 // Perturbation Theory
559 else if (!strcmp(method, "MP2-R12/A")) {
560 need_cints = true;
561 auxbasis_key = "aux_basis";
562 guess_method = 0;
563 method_object = "MBPT2_R12";
564 reference_method = "HF";
565 o_extra << indent << "stdapprox = \"A\"" << endl;
566 if (method_absmethod_.val() != 0)
567 o_extra << indent
568 << "abs_method = " << method_absmethod_.val() << endl;
569 if (method_ebc_.val() != 0)
570 o_extra << indent
571 << "ebc = " << method_ebc_.val() << endl;
572 if (method_gbc_.val() != 0)
573 o_extra << indent
574 << "gbc = " << method_gbc_.val() << endl;
575 if (mult_.val() != 1) {
576 error("MP2-R12 can only be used with multiplicity 1");
577 }
578 }
579 // Perturbation Theory
580 else if (!strcmp(method, "MP2-R12/A'")) {
581 need_cints = true;
582 auxbasis_key = "aux_basis";
583 guess_method = 0;
584 method_object = "MBPT2_R12";
585 reference_method = "HF";
586 o_extra << indent << "stdapprox = \"A'\"" << endl;
587 if (method_absmethod_.val() != 0)
588 o_extra << indent
589 << "abs_method = " << method_absmethod_.val() << endl;
590 if (method_ebc_.val() != 0)
591 o_extra << indent
592 << "ebc = " << method_ebc_.val() << endl;
593 if (method_gbc_.val() != 0)
594 o_extra << indent
595 << "gbc = " << method_gbc_.val() << endl;
596 if (mult_.val() != 1) {
597 error("MP2-R12 can only be used with multiplicity 1");
598 }
599 }
600 else if (!strcmp(method, "ZAPT2")) {
601 guess_method = 0;
602 method_object = "MBPT2";
603 reference_method = "RHF";
604 if (mult_.val() == 1) {
605 error("ZAPT2 can only be used with multiplicity != 1: try MP2");
606 }
607 if (optimize_.val() || gradient_.val() || frequencies_.val()) {
608 error("cannot do a gradient or optimization with ZAPT2");
609 }
610 }
611 else error2("invalid method: ", method);
612 }
613 else error("no method given");
614 ostrs << indent << keyword << "<" << method_object << ">: (" << endl;
615 ostrs << incindent;
616 ostrs << o_extra.str();
617 if (auxbasis_key
618 && auxbasis_.val() != 0
619 && strcmp(auxbasis_.val(),basis_.val()) != 0) write_basis_object(ostrs, auxbasis_key, auxbasis_.val());
620 if (need_cints) ostrs << indent << "integrals<IntegralCints>: ()" << endl;
621 ostrs << indent << "total_charge = " << charge_.val() << endl;
622 ostrs << indent << "multiplicity = " << mult_.val() << endl;
623 ostrs << indent << "molecule = $:molecule" << endl;
624 if (memory_.val()) ostrs << indent << "memory = " << memory_.val() << endl;
625 if (!strcmp(keyword, "mole") && !reference_method) {
626 ostrs << indent << "print_npa = 1" << endl;
627 }
628 if (reference_method) {
629 write_vector(ostrs, "nfzc", "frozen_docc", frozen_docc_, 0);
630 write_vector(ostrs, "nfzv", "frozen_uocc", frozen_uocc_, 0);
631 }
632 else {
633 if (uscf && (docc_.set() || socc_.set())) {
634 error("cannot set docc or socc for unrestricted methods"
635 " (use alpha and beta)");
636 }
637 else if (uscf) {
638 write_vector(ostrs, "alpha", "alpha", alpha_, 1);
639 write_vector(ostrs, "beta", "beta", beta_, 1);
640 }
641 else if (alpha_.set() || beta_.set()) {
642 error("cannot set alpha or beta for restricted methods"
643 " (use docc and socc)");
644 }
645 else {
646 write_vector(ostrs, "docc", "docc", docc_, 1);
647 write_vector(ostrs, "socc", "socc", socc_, 1);
648 }
649 }
650 if (coor) ostrs << indent << "coor = $:mpqc:coor" << endl;
651 if (basis) {
652 write_basis_object(ostrs, "basis", basis);
653 }
654 else {
655 ostrs << indent << "basis = $:basis" << endl;
656 }
657 if (dft) {
658 if (method_xc_.set()) {
659 ostrs << indent << "functional<StdDenFunctional>: ( name = \""
660 << method_xc_.val()
661 << "\" )" << endl;
662 }
663 else error("no exchange-correlation functional given");
664 if (method_grid_.set()) {
665 ostrs << indent << "integrator<RadialAngularIntegrator>: ( grid = \""
666 << method_grid_.val()
667 << "\" )" << endl;
668 }
669 }
670 if (dft || (!(basis
671 && !strncmp("STO",basis,3))
672 && !(basis
673 && !strncmp("DZ",basis,2))
674 && strncmp("STO",basis_.val(),3)
675 && guess_method)) {
676 if (frequencies_.val()) {
677 ostrs << indent << "keep_guess_wavefunction = 1" << endl;;
678 }
679 const char *guess_basis;
680 if (need_cints) guess_basis = "DZ (Dunning)";
681 else guess_basis = "STO-3G";
682 write_energy_object(ostrs, "guess_wavefunction",
683 guess_method, guess_basis, 0, need_cints);
684 }
685 if (reference_method) {
686 ostrs << indent << "nfzc = auto" << endl;;
687 write_energy_object(ostrs, "reference",
688 reference_method, 0, 0, need_cints);
689 }
690 ostrs << decindent;
691 ostrs << indent << ")" << endl;
692}
693
694void
695MPQCIn::write_basis_object(ostream &ostrs,
696 const char *keyword,
697 const char *basis)
698{
699 if (!basis) error("no basis given");
700 ostrs << indent << keyword << "<GaussianBasisSet>: (" << endl;
701 ostrs << incindent;
702 ostrs << indent << "molecule = $:molecule" << endl;
703 ostrs << indent << "name = \"" << basis << "\"" << endl;
704 ostrs << decindent;
705 ostrs << indent << ")" << endl;
706}
Note: See TracBrowser for help on using the repository browser.