OpenRDM
An open-source library for reduced-density matrix-based analysis and computation
mcpdft_solver.h
Go to the documentation of this file.
1 /*
2  * @BEGIN LICENSE
3  *
4  * mcpdft by Psi4 Developer, a plugin to:
5  *
6  * Psi4: an open-source quantum chemistry software package
7  *
8  * Copyright (c) 2007-2016 The Psi4 Developers.
9  *
10  * The copyrights for code used from other parties are included in
11  * the corresponding files.
12  *
13  * This program is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * This program is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License along
24  * with this program; if not, write to the Free Software Foundation, Inc.,
25  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26  *
27  * @END LICENSE
28  */
29 
30 
31 #ifndef MCPDFT_SOLVER_H
32 #define MCPDFT_SOLVER_H
33 
34 #define PSIF_DCC_QMO 268
35 #define PSIF_V2RDM_CHECKPOINT 269
36 #define PSIF_V2RDM_D2AA 270
37 #define PSIF_V2RDM_D2AB 271
38 #define PSIF_V2RDM_D2BB 272
39 #define PSIF_V2RDM_D3AAA 273
40 #define PSIF_V2RDM_D3AAB 274
41 #define PSIF_V2RDM_D3BBA 275
42 #define PSIF_V2RDM_D3BBB 276
43 #define PSIF_V2RDM_D1A 277
44 #define PSIF_V2RDM_D1B 278
45 
46 #include<stdio.h>
47 #include<stdlib.h>
48 #include<math.h>
49 #include<string>
50 
51 #include "psi4/libmints/wavefunction.h"
52 
53 // for reading integrals from disk
54 #include <psi4/libiwl/iwl.h>
55 
56 // for dft
57 #include "psi4/libfock/v.h"
58 #include "psi4/libfunctional/superfunctional.h"
59 
60 // for grid
61 #include "psi4/libfock/points.h"
62 #include "psi4/libfock/cubature.h"
63 
64 #include "psi4/psi4-dec.h"
65 #include <psi4/psifiles.h>
66 #include <psi4/libpsio/psio.hpp>
67 #include <psi4/libpsi4util/PsiOutStream.h>
68 
69 namespace psi{ namespace mcpdft{
70 
71 struct tpdm {
72  int i;
73  int j;
74  int k;
75  int l;
76  double val;
77 };
78 
79 struct opdm {
80  int i;
81  int j;
82  double val;
83 };
84 
85 class MCPDFTSolver: public Wavefunction{
86 
87  public:
88 
89  MCPDFTSolver(std::shared_ptr<psi::Wavefunction> reference_wavefunction,Options & options);
90  ~MCPDFTSolver();
91  void common_init();
92  double compute_energy();
93  void polyradical_analysis();
94  virtual bool same_a_b_orbs() const { return same_a_b_orbs_; }
95  virtual bool same_a_b_dens() const { return same_a_b_dens_; }
96 
98  void print_banner();
99 
100 // protected:
101 
107 
108 
111 
114 
117 
119  std::shared_ptr<VBase> potential_;
120 
122  std::shared_ptr<PointFunctions> points_func_;
123 
126 
129 
131  int * symmetry_;
132 
134  std::vector < std::vector < std::pair<int,int> > > gems_;
135 
136  // geminal / orbital maps
137 
138  int *** bas_;
139  int *** ibas_;
140 
141  // read erfc integrals from disk
142  void ReadAllIntegrals(iwlbuf *Buf);
143 
145  void ReadRangeSeparatedIntegrals();
146 
148  double * erfc_tei_;
149 
152 
155 
158 
161 
164 
167 
169  int deriv_;
170 
172  bool is_gga_;
173 
175  bool is_df_;
176 
178  bool is_meta_;
179 
182 
185 
187  long int phi_points_;
188 
190  int nQ_;
191 
193 // std::shared_ptr<Matrix> super_phi_ao_;
194 //
195 // /// d phi / dx matrix (AO)
196 // std::shared_ptr<Matrix> super_phi_x_ao_;
197 //
198 // /// d phi / dy matrix (AO)
199 // std::shared_ptr<Matrix> super_phi_y_ao_;
200 //
201 // /// d phi / dz matrix (AO)
202 // std::shared_ptr<Matrix> super_phi_z_ao_;
203 
205  std::shared_ptr<Matrix> super_phi_;
206 
208  std::shared_ptr<Matrix> super_phi_x_;
209 
211  std::shared_ptr<Matrix> super_phi_y_;
212 
214  std::shared_ptr<Matrix> super_phi_z_;
215 
217  std::shared_ptr<Matrix> super_phi_xx_;
218 
220  std::shared_ptr<Matrix> super_phi_xy_;
221 
223  std::shared_ptr<Matrix> super_phi_xz_;
224 
226  std::shared_ptr<Matrix> super_phi_yy_;
227 
229  std::shared_ptr<Matrix> super_phi_yz_;
230 
232  std::shared_ptr<Matrix> super_phi_zz_;
233 
235  std::shared_ptr<Matrix> super_gamma_aa_;
236 
238  std::shared_ptr<Matrix> super_gamma_ab_;
239 
241  std::shared_ptr<Matrix> super_gamma_bb_;
242 
244  std::shared_ptr<Matrix> super_tau_a_;
245 
247  std::shared_ptr<Matrix> super_tau_b_;
248 
250  std::shared_ptr<Vector> grid_x_;
251 
253  std::shared_ptr<Vector> grid_y_;
254 
256  std::shared_ptr<Vector> grid_z_;
257 
259  std::shared_ptr<Vector> grid_w_;
260 
262  std::shared_ptr<Vector> sigma_aa_;
263 
265  std::shared_ptr<Vector> sigma_bb_;
266 
268  std::shared_ptr<Vector> sigma_ab_;
269 
271  void GetGridInfo();
272 
274  void BuildPhiMatrixAO(std::string phi_type, std::shared_ptr<Matrix> myphi);
275 
277  void TransformPhiMatrixAOMO(std::shared_ptr<Matrix> phi_in, std::shared_ptr<Matrix> phi_out);
278 
280  void PrintTPDM(double* D);
281 
283  void ReadTPDM();
284 
286  void ReadOPDM();
287 
289  void ReadCITPDM(double* D, const char* fileName);
290 
292  void ReadCIOPDM(std::shared_ptr<Matrix> D, const char* fileName);
293 
295  std::vector< std::shared_ptr<Matrix> > BuildJK();
296 
299  double RangeSeparatedTEE(std::string range_separation_type);
300 
302  std::shared_ptr<Vector> rho_a_;
303 
305  std::shared_ptr<Vector> rho_b_;
306 
308  std::shared_ptr<Vector> rho_;
309 
311  std::shared_ptr<Vector> rho_a_x_;
312 
314  std::shared_ptr<Vector> rho_b_x_;
315 
317  std::shared_ptr<Vector> rho_a_y_;
318 
320  std::shared_ptr<Vector> rho_b_y_;
321 
323  std::shared_ptr<Vector> rho_a_z_;
324 
326  std::shared_ptr<Vector> rho_b_z_;
327 
329  std::shared_ptr<Vector> tr_rho_a_;
330 
332  std::shared_ptr<Vector> tr_rho_b_;
333 
335  std::shared_ptr<Vector> tr_rho_a_x_;
336 
338  std::shared_ptr<Vector> tr_rho_b_x_;
339 
341  std::shared_ptr<Vector> tr_rho_a_y_;
342 
344  std::shared_ptr<Vector> tr_rho_b_y_;
345 
347  std::shared_ptr<Vector> tr_rho_a_z_;
348 
350  std::shared_ptr<Vector> tr_rho_b_z_;
351 
353  std::shared_ptr<Vector> tr_sigma_aa_;
354 
356  std::shared_ptr<Vector> tr_sigma_bb_;
357 
359  std::shared_ptr<Vector> tr_sigma_ab_;
360 
362  std::shared_ptr<Vector> tr_sigma_;
363 
365  std::shared_ptr<Vector> R_;
366 
368  std::shared_ptr<Vector> zeta_;
369 
371  std::shared_ptr<Vector> m_;
372 
374  std::shared_ptr<Vector> tr_m_;
375 
377  std::shared_ptr<Vector> tr_zeta_;
378 
380  std::shared_ptr<Vector> rs_;
381 
383  std::shared_ptr<Vector> tr_rs_;
384 
386  std::shared_ptr<Vector> tau_a_;
387 
389  std::shared_ptr<Vector> tau_b_;
390 
392  std::shared_ptr<Vector> pi_;
393 
395  std::shared_ptr<Vector> pi_x_;
396 
398  std::shared_ptr<Vector> pi_y_;
399 
401  std::shared_ptr<Vector> pi_z_;
402 
404  void BuildRho();
405 
407  void BuildRhoFast(int na, int nb);
408 
410  void BuildPi(double * D2ab);
411 
413  void BuildPiFast(tpdm * D2ab, int nab);
414 
416  void BuildPiLowMemory(tpdm * D2ab, int nab);
417 
419  void BuildExchangeCorrelationHole(int p, tpdm * D2ab, int nab, tpdm * D2aa, int naa, tpdm * D2bb, int nbb);
420 
422  std::shared_ptr<Vector> Dr_;
423 
425  std::shared_ptr<Vector> Ur_;
426 
428  std::shared_ptr<Vector> Nr_;
429 
431  void Build_Grad_Pi();
432 
434  void Build_R();
435 
437  void Translate();
438 
440  void Fully_Translate();
441 
443  double Gfunction(double r, double A, double a1, double b1, double b2, double b3, double b4, double p);
444 
445  //############################################################
446  //############ Exchange functions' declarations ##############
447  //############################################################
448 
450  double EX_LDA(std::shared_ptr<Vector> rho_a, std::shared_ptr<Vector> rho_b);
451 
453  double EX_LSDA_Sigma(std::shared_ptr<Vector> rho_sigma);
454 
456  double EX_LSDA(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> ZETA);
457 
459  double EX_LSDA(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B);
460 
462  double EX_B86_MGC();
463 
465  double EX_B88(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);
466  double EX_B88_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);
467 
469  double EX_PBE(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);
470  double EX_PBE_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);
471 
473  double EX_wPBE_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);
474 
476  double EX_RPBE();
477 
479  double EX_UPBE();
480 
481  //############################################################
482  //########### Correlation functions' declarations ############
483  //############################################################
484 
486  double EC_B88_OP(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);
487 
489  double EC_LYP_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B,
490  std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_AB, std::shared_ptr<Vector> SIGMA_BB);
491 
493  double EC_PBE(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B,
494  std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_AB, std::shared_ptr<Vector> SIGMA_BB);
495  double EC_PBE_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B,
496  std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_AB, std::shared_ptr<Vector> SIGMA_BB);
497 
499  double EC_VWN3_RPA(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> ZETA, std::shared_ptr<Vector> RS);
500  double EC_VWN3_RPA_III(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B);
501  // double EC_VWN3_RPA_III(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> ZETTA);
502 
504  double EC_PW92_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B);
505 
506 };
507 
508 }} // end of namespaces
509 
510 #endif
std::shared_ptr< Vector > grid_x_
grid x values
Definition: mcpdft_solver.h:250
std::shared_ptr< Vector > Ur_
effective unpaired electron density function (Head-Gordon)
Definition: mcpdft_solver.h:425
double lr_ex_energy_
Definition: mcpdft_solver.h:105
std::shared_ptr< Matrix > super_gamma_bb_
gamma_bb matrix
Definition: mcpdft_solver.h:241
std::shared_ptr< Matrix > super_phi_xx_
d2 phi / dx2 matrix
Definition: mcpdft_solver.h:217
int k
Definition: mcpdft_solver.h:74
int deriv_
level of differentiation
Definition: mcpdft_solver.h:169
virtual bool same_a_b_orbs() const
Definition: mcpdft_solver.h:94
bool is_unpolarized_
is unpolarized (restricted)?
Definition: mcpdft_solver.h:184
std::shared_ptr< Vector > tr_sigma_aa_
translated inner product of alpha density gradient with itself (gamma_aa)
Definition: mcpdft_solver.h:353
std::shared_ptr< Vector > rho_a_z_
alpha-spin density gradient z
Definition: mcpdft_solver.h:323
int * symmetry_
orbital symmetry
Definition: mcpdft_solver.h:131
std::shared_ptr< Vector > tr_rho_b_z_
translated beta-spin density gradient z
Definition: mcpdft_solver.h:350
std::shared_ptr< Vector > rs_
effective radius of density
Definition: mcpdft_solver.h:380
std::shared_ptr< Matrix > super_gamma_aa_
gamma_aa matrix
Definition: mcpdft_solver.h:235
std::shared_ptr< Matrix > super_tau_a_
tau_a matrix
Definition: mcpdft_solver.h:244
std::shared_ptr< Vector > Dr_
effective unpaired electron density function (Yamaguchi et al.)
Definition: mcpdft_solver.h:422
std::shared_ptr< Matrix > super_phi_
phi matrix (AO)
Definition: mcpdft_solver.h:205
double val
Definition: mcpdft_solver.h:76
std::shared_ptr< Vector > pi_y_
y-component of the gradient of the on-top pair density
Definition: mcpdft_solver.h:398
double coulomb_energy_
======================
Definition: mcpdft_solver.h:103
int l
Definition: mcpdft_solver.h:75
Definition: mcpdft_solver.h:71
bool is_gga_
is gga?
Definition: mcpdft_solver.h:172
std::shared_ptr< Vector > rho_b_y_
beta-spin density gradient y
Definition: mcpdft_solver.h:320
std::shared_ptr< Matrix > super_phi_yy_
d2 phi / dy2 matrix
Definition: mcpdft_solver.h:226
std::shared_ptr< Vector > rho_b_z_
beta-spin density gradient z
Definition: mcpdft_solver.h:326
std::vector< std::vector< std::pair< int, int > > > gems_
geminals, by symmetry
Definition: mcpdft_solver.h:134
double sr_Vee_energy_
short-range part of the two-electron reference energy
Definition: mcpdft_solver.h:160
double val
Definition: mcpdft_solver.h:82
Definition: mcpdft_solver.h:79
std::shared_ptr< Vector > tau_b_
tau kinetic energy of beta electrons
Definition: mcpdft_solver.h:389
std::shared_ptr< Vector > tr_m_
translated spin magnetization density
Definition: mcpdft_solver.h:374
Definition: mcpdft_solver.h:85
std::shared_ptr< Vector > Nr_
effective unpaired electron density function (Head-Gordon&#39;s quartic proposal)
Definition: mcpdft_solver.h:428
std::shared_ptr< Vector > tr_rho_a_
translated alpha-spin density
Definition: mcpdft_solver.h:329
double * erfc_tei_
erfc integrals
Definition: mcpdft_solver.h:148
virtual bool same_a_b_dens() const
Definition: mcpdft_solver.h:95
std::shared_ptr< Vector > sigma_bb_
inner product of beta density gradient with itself (gamma_bb)
Definition: mcpdft_solver.h:265
std::shared_ptr< Vector > rho_
total density
Definition: mcpdft_solver.h:308
Definition: functional.cc:5
long int available_memory_
available memory
Definition: mcpdft_solver.h:181
Definition: blas.h:40
std::shared_ptr< Matrix > super_phi_z_
d phi / dz matrix (MO/NO)
Definition: mcpdft_solver.h:214
std::shared_ptr< Vector > tr_zeta_
translated zeta factor
Definition: mcpdft_solver.h:377
int *** ibas_
Definition: mcpdft_solver.h:139
std::shared_ptr< Vector > pi_z_
z-component of the gradient of the on-top pair density
Definition: mcpdft_solver.h:401
std::shared_ptr< Matrix > super_phi_yz_
d2 phi / dydz matrix
Definition: mcpdft_solver.h:229
std::shared_ptr< Vector > tr_rho_b_
translated beta-spin density
Definition: mcpdft_solver.h:332
std::shared_ptr< Vector > m_
spin magnetization density
Definition: mcpdft_solver.h:371
std::shared_ptr< Vector > pi_x_
x-component of the gradient of the on-top pair density
Definition: mcpdft_solver.h:395
double mp2_corr_energy_
mp2 correlation energy for double-hybrids
Definition: mcpdft_solver.h:163
std::shared_ptr< VBase > potential_
dft potential object
Definition: mcpdft_solver.h:119
std::shared_ptr< Vector > tr_sigma_ab_
translated inner product of alpha density gradient with beta density gradient (gamma_ab) ...
Definition: mcpdft_solver.h:359
std::shared_ptr< Vector > tr_rho_b_y_
translated beta-spin density gradient y
Definition: mcpdft_solver.h:344
bool is_meta_
is meta?
Definition: mcpdft_solver.h:178
std::shared_ptr< Vector > zeta_
zeta factor zeta = ( rho_a(r) - rho_b(r) ) / ( rho_a(r) + rho_b(r) )
Definition: mcpdft_solver.h:368
int nQ_
number of auxilliary basis functions
Definition: mcpdft_solver.h:190
std::shared_ptr< Matrix > super_gamma_ab_
gamma_ab matrix
Definition: mcpdft_solver.h:238
double lr_Vee_energy_
long-range part of the two-electron reference energy
Definition: mcpdft_solver.h:157
std::shared_ptr< Vector > sigma_aa_
inner product of alpha density gradient with itself (gamma_aa)
Definition: mcpdft_solver.h:262
double two_electron_energy_
two-electron part of the reference energy
Definition: mcpdft_solver.h:154
double reference_energy_
reference energy
Definition: mcpdft_solver.h:151
std::shared_ptr< Vector > tr_rho_a_z_
translated alpha-spin density gradient z
Definition: mcpdft_solver.h:347
std::shared_ptr< Vector > pi_
the on-top pair density
Definition: mcpdft_solver.h:392
double hf_ex_energy_
Definition: mcpdft_solver.h:104
std::shared_ptr< Vector > rho_a_
alpha-spin density
Definition: mcpdft_solver.h:302
int * pitzer_offset_
offset for orbitals in each irrep
Definition: mcpdft_solver.h:166
std::shared_ptr< PointFunctions > points_func_
points function
Definition: mcpdft_solver.h:122
std::shared_ptr< Vector > grid_w_
grid weights
Definition: mcpdft_solver.h:259
std::shared_ptr< Vector > rho_b_x_
beta-spin density gradient x
Definition: mcpdft_solver.h:314
std::shared_ptr< Matrix > super_phi_y_
d phi / dy matrix (MO/NO)
Definition: mcpdft_solver.h:211
std::shared_ptr< Matrix > super_phi_x_
d phi / dx matrix (MO/NO)
Definition: mcpdft_solver.h:208
std::shared_ptr< Vector > tr_rs_
translated effective radius of density
Definition: mcpdft_solver.h:383
std::shared_ptr< Vector > rho_b_
beta-spin density
Definition: mcpdft_solver.h:305
std::shared_ptr< Matrix > super_tau_b_
tau_b matrix
Definition: mcpdft_solver.h:247
std::shared_ptr< Matrix > super_phi_zz_
d2 phi / dz2 matrix
Definition: mcpdft_solver.h:232
int j
Definition: mcpdft_solver.h:81
std::shared_ptr< Vector > R_
R(r) = 4 * Pi(r) / rho(r) ^ 2.
Definition: mcpdft_solver.h:365
std::shared_ptr< Vector > grid_y_
grid y values
Definition: mcpdft_solver.h:253
std::shared_ptr< Vector > tr_rho_a_x_
translated alpha-spin density gradient x
Definition: mcpdft_solver.h:335
std::shared_ptr< Vector > tr_sigma_bb_
translated inner product of beta density gradient with itself (gamma_bb)
Definition: mcpdft_solver.h:356
std::shared_ptr< Vector > rho_a_y_
alpha-spin density gradient y
Definition: mcpdft_solver.h:317
std::shared_ptr< Vector > sigma_ab_
inner product of alpha density gradient with beta density gradient (gamma_ab)
Definition: mcpdft_solver.h:268
std::shared_ptr< Matrix > super_phi_xy_
d phi / dxdy matrix
Definition: mcpdft_solver.h:220
int i
Definition: mcpdft_solver.h:80
std::shared_ptr< Vector > tr_rho_a_y_
translated alpha-spin density gradient y
Definition: mcpdft_solver.h:341
int *** bas_
Definition: mcpdft_solver.h:138
int max_functions_
maximum number of functions in a block of the phi matrix
Definition: mcpdft_solver.h:125
int i
Definition: mcpdft_solver.h:72
std::shared_ptr< Vector > tr_sigma_
translated inner product of total density gradient
Definition: mcpdft_solver.h:362
std::shared_ptr< Vector > grid_z_
grid z values
Definition: mcpdft_solver.h:256
std::shared_ptr< Vector > rho_a_x_
alpha-spin density gradient x
Definition: mcpdft_solver.h:311
std::shared_ptr< Vector > tau_a_
tau kinetic energy of alpha electrons
Definition: mcpdft_solver.h:386
opdm * opdm_a_
nonzero elements of alpha opdm
Definition: mcpdft_solver.h:113
opdm * opdm_b_
nonzero elements of beta opdm
Definition: mcpdft_solver.h:116
int max_points_
maximum number of grid points in a block of the phi matrix
Definition: mcpdft_solver.h:128
int j
Definition: mcpdft_solver.h:73
bool is_low_memory_
======================
Definition: mcpdft_solver.h:110
bool is_df_
is df?
Definition: mcpdft_solver.h:175
long int phi_points_
number of grid points_
Definition: mcpdft_solver.h:187
std::shared_ptr< Matrix > super_phi_xz_
d2 phi / dxdz matrix
Definition: mcpdft_solver.h:223
std::shared_ptr< Vector > tr_rho_b_x_
translated beta-spin density gradient x
Definition: mcpdft_solver.h:338