MPQC 2.3.1

- Main Page
- Related Pages
- Classes
- Files

- File List

functional.h
1 //
2 // functional.h --- definition of the dft functional
3 //
4 // Copyright (C) 1997 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #ifndef _chemistry_qc_dft_functional_h
29 #define _chemistry_qc_dft_functional_h
30
31 #ifdef __GNUC__
32 #pragma interface
33 #endif
34
35 #include <util/state/state.h>
36 #include <math/scmat/vector3.h>
37 #include <chemistry/qc/wfn/wfn.h>
38 #include <chemistry/qc/wfn/density.h>
39
40 namespace sc {
41
43 struct PointInputData {
44 enum {X=BatchElectronDensity::X,
45 Y=BatchElectronDensity::Y,
46 Z=BatchElectronDensity::Z};
47 enum {XX=BatchElectronDensity::XX,
48 YX=BatchElectronDensity::YX,
49 YY=BatchElectronDensity::YY,
50 ZX=BatchElectronDensity::ZX,
51 ZY=BatchElectronDensity::ZY,
52 ZZ=BatchElectronDensity::ZZ};
53 struct SpinData {
54 double rho;
55 // rho^(1/3)
56 double rho_13;
57
58 double del_rho[3];
59 // gamma = (del rho).(del rho)
60 double gamma;
61
62 // hessian of rho
63 double hes_rho[6];
64 // del^2 rho
65 double lap_rho;
66 };
67 SpinData a, b;
68
69 // gamma_ab = (del rho_a).(del rho_b)
70 double gamma_ab;
71
72 const SCVector3 &r;
73
74 // fill in derived quantities
75 void compute_derived(int spin_polarized,
76 int need_gradient,
77 int need_hessian);
78
79 PointInputData(const SCVector3& r_): r(r_) {}
80 };
81
83 struct PointOutputData {
84 // energy at r
85 double energy;
86
87 // derivative of functional wrt density
88 double df_drho_a;
89 double df_drho_b;
90
91 // derivative of functional wrt density gradient
92 double df_dgamma_aa;
93 double df_dgamma_bb;
94 double df_dgamma_ab;
95
96 void zero(){energy=df_drho_a=df_drho_b=df_dgamma_aa=df_dgamma_bb=df_dgamma_ab=0.0;}
97
98 };
99
101 class DenFunctional: virtual public SavableState {
102 protected:
103 int spin_polarized_;
104 int compute_potential_;
105 double a0_; // for ACM functionals
106
107 void do_fd_point(PointInputData&id,double&in,double&out,
108 double lower_bound, double upper_bound);
109 public:
110 DenFunctional();
111 DenFunctional(const Ref<KeyVal> &);
112 DenFunctional(StateIn &);
113 ~DenFunctional();
114 void save_data_state(StateOut &);
115
116 // Set to zero if dens_alpha == dens_beta everywhere.
117 // The default is false.
118 virtual void set_spin_polarized(int i);
119 // Set to nonzero if the potential should be computed.
120 // The default is false.
121 virtual void set_compute_potential(int i);
122
123 // Must return 1 if the density gradient must also be provided.
124 // The default implementation returns 0.
125 virtual int need_density_gradient();
126 // Must return 1 if the density hessian must also be provided.
127 // The default implementation returns 0.
128 virtual int need_density_hessian();
129
130 virtual void point(const PointInputData&, PointOutputData&) = 0;
131 void gradient(const PointInputData&, PointOutputData&,
132 double *gradient, int acenter,
133 GaussianBasisSet *basis,
134 const double *dmat_a, const double *dmat_b,
135 int ncontrib, const int *contrib,
136 int ncontrib_bf, const int *contrib_bf,
137 const double *bs_values, const double *bsg_values,
138 const double *bsh_values);
139
141 virtual double a0() const;
142
143 void fd_point(const PointInputData&, PointOutputData&);
144 int test(const PointInputData &);
145 int test();
146 };
147
148
151 class NElFunctional: public DenFunctional {
152 public:
153 NElFunctional();
154 NElFunctional(const Ref<KeyVal> &);
155 NElFunctional(StateIn &);
156 ~NElFunctional();
157 void save_data_state(StateOut &);
158
159 void point(const PointInputData&, PointOutputData&);
160 };
161
164 class SumDenFunctional: public DenFunctional {
165 protected:
166 int n_;
167 Ref<DenFunctional> *funcs_;
168 double *coefs_;
169 public:
170 SumDenFunctional();
198 SumDenFunctional(const Ref<KeyVal> &);
199 SumDenFunctional(StateIn &);
200 ~SumDenFunctional();
201 void save_data_state(StateOut &);
202
203 void set_spin_polarized(int);
204 void set_compute_potential(int);
205 int need_density_gradient();
206
207 void point(const PointInputData&, PointOutputData&);
208
209 void print(std::ostream& =ExEnv::out0()) const;
210
213 double a0() const;
214 };
215
288 class StdDenFunctional: public SumDenFunctional {
289 protected:
290 char *name_;
291 void init_arrays(int n);
292 public:
293 StdDenFunctional();
298 StdDenFunctional(const Ref<KeyVal> &);
299 StdDenFunctional(StateIn &);
300 ~StdDenFunctional();
301 void save_data_state(StateOut &);
302
303 void print(std::ostream& =ExEnv::out0()) const;
304 };
305
307 class LSDACFunctional: public DenFunctional {
308 protected:
309 public:
310 LSDACFunctional();
311 LSDACFunctional(const Ref<KeyVal> &);
312 LSDACFunctional(StateIn &);
313 ~LSDACFunctional();
314 void save_data_state(StateOut &);
315
316 void point(const PointInputData&, PointOutputData&);
317 virtual
318 void point_lc(const PointInputData&, PointOutputData&,
319 double &ec_local, double &decrs, double &deczeta) = 0;
320
321 };
322
323
332 class PBECFunctional: public DenFunctional {
333 protected:
334 Ref<LSDACFunctional> local_;
335 double gamma;
336 double beta;
337 void init_constants();
338 double rho_deriv(double rho_a, double rho_b, double mdr,
339 double ec_local, double ec_local_dra);
340 double gab_deriv(double rho, double phi, double mdr, double ec_local);
341 public:
342 PBECFunctional();
343 PBECFunctional(const Ref<KeyVal> &);
344 PBECFunctional(StateIn &);
345 ~PBECFunctional();
346 void save_data_state(StateOut &);
347 int need_density_gradient();
348 void point(const PointInputData&, PointOutputData&);
349 void set_spin_polarized(int);
350
351 };
352
363 class PW91CFunctional: public DenFunctional {
364 protected:
365 Ref<LSDACFunctional> local_;
366 double a;
367 double b;
368 double c;
369 double d;
370 double alpha;
371 double c_c0;
372 double c_x;
373 double nu;
374 void init_constants();
375 double limit_df_drhoa(double rhoa, double gamma,
376 double ec, double decdrhoa);
377
378 public:
379 PW91CFunctional();
380 PW91CFunctional(const Ref<KeyVal> &);
381 PW91CFunctional(StateIn &);
382 ~PW91CFunctional();
383 void save_data_state(StateOut &);
384 int need_density_gradient();
385
386 void point(const PointInputData&, PointOutputData&);
387 void set_spin_polarized(int);
388
389 };
390
397 class P86CFunctional: public DenFunctional {
398 protected:
399 double a_;
400 double C1_;
401 double C2_;
402 double C3_;
403 double C4_;
404 double C5_;
405 double C6_;
406 double C7_;
407 void init_constants();
408 public:
409 P86CFunctional();
410 P86CFunctional(const Ref<KeyVal> &);
411 P86CFunctional(StateIn &);
412 ~P86CFunctional();
413 void save_data_state(StateOut &);
414 int need_density_gradient();
415 void point(const PointInputData&, PointOutputData&);
416
417 };
418
419
420 // The Perdew 1986 (P86) Correlation Functional computes energies and densities
421 // using the designated local correlation functional.
422 class NewP86CFunctional: public DenFunctional {
423 protected:
424 double a_;
425 double C1_;
426 double C2_;
427 double C3_;
428 double C4_;
429 double C5_;
430 double C6_;
431 double C7_;
432 void init_constants();
433 double rho_deriv(double rho_a, double rho_b, double mdr);
434 double gab_deriv(double rho_a, double rho_b, double mdr);
435
436 public:
437 NewP86CFunctional();
438 NewP86CFunctional(const Ref<KeyVal> &);
439 NewP86CFunctional(StateIn &);
440 ~NewP86CFunctional();
441 void save_data_state(StateOut &);
442 int need_density_gradient();
443 void point(const PointInputData&, PointOutputData&);
444 };
445
449 class SlaterXFunctional: public DenFunctional {
450 protected:
451 public:
452 SlaterXFunctional();
453 SlaterXFunctional(const Ref<KeyVal> &);
454 SlaterXFunctional(StateIn &);
455 ~SlaterXFunctional();
456 void save_data_state(StateOut &);
457 void point(const PointInputData&, PointOutputData&);
458 };
459
467 class VWNLCFunctional: public LSDACFunctional {
468 protected:
469 double Ap_, Af_, A_alpha_;
470 double x0p_mc_, bp_mc_, cp_mc_, x0f_mc_, bf_mc_, cf_mc_;
471 double x0p_rpa_, bp_rpa_, cp_rpa_, x0f_rpa_, bf_rpa_, cf_rpa_;
472 double x0_alpha_mc_, b_alpha_mc_, c_alpha_mc_;
473 double x0_alpha_rpa_, b_alpha_rpa_, c_alpha_rpa_;
474 void init_constants();
475
476 double F(double x, double A, double x0, double b, double c);
477 double dFdr_s(double x, double A, double x0, double b, double c);
478 public:
479 VWNLCFunctional();
480 VWNLCFunctional(const Ref<KeyVal> &);
481 VWNLCFunctional(StateIn &);
482 ~VWNLCFunctional();
483 void save_data_state(StateOut &);
484
485 virtual
486 void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
487 };
488
491 class VWN1LCFunctional: public VWNLCFunctional {
492 protected:
493 double x0p_, bp_, cp_, x0f_, bf_, cf_;
494 public:
496 VWN1LCFunctional();
498 VWN1LCFunctional(int use_rpa);
504 VWN1LCFunctional(const Ref<KeyVal> &);
505 VWN1LCFunctional(StateIn &);
506 ~VWN1LCFunctional();
507 void save_data_state(StateOut &);
508
509 void point_lc(const PointInputData&, PointOutputData&,
510 double &, double &, double &);
511 };
512
515 class VWN2LCFunctional: public VWNLCFunctional {
516 protected:
517 public:
519 VWN2LCFunctional();
521 VWN2LCFunctional(const Ref<KeyVal> &);
522 VWN2LCFunctional(StateIn &);
523 ~VWN2LCFunctional();
524 void save_data_state(StateOut &);
525
526 void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
527 };
528
529
532 class VWN3LCFunctional: public VWNLCFunctional {
533 protected:
534 int monte_carlo_prefactor_;
535 int monte_carlo_e0_;
536 public:
537 VWN3LCFunctional(int mcp = 1, int mce0 = 1);
538 VWN3LCFunctional(const Ref<KeyVal> &);
539 VWN3LCFunctional(StateIn &);
540 ~VWN3LCFunctional();
541 void save_data_state(StateOut &);
542
543 void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
544 };
545
548 class VWN4LCFunctional: public VWNLCFunctional {
549 protected:
550 int monte_carlo_prefactor_;
551 public:
552 VWN4LCFunctional();
553 VWN4LCFunctional(const Ref<KeyVal> &);
554 VWN4LCFunctional(StateIn &);
555 ~VWN4LCFunctional();
556 void save_data_state(StateOut &);
557
558 void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
559 };
560
563 class VWN5LCFunctional: public VWNLCFunctional {
564 protected:
565 public:
566 VWN5LCFunctional();
567 VWN5LCFunctional(const Ref<KeyVal> &);
568 VWN5LCFunctional(StateIn &);
569 ~VWN5LCFunctional();
570 void save_data_state(StateOut &);
571
572 void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
573 };
574
580 class PW92LCFunctional: public LSDACFunctional {
581 protected:
582 double F(double x, double A, double alpha_1, double beta_1, double beta_2,
583 double beta_3, double beta_4, double p);
584 double dFdr_s(double x, double A, double alpha_1, double beta_1, double beta_2,
585 double beta_3, double beta_4, double p);
586 public:
587 PW92LCFunctional();
588 PW92LCFunctional(const Ref<KeyVal> &);
589 PW92LCFunctional(StateIn &);
590 ~PW92LCFunctional();
591 void save_data_state(StateOut &);
592
593 void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
594 };
595
601 class PZ81LCFunctional: public LSDACFunctional {
602 protected:
603 double Fec_rsgt1(double rs, double beta_1, double beta_2, double gamma);
604 double dFec_rsgt1_drho(double rs, double beta_1, double beta_2, double gamma,
605 double &dec_drs);
606 double Fec_rslt1(double rs, double A, double B, double C, double D);
607 double dFec_rslt1_drho(double rs, double A, double B, double C, double D,
608 double &dec_drs);
609 public:
610 PZ81LCFunctional();
611 PZ81LCFunctional(const Ref<KeyVal> &);
612 PZ81LCFunctional(StateIn &);
613 ~PZ81LCFunctional();
614 void save_data_state(StateOut &);
615
616 void point_lc(const PointInputData&, PointOutputData&, double &, double &, double &);
617 };
618
620 class XalphaFunctional: public DenFunctional {
621 protected:
622 double alpha_;
623 double factor_;
624 public:
625 XalphaFunctional();
626 XalphaFunctional(const Ref<KeyVal> &);
627 XalphaFunctional(StateIn &);
628 ~XalphaFunctional();
629 void save_data_state(StateOut &);
630
631 void point(const PointInputData&, PointOutputData&);
632
633 void print(std::ostream& =ExEnv::out0()) const;
634 };
635
640 class Becke88XFunctional: public DenFunctional {
641 protected:
642 double beta_;
643 double beta6_;
644 public:
645 Becke88XFunctional();
646 Becke88XFunctional(const Ref<KeyVal> &);
647 Becke88XFunctional(StateIn &);
648 ~Becke88XFunctional();
649 void save_data_state(StateOut &);
650
651 int need_density_gradient();
652
653 void point(const PointInputData&, PointOutputData&);
654 };
655
664 class LYPCFunctional: public DenFunctional {
665 protected:
666 double a_;
667 double b_;
668 double c_;
669 double d_;
670 void init_constants();
671 public:
672 LYPCFunctional();
673 LYPCFunctional(const Ref<KeyVal> &);
674 LYPCFunctional(StateIn &);
675 ~LYPCFunctional();
676 void save_data_state(StateOut &);
677
678 int need_density_gradient();
679
680 void point(const PointInputData&, PointOutputData&);
681 };
682
687 class PW86XFunctional: public DenFunctional {
688 protected:
689 double a_;
690 double b_;
691 double c_;
692 double m_;
693 void init_constants();
694 public:
695 PW86XFunctional();
696 PW86XFunctional(const Ref<KeyVal> &);
697 PW86XFunctional(StateIn &);
698 ~PW86XFunctional();
699 void save_data_state(StateOut &);
700
701 int need_density_gradient();
702
703 void point(const PointInputData&, PointOutputData&);
704 };
705
721 class PBEXFunctional: public DenFunctional {
722 protected:
723 double mu;
724 double kappa;
725 void spin_contrib(const PointInputData::SpinData &,
726 double &mpw, double &dmpw_dr, double &dmpw_dg);
727 void init_constants();
728 public:
729 PBEXFunctional();
730 PBEXFunctional(const Ref<KeyVal> &);
731 PBEXFunctional(StateIn &);
732 ~PBEXFunctional();
733 void save_data_state(StateOut &);
734
735 int need_density_gradient();
736
737 void point(const PointInputData&, PointOutputData&);
738 };
739
750 class PW91XFunctional: public DenFunctional {
751 protected:
752 double a;
753 double b;
754 double c;
755 double d;
756 double a_x;
757 void spin_contrib(const PointInputData::SpinData &,
758 double &mpw, double &dmpw_dr, double &dmpw_dg);
759 void init_constants();
760 public:
761 PW91XFunctional();
762 PW91XFunctional(const Ref<KeyVal> &);
763 PW91XFunctional(StateIn &);
764 ~PW91XFunctional();
765 void save_data_state(StateOut &);
766
767 int need_density_gradient();
768
769 void point(const PointInputData&, PointOutputData&);
770 };
771
776 class mPW91XFunctional: public DenFunctional {
777 protected:
778 double b;
779 double beta;
780 double c;
781 double d;
782 double a_x;
783 double x_d_coef;
784
785 void spin_contrib(const PointInputData::SpinData &,
786 double &mpw, double &dmpw_dr, double &dmpw_dg);
787 public:
788 enum Func { B88, PW91, mPW91 };
789
791 mPW91XFunctional();
794 mPW91XFunctional(Func variant);
813 mPW91XFunctional(const Ref<KeyVal> &);
814 mPW91XFunctional(StateIn &);
815 ~mPW91XFunctional();
816 void save_data_state(StateOut &);
817
818 int need_density_gradient();
819
820 void point(const PointInputData&, PointOutputData&);
821
822 void init_constants(Func);
823 };
824
829 class G96XFunctional: public DenFunctional {
830 protected:
831 double b_;
832 void init_constants();
833 public:
834 G96XFunctional();
835 G96XFunctional(const Ref<KeyVal> &);
836 G96XFunctional(StateIn &);
837 ~G96XFunctional();
838 void save_data_state(StateOut &);
839
840 int need_density_gradient();
841
842 void point(const PointInputData&, PointOutputData&);
843 };
844
845 }
846
847 #endif
848
849 // Local Variables:
850 // mode: c++
851 // c-file-style: "CLJ"
852 // End: