Skip to content

Commit 63e402b

Browse files
lcyyorklcyyork
authoredSep 5, 2024
Add genuine analytic DF-MCSCF gradients (evangelistalab#405)
* df algorithm for tei in df-mcscf * d3 dgemm seems work * add nonfrozen df-mcscf gradient * change MCSCF_FREEZE_CORE to MCSCF_IGNORE_FROZEN_ORBS * add frozen-core mcscf gradient; not yet working * frozen-core mcscf * clean up * consider zero dim * fix another zero dim * skip mcscf grad for post-mcscf grad * fix last commit * minor change to copy data * update doc * fix logic of semicanonical mix inactive * oops
1 parent 4d3fe8e commit 63e402b

File tree

28 files changed

+3408
-168
lines changed

28 files changed

+3408
-168
lines changed
 

‎docs/source/methods/mcscf.rst

+31-15
Original file line numberDiff line numberDiff line change
@@ -257,7 +257,7 @@ computation with two doubly occupied orbitals frozen ::
257257

258258

259259
To override this behavior and freeze the core orbitals in the MCSCF procedure, the user can set the option
260-
:code:`MCSCF_FREEZE_CORE` to :code:`True`.
260+
:code:`MCSCF_IGNORE_FROZEN_ORBS` to :code:`False`.
261261
This might be necessary in certain cases, if orbital rotations involving core orbitals cause convergence issues.
262262
The following input (see :code:`tests/manual/mcscf-5/input.dat`) performs an MCSCF calculation on CO molecule and
263263
freezes the 1s orbital of the carbon and oxygen atoms
@@ -277,7 +277,7 @@ freezes the 1s orbital of the carbon and oxygen atoms
277277
mcscf_e_convergence 8
278278
mcscf_g_convergence 6
279279
mcscf_micro_maxiter 4
280-
mcscf_freeze_core true # enables freezing the MCSCF core orbitals
280+
mcscf_ignore_frozen_orbs false # enables freezing the MCSCF core orbitals
281281
}
282282

283283
energy('forte')
@@ -321,17 +321,17 @@ The maximum number of macro iterations.
321321

322322
**MCSCF_MICRO_MAXITER**
323323

324-
The maximum number of micro iterations.
324+
The maximum number of micro iterations (orbital optimization) for a given CI.
325325

326326
* Type: int
327-
* Default: 40
327+
* Default: 6
328328

329-
**MCSCF_MICRO_MINITER**
329+
**MCSCF_MCI_MAXITER**
330330

331-
The minimum number of micro iterations.
331+
The maximum number of micro CI iterations in every macro iteration.
332332

333333
* Type: int
334-
* Default: 6
334+
* Default: 12
335335

336336
**MCSCF_E_CONVERGENCE**
337337

@@ -389,14 +389,6 @@ For example, 1 means do DIIS every iteration and 2 is for every other iteration,
389389
* Type: int
390390
* Default: 1
391391

392-
**MCSCF_CI_SOLVER**
393-
394-
Which active space solver to be used.
395-
396-
* Type: string
397-
* Options: CAS, FCI, ACI, PCI
398-
* Default: CAS
399-
400392
**MCSCF_DEBUG_PRINTING**
401393

402394
Whether to enable debug printing.
@@ -426,6 +418,14 @@ Stop Forte if MCSCF did not converge.
426418
* Type: Boolean
427419
* Default: True
428420

421+
**MCSCF_IGNORE_FROZEN_ORBS**
422+
423+
Whether to ignore frozen orbitals in the input file or not.
424+
By default, all orbitals will be optimized in MCSCF.
425+
426+
* Type: Boolean
427+
* Default: True
428+
429429
Expert Options
430430
~~~~~~~~~~~~~~~
431431

@@ -459,6 +459,22 @@ This option is useful when doing core-excited state computations.
459459
* Type: array
460460
* Default: No Default
461461

462+
**MCSCF_ORB_ORTHO_TRANS**
463+
464+
Ways to compute the orthogonal transformation U from orbital rotation R
465+
466+
* Type: string
467+
* Options: CAYLEY, POWER, PADE
468+
* Default: CAYLEY
469+
470+
**MCSCF_DF_TEIALG**
471+
472+
Algorithm to build (pu|xy) integrals in DF-MCSCF.
473+
Use (Q|pu) if True; Otherwise use JK build for every xy pair
474+
475+
* Type: Boolean
476+
* Fefault: True
477+
462478
CPSCF Options
463479
~~~~~~~~~~~~~
464480

‎forte/helpers/helpers.cc

+4
Original file line numberDiff line numberDiff line change
@@ -149,6 +149,10 @@ std::pair<double, std::string> to_xb(size_t nele, size_t type_size) {
149149
}
150150

151151
void matrix_transpose_in_place(std::vector<double>& data, const size_t m, const size_t n) {
152+
matrix_transpose_in_place(data.data(), m, n);
153+
}
154+
155+
void matrix_transpose_in_place(double* data, const size_t m, const size_t n) {
152156
int nthreads = std::min(omp_get_max_threads(), int(m > n ? n : m));
153157
std::vector<double> tmp(nthreads * (m > n ? m : n));
154158
auto tmp_begin = tmp.begin();

‎forte/helpers/helpers.h

+1
Original file line numberDiff line numberDiff line change
@@ -197,6 +197,7 @@ void apply_permutation_in_place(std::vector<T>& vec, const std::vector<std::size
197197
* Also see https://github.com/bryancatanzaro/inplace
198198
*/
199199
void matrix_transpose_in_place(std::vector<double>& data, const size_t m, const size_t n);
200+
void matrix_transpose_in_place(double* data, const size_t m, const size_t n);
200201

201202
void push_to_psi4_env_globals(double value, const std::string& label);
202203

‎forte/integrals/integrals.cc

+21
Original file line numberDiff line numberDiff line change
@@ -488,6 +488,27 @@ void ForteIntegrals::print_ints() {
488488
}
489489
}
490490

491+
std::shared_ptr<psi::Matrix> ForteIntegrals::Ca_SO2AO(std::shared_ptr<psi::Matrix> Ca_SO) const {
492+
auto aotoso = wfn_->aotoso();
493+
auto nao = nso_;
494+
495+
auto Ca_ao = std::make_shared<psi::Matrix>("Ca_AO", nao, nmo_);
496+
497+
// Transform from the SO to the AO basis
498+
for (int h = 0, index = 0; h < nirrep_; ++h) {
499+
auto nso = nsopi_[h];
500+
if (nso == 0)
501+
continue;
502+
503+
for (int i = 0, nmo_this = nmopi_[h]; i < nmo_this; ++i) {
504+
// notes: LDA value is nso (not nao, see libqt/blas_intfc23.cc)
505+
C_DGEMV('N', nao, nso, 1.0, aotoso->pointer(h)[0], nso, &Ca_SO->pointer(h)[0][i],
506+
nmo_this, 0.0, &Ca_ao->pointer()[0][index++], nmo_);
507+
}
508+
}
509+
return Ca_ao;
510+
}
511+
491512
void ForteIntegrals::rotate_orbitals(std::shared_ptr<psi::Matrix> Ua,
492513
std::shared_ptr<psi::Matrix> Ub, bool re_transform) {
493514
// 1. Rotate the orbital coefficients and store them in the ForteIntegral object

‎forte/integrals/integrals.h

+2
Original file line numberDiff line numberDiff line change
@@ -391,6 +391,8 @@ class ForteIntegrals {
391391

392392
/// Orbital coefficients in AO x MO basis where MO is Pitzer order
393393
virtual std::shared_ptr<psi::Matrix> Ca_AO() const = 0;
394+
/// Transform SO orbital coefficients to AO x MO basis where MO is Pitzer order
395+
std::shared_ptr<psi::Matrix> Ca_SO2AO(std::shared_ptr<psi::Matrix> Ca_SO) const;
394396

395397
/// Obtain AO dipole integrals [X, Y, Z]
396398
/// Each direction is a std::shared_ptr<psi::Matrix> of dimension nao * nao

‎forte/integrals/integrals_psi4_interface.cc

+1-20
Original file line numberDiff line numberDiff line change
@@ -376,26 +376,7 @@ void Psi4Integrals::rotate_mos() {
376376
wfn_->epsilon_b()->copy(eps_a_new);
377377
}
378378

379-
std::shared_ptr<psi::Matrix> Psi4Integrals::Ca_AO() const {
380-
auto aotoso = wfn_->aotoso();
381-
auto nao = nso_;
382-
383-
auto Ca_ao = std::make_shared<psi::Matrix>("Ca_AO", nao, nmo_);
384-
385-
// Transform from the SO to the AO basis
386-
for (int h = 0, index = 0; h < nirrep_; ++h) {
387-
auto nso = nsopi_[h];
388-
if (nso == 0)
389-
continue;
390-
391-
for (int i = 0, nmo_this = nmopi_[h]; i < nmo_this; ++i) {
392-
// notes: LDA value is nso (not nao, see libqt/blas_intfc23.cc)
393-
C_DGEMV('N', nao, nso, 1.0, aotoso->pointer(h)[0], nso, &Ca_->pointer(h)[0][i],
394-
nmo_this, 0.0, &Ca_ao->pointer()[0][index++], nmo_);
395-
}
396-
}
397-
return Ca_ao;
398-
}
379+
std::shared_ptr<psi::Matrix> Psi4Integrals::Ca_AO() const { return Ca_SO2AO(Ca_); }
399380

400381
void Psi4Integrals::build_multipole_ints_ao() {
401382
std::shared_ptr<psi::BasisSet> basisset = wfn_->basisset();

‎forte/mcscf/mcscf_2step.cc

+10-11
Original file line numberDiff line numberDiff line change
@@ -168,8 +168,8 @@ double MCSCF_2STEP::compute_energy() {
168168
};
169169

170170
// prepare for orbital gradients
171-
const bool freeze_core = options_->get_bool("MCSCF_FREEZE_CORE");
172-
MCSCF_ORB_GRAD cas_grad(options_, mo_space_info_, ints_, freeze_core);
171+
const bool ignore_frozen = options_->get_bool("MCSCF_IGNORE_FROZEN_ORBS");
172+
MCSCF_ORB_GRAD cas_grad(options_, mo_space_info_, ints_, ignore_frozen);
173173
auto nrot = cas_grad.nrot();
174174
auto dG = std::make_shared<psi::Vector>("dG", nrot);
175175

@@ -205,7 +205,7 @@ double MCSCF_2STEP::compute_energy() {
205205
if (no_orb_opt) {
206206
energy_ = e_c;
207207
pass_energy_to_psi4();
208-
if (der_type_ == "FIRST") {
208+
if (der_type_ == "FIRST" and options_->get_str("CORRELATION_SOLVER") == "NONE") {
209209
cas_grad.compute_nuclear_gradient();
210210
}
211211
return energy_;
@@ -441,14 +441,13 @@ double MCSCF_2STEP::compute_energy() {
441441
auto F = cas_grad.fock(rdms);
442442
ints_->set_fock_matrix(F, F);
443443

444-
auto inactive_mix = options_->get_bool("SEMI_CANONICAL_MIX_INACTIVE");
445-
auto active_mix = options_->get_bool("SEMI_CANONICAL_MIX_ACTIVE");
444+
// if we do not freeze orbitals, we need to set the inactive_mix flag to make sure
445+
// the frozen and non-frozen core/virtual orbitals are canonicalized together.
446+
auto inactive_mix = ignore_frozen;
446447

447-
// if we do not freeze the core, we need to set the inactive_mix flag to make sure
448-
// the core orbitals are canonicalized together with the active orbitals
449-
if (not freeze_core) {
450-
inactive_mix = true;
451-
}
448+
if (!ignore_frozen)
449+
inactive_mix = options_->get_bool("SEMI_CANONICAL_MIX_INACTIVE");
450+
auto active_mix = options_->get_bool("SEMI_CANONICAL_MIX_ACTIVE");
452451

453452
psi::outfile->Printf("\n Canonicalizing final MCSCF orbitals");
454453
SemiCanonical semi(mo_space_info_, ints_, options_, inactive_mix, active_mix);
@@ -467,7 +466,7 @@ double MCSCF_2STEP::compute_energy() {
467466
throw_convergence_error();
468467

469468
// for nuclear gradient
470-
if (der_type_ == "FIRST") {
469+
if (der_type_ == "FIRST" and options_->get_str("CORRELATION_SOLVER") == "NONE") {
471470
// TODO: remove this re-diagonalization if CI transformation is impelementd
472471
if (not is_single_reference()) {
473472
diagonalize_hamiltonian(

‎forte/mcscf/mcscf_orb_grad.cc

+116-64
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,7 @@
4242
#include "psi4/libpsi4util/process.h"
4343
#include "psi4/libiwl/iwl.hpp"
4444
#include "psi4/libpsio/psio.hpp"
45+
#include "psi4/lib3index/dfhelper.h"
4546

4647
#include "helpers/helpers.h"
4748
#include "helpers/lbfgs/lbfgs.h"
@@ -59,9 +60,9 @@ using namespace ambit;
5960
namespace forte {
6061

6162
MCSCF_ORB_GRAD::MCSCF_ORB_GRAD(std::shared_ptr<ForteOptions> options,
62-
std::shared_ptr<MOSpaceInfo> mo_space_info,
63-
std::shared_ptr<ForteIntegrals> ints, bool freeze_core)
64-
: options_(options), mo_space_info_(mo_space_info), ints_(ints), freeze_core_(freeze_core) {
63+
std::shared_ptr<MOSpaceInfo> mo_space_info,
64+
std::shared_ptr<ForteIntegrals> ints, bool ignore_frozen)
65+
: options_(options), mo_space_info_(mo_space_info), ints_(ints), ignore_frozen_(ignore_frozen) {
6566
startup();
6667
}
6768

@@ -77,6 +78,33 @@ void MCSCF_ORB_GRAD::startup() {
7778

7879
// setup JK
7980
JK_ = ints_->jk();
81+
if (ints_->integral_type() == IntegralType::DF or
82+
ints_->integral_type() == IntegralType::DiskDF) {
83+
tei_alg_ = options_->get_bool("MCSCF_DF_TEIALG") ? TEIALG::DF : TEIALG::JK;
84+
85+
outfile->Printf("\n\n DF-MCSCF adopts integrals from DFHelper of Psi4.\n");
86+
87+
auto wfn = ints_->wfn();
88+
auto bs = wfn->basisset();
89+
auto bs_aux = wfn->get_basisset("DF_BASIS_SCF");
90+
91+
df_helper_ = std::make_shared<psi::DFHelper>(bs, bs_aux);
92+
size_t mem_sys = psi::Process::environment.get_memory() / sizeof(double);
93+
int64_t mem = mem_sys - JK_->memory_estimate();
94+
if (mem < 0) {
95+
auto xb = to_xb(static_cast<size_t>(-mem), sizeof(double));
96+
std::string msg = "Not enough memory! Need at least ";
97+
msg += std::to_string(xb.first) + " " + xb.second.c_str() + " more.";
98+
outfile->Printf("\n %s", msg.c_str());
99+
throw std::runtime_error(msg);
100+
}
101+
df_helper_->set_schwarz_cutoff(options_->get_double("INTS_TOLERANCE"));
102+
df_helper_->set_fitting_condition(options_->get_double("DF_FITTING_CONDITION"));
103+
df_helper_->set_memory(static_cast<size_t>(mem));
104+
df_helper_->set_nthreads(omp_get_max_threads());
105+
df_helper_->set_print_lvl(0);
106+
df_helper_->initialize();
107+
}
80108

81109
// allocate memory for tensors and matrices
82110
init_tensors();
@@ -86,62 +114,49 @@ void MCSCF_ORB_GRAD::startup() {
86114
}
87115

88116
void MCSCF_ORB_GRAD::setup_mos() {
117+
label_to_mos_.clear();
118+
label_to_cmos_.clear();
89119

90120
nirrep_ = mo_space_info_->nirrep();
91121

92122
nsopi_ = ints_->nsopi();
93123
nmopi_ = mo_space_info_->dimension("ALL");
94-
95-
if (freeze_core_) {
96-
ncmopi_ = mo_space_info_->dimension("CORRELATED");
97-
} else {
98-
ncmopi_ =
99-
mo_space_info_->dimension("FROZEN_DOCC") + mo_space_info_->dimension("CORRELATED");
100-
}
101-
102124
ndoccpi_ = mo_space_info_->dimension("INACTIVE_DOCC");
103-
nfrzvpi_ = mo_space_info_->dimension("FROZEN_UOCC");
125+
nactvpi_ = mo_space_info_->dimension("ACTIVE");
126+
actv_mos_ = mo_space_info_->absolute_mo("ACTIVE");
104127

105-
if (freeze_core_) {
106-
nfrzcpi_ = mo_space_info_->dimension("FROZEN_DOCC");
107-
core_mos_ = mo_space_info_->absolute_mo("RESTRICTED_DOCC");
108-
} else {
128+
if (ignore_frozen_) {
129+
ncmopi_ = mo_space_info_->dimension("ALL");
130+
nfrzvpi_ = psi::Dimension(nirrep_);
109131
nfrzcpi_ = psi::Dimension(nirrep_);
110132
core_mos_ = mo_space_info_->absolute_mo("INACTIVE_DOCC");
133+
label_to_mos_["f"] = std::vector<size_t>();
134+
label_to_mos_["u"] = std::vector<size_t>();
135+
label_to_mos_["v"] = mo_space_info_->absolute_mo("INACTIVE_UOCC");
136+
label_to_cmos_["c"] = mo_space_info_->corr_absolute_mo("INACTIVE_DOCC");
137+
label_to_cmos_["v"] = mo_space_info_->corr_absolute_mo("INACTIVE_UOCC");
138+
} else {
139+
ncmopi_ = mo_space_info_->dimension("CORRELATED");
140+
nfrzvpi_ = mo_space_info_->dimension("FROZEN_UOCC");
141+
nfrzcpi_ = mo_space_info_->dimension("FROZEN_DOCC");
142+
core_mos_ = mo_space_info_->absolute_mo("RESTRICTED_DOCC");
143+
label_to_mos_["f"] = mo_space_info_->absolute_mo("FROZEN_DOCC");
144+
label_to_mos_["u"] = mo_space_info_->absolute_mo("FROZEN_UOCC");
145+
label_to_mos_["v"] = mo_space_info_->absolute_mo("RESTRICTED_UOCC");
146+
label_to_cmos_["c"] = mo_space_info_->corr_absolute_mo("RESTRICTED_DOCC");
147+
label_to_cmos_["v"] = mo_space_info_->corr_absolute_mo("RESTRICTED_UOCC");
111148
}
112149

113-
nactvpi_ = mo_space_info_->dimension("ACTIVE");
114-
actv_mos_ = mo_space_info_->absolute_mo("ACTIVE");
150+
label_to_mos_["c"] = core_mos_;
151+
label_to_mos_["a"] = actv_mos_;
152+
label_to_cmos_["a"] = mo_space_info_->corr_absolute_mo("ACTIVE");
115153

116154
nso_ = nsopi_.sum();
117155
nmo_ = nmopi_.sum();
118156
ncmo_ = ncmopi_.sum();
119157
nactv_ = nactvpi_.sum();
120158
nfrzc_ = nfrzcpi_.sum();
121159

122-
label_to_mos_.clear();
123-
124-
if (freeze_core_) {
125-
label_to_mos_["f"] = mo_space_info_->absolute_mo("FROZEN_DOCC");
126-
} else {
127-
label_to_mos_["f"] = std::vector<size_t>();
128-
}
129-
label_to_mos_["c"] = core_mos_;
130-
label_to_mos_["a"] = actv_mos_;
131-
label_to_mos_["v"] = mo_space_info_->absolute_mo("RESTRICTED_UOCC");
132-
label_to_mos_["u"] = mo_space_info_->absolute_mo("FROZEN_UOCC");
133-
134-
label_to_cmos_.clear();
135-
if (freeze_core_) {
136-
label_to_cmos_["c"] = mo_space_info_->corr_absolute_mo("RESTRICTED_DOCC");
137-
label_to_cmos_["a"] = mo_space_info_->corr_absolute_mo("ACTIVE");
138-
label_to_cmos_["v"] = mo_space_info_->corr_absolute_mo("RESTRICTED_UOCC");
139-
} else {
140-
label_to_cmos_["c"] = mo_space_info_->absolute_mo("INACTIVE_DOCC");
141-
label_to_cmos_["a"] = mo_space_info_->absolute_mo("ACTIVE");
142-
label_to_cmos_["v"] = mo_space_info_->absolute_mo("RESTRICTED_UOCC");
143-
}
144-
145160
// in Pitzer ordering
146161
mos_rel_.resize(nmo_);
147162
for (int h = 0, offset = 0; h < nirrep_; ++h) {
@@ -456,7 +471,11 @@ void MCSCF_ORB_GRAD::build_mo_integrals() {
456471
if (ints_->integral_type() == Custom) {
457472
fill_tei_custom(V_);
458473
} else {
459-
build_tei_from_ao();
474+
if (tei_alg_ == TEIALG::JK) {
475+
build_tei_jk();
476+
} else {
477+
build_tei_df();
478+
}
460479
}
461480
}
462481

@@ -475,35 +494,68 @@ void MCSCF_ORB_GRAD::fill_tei_custom(ambit::BlockedTensor V) {
475494
}
476495
}
477496

478-
void MCSCF_ORB_GRAD::build_tei_from_ao() {
497+
void MCSCF_ORB_GRAD::build_tei_df() {
479498
if (nactv_ == 0)
480499
return;
481500

482-
// This function will do an integral transformation using the JK builder,
483-
// and return the integrals of type <px|uy> = (pu|xy).
484501
timer_on("Build (pu|xy) integrals");
485502

486-
// Transform C matrix to C1 symmetry
487-
// JK does not support mixed symmetry needed for 4-index integrals (York 09/09/2020)
488-
auto aotoso = ints_->wfn()->aotoso();
489-
auto C_nosym = std::make_shared<psi::Matrix>(nso_, nmo_);
490-
491-
// Transform from the SO to the AO basis for the C matrix
492-
// MO in Pitzer ordering and only keep the non-frozen MOs
493-
for (int h = 0, index = 0; h < nirrep_; ++h) {
494-
for (int i = 0; i < nmopi_[h]; ++i) {
495-
int nao = nso_, nso_h = nsopi_[h];
503+
auto C_nosym = ints_->Ca_SO2AO(C_);
496504

497-
if (!nso_h)
498-
continue;
505+
auto Cact = std::make_shared<psi::Matrix>("Cact", nso_, nactv_);
506+
for (size_t x = 0; x < nactv_; ++x) {
507+
Cact->set_column(0, x, C_nosym->get_column(0, actv_mos_[x]));
508+
}
499509

500-
C_DGEMV('N', nao, nso_h, 1.0, aotoso->pointer(h)[0], nso_h, &C_->pointer(h)[0][i],
501-
nmopi_[h], 0.0, &C_nosym->pointer()[0][index], nmo_);
510+
df_helper_->add_space("ALL", C_nosym);
511+
df_helper_->add_space("ACT", Cact);
512+
df_helper_->add_transformation("B", "ALL", "ACT", "Qpq");
513+
df_helper_->transform();
514+
515+
auto B = df_helper_->get_tensor("B"); // naux * (nmo * nact)
516+
517+
size_t naux = df_helper_->get_naux();
518+
auto Baa = std::make_shared<psi::Matrix>("Baa", naux, nactv_ * nactv_);
519+
for (size_t Q = 0; Q < naux; ++Q) {
520+
for (size_t u = 0, nua = 0; u < nactv_; ++u) {
521+
nua = label_to_mos_["a"][u] * nactv_;
522+
Baa->set(Q, u * nactv_ + u, B->get(Q, nua + u));
523+
for (size_t v = u + 1; v < nactv_; ++v) {
524+
Baa->set(Q, u * nactv_ + v, B->get(Q, nua + v));
525+
Baa->set(Q, v * nactv_ + u, B->get(Q, nua + v));
526+
}
527+
}
528+
}
502529

503-
index += 1;
530+
auto puxy = psi::linalg::doublet(B, Baa, true, false); // (nmo * nactv) * (nactv * nactv)
531+
auto puxy_ptr = puxy->get_pointer();
532+
for (const auto& [p_label, p_mos] : label_to_mos_) {
533+
std::string block = p_label + "aaa";
534+
auto& data = V_.block(block).data();
535+
for (size_t p = 0, psize = p_mos.size(), nactv2 = nactv_ * nactv_; p < psize; ++p) {
536+
for (size_t u = 0; u < nactv_; ++u) {
537+
C_DCOPY(nactv2, puxy_ptr + (p_mos[p] * nactv_ + u) * nactv2, 1,
538+
&data[(p * nactv_ + u) * nactv2], 1);
539+
}
504540
}
505541
}
506542

543+
df_helper_->clear_all();
544+
timer_off("Build (pu|xy) integrals");
545+
}
546+
547+
void MCSCF_ORB_GRAD::build_tei_jk() {
548+
if (nactv_ == 0)
549+
return;
550+
551+
// This function will do an integral transformation using the JK builder,
552+
// and return the integrals of type <px|uy> = (pu|xy).
553+
timer_on("Build (pu|xy) integrals");
554+
555+
// Transform C matrix to C1 symmetry
556+
// JK does not support mixed symmetry needed for 4-index integrals (York 09/09/2020)
557+
auto C_nosym = ints_->Ca_SO2AO(C_);
558+
507559
// set up the active part of the C matrix
508560
auto Cact = std::make_shared<psi::Matrix>("Cact", nso_, nactv_);
509561
std::vector<std::shared_ptr<psi::Matrix>> Cact_vec(nactv_);
@@ -693,7 +745,7 @@ std::shared_ptr<psi::Matrix> MCSCF_ORB_GRAD::fock(std::shared_ptr<RDMs> rdms) {
693745
}
694746

695747
double MCSCF_ORB_GRAD::evaluate(std::shared_ptr<psi::Vector> x, std::shared_ptr<psi::Vector> g,
696-
bool do_g) {
748+
bool do_g) {
697749
// if need to update orbitals and integrals
698750
if (update_orbitals(x)) {
699751
build_mo_integrals();
@@ -812,7 +864,7 @@ std::shared_ptr<psi::Matrix> MCSCF_ORB_GRAD::cayley_trans(const std::shared_ptr<
812864

813865
std::vector<std::tuple<int, int, int>>
814866
MCSCF_ORB_GRAD::test_orbital_rotations(const std::shared_ptr<psi::Matrix>& U,
815-
const std::string& warning_msg) {
867+
const std::string& warning_msg) {
816868
// the overlap between new and old orbitals is simply U
817869
// O = Cold^T S Cnew = Cold^T S Cold U = U
818870
// MOM projection index: Pj = sum_{i}^{actv} O_ij
@@ -900,7 +952,7 @@ void MCSCF_ORB_GRAD::compute_orbital_grad() {
900952
}
901953

902954
void MCSCF_ORB_GRAD::hess_diag(std::shared_ptr<psi::Vector>,
903-
const std::shared_ptr<psi::Vector>& h0) {
955+
const std::shared_ptr<psi::Vector>& h0) {
904956
compute_orbital_hess_diag();
905957
h0->copy(*hess_diag_);
906958
}
@@ -1006,7 +1058,7 @@ void MCSCF_ORB_GRAD::compute_orbital_hess_diag() {
10061058
}
10071059

10081060
void MCSCF_ORB_GRAD::reshape_rot_ambit(ambit::BlockedTensor bt,
1009-
const std::shared_ptr<psi::Vector>& sv) {
1061+
const std::shared_ptr<psi::Vector>& sv) {
10101062
size_t vec_size = sv->dimpi().sum();
10111063
if (vec_size != nrot_) {
10121064
throw std::runtime_error(

‎forte/mcscf/mcscf_orb_grad.h

+19-5
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@
3737
#include "psi4/libdiis/diismanager.h"
3838
#include "psi4/libfock/jk.h"
3939
#include "psi4/libmints/matrix.h"
40+
#include "psi4/lib3index/dfhelper.h"
4041

4142
#include "ambit/blocked_tensor.h"
4243

@@ -59,8 +60,8 @@ class MCSCF_ORB_GRAD {
5960
* See J. Chem. Phys. 142, 224103 (2015) and Theor. Chem. Acc. 97, 88-95 (1997)
6061
*/
6162
MCSCF_ORB_GRAD(std::shared_ptr<ForteOptions> options,
62-
std::shared_ptr<MOSpaceInfo> mo_space_info,
63-
std::shared_ptr<ForteIntegrals> ints, bool freeze_core);
63+
std::shared_ptr<MOSpaceInfo> mo_space_info, std::shared_ptr<ForteIntegrals> ints,
64+
bool ignore_frozen);
6465

6566
/// Evaluate the energy and orbital gradient
6667
double evaluate(std::shared_ptr<psi::Vector> x, std::shared_ptr<psi::Vector> g,
@@ -124,6 +125,13 @@ class MCSCF_ORB_GRAD {
124125
/// The JK object of Psi4
125126
std::shared_ptr<psi::JK> JK_;
126127

128+
/// The DFHelper object of Psi4
129+
std::shared_ptr<psi::DFHelper> df_helper_;
130+
131+
/// Algorithm for computing (pu|xy) integrals
132+
enum TEIALG { JK, DF };
133+
TEIALG tei_alg_ = JK;
134+
127135
// => MO spaces related <=
128136

129137
/// The number of irreps
@@ -155,8 +163,8 @@ class MCSCF_ORB_GRAD {
155163
/// The number of frozen-core orbitals
156164
size_t nfrzc_;
157165

158-
/// Freeze core or not
159-
bool freeze_core_ = false;
166+
/// Ignore frozen orbitals in the input
167+
bool ignore_frozen_ = true;
160168

161169
/// List of core MOs (Absolute)
162170
std::vector<size_t> core_mos_;
@@ -257,7 +265,8 @@ class MCSCF_ORB_GRAD {
257265
void build_mo_integrals();
258266

259267
/// Build two-electron integrals
260-
void build_tei_from_ao();
268+
void build_tei_jk();
269+
void build_tei_df();
261270

262271
/// Fill two-electron integrals for custom integrals
263272
void fill_tei_custom(ambit::BlockedTensor V);
@@ -304,6 +313,11 @@ class MCSCF_ORB_GRAD {
304313
void dump_tpdm_iwl();
305314
/// Dump the Hartree-Fock MO 2-RDM to file using IWL
306315
void dump_tpdm_iwl_hf();
316+
/// Dump AO 2-RDM to file for DF-MCSCF
317+
void dump_tpdm_df();
318+
/// Dump the Hartree-Fock AO 2-RDM to file for DF-MCSCF
319+
void dump_tpdm_df_hf(std::shared_ptr<psi::Matrix> Jm12, std::shared_ptr<psi::Matrix> d3,
320+
std::shared_ptr<psi::Matrix> d2);
307321

308322
/// Are there any frozen orbitals?
309323
bool is_frozen_orbs_;

‎forte/mcscf/mcscf_orb_grad_deriv.cc

+372-26
Large diffs are not rendered by default.

‎forte/options.yaml

+9-5
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,10 @@ Driver:
5252
default: "SS"
5353
choices: ["SS", "SA", "MS", "DWMS"]
5454
help: "The type of computation"
55+
MCSCF_REFERENCE:
56+
type: bool
57+
default: True
58+
help: "Whether to use MCSCF reference in Forte or not"
5559
NEL:
5660
type: int
5761
default: None
@@ -1370,14 +1374,14 @@ MCSCF:
13701374
type: int
13711375
default: 1
13721376
help: "How often to solve CI?\n< 1: do CI in the first macro iteration ONLY\n= n: do CI every n macro iteration"
1373-
MCSCF_REFERENCE:
1377+
MCSCF_DF_TEIALG:
13741378
type: bool
13751379
default: True
1376-
help: "Run a FCI followed by MCSCF computation?"
1377-
MCSCF_FREEZE_CORE:
1380+
help: "Algorithm to build (pu|xy) integrals in DF-MCSCF. Use (Q|pu) if True; Otherwise use JK build for every xy pair"
1381+
MCSCF_IGNORE_FROZEN_ORBS:
13781382
type: bool
1379-
default: False
1380-
help: "Freeze core orbitals in MCSCF?"
1383+
default: True
1384+
help: "Ignore frozen orbitals in the input file for MCSCF or not"
13811385
MCSCF_MULTIPLICITY:
13821386
type: int
13831387
default: 0

‎forte/pymodule.py

+9-9
Original file line numberDiff line numberDiff line change
@@ -189,14 +189,14 @@ def energy_forte(name, **kwargs):
189189
psi4.core.print_out("\n\n Skipping MCSCF computation. Using HF or orbitals passed via ref_wfn\n")
190190
else:
191191
active_space_solver_type = data.options.get_str("ACTIVE_SPACE_SOLVER")
192-
mcscf_freeze_core = data.options.get_bool("MCSCF_FREEZE_CORE")
193-
194-
# freeze core orbitals check
195-
frozen_docc_set = data.mo_space_info.size("FROZEN_DOCC") > 0
196-
if not mcscf_freeze_core and frozen_docc_set and data.options.get_str("CORRELATION_SOLVER") == "NONE":
197-
msg = "\n WARNING: By default, Forte will not freeze core orbitals in MCSCF,\n unless the option MCSCF_FREEZE_CORE is set to True.\n"
198-
msg += f"\n Your input file specifies the FROZEN_DOCC array ({data.mo_space_info.size('FROZEN_DOCC')} MOs) in the\n MO_SPACE_INFO block, but the option MCSCF_FREEZE_CORE is set to False.\n"
199-
msg += "\n If you want to freeze the core orbitals in MCSCF, set MCSCF_FREEZE_CORE to True,\n otherwise change the FROZEN_DOCC array to zero(s) and update the RESTRICTED_DOCC array.\n"
192+
mcscf_ignore_frozen = data.options.get_bool("MCSCF_IGNORE_FROZEN_ORBS")
193+
194+
# freeze core/virtual orbitals check
195+
frozen_set = data.mo_space_info.size("FROZEN_DOCC") > 0 or data.mo_space_info.size("FROZEN_UOCC") > 0
196+
if mcscf_ignore_frozen and frozen_set and data.options.get_str("CORRELATION_SOLVER") == "NONE":
197+
msg = "\n WARNING: By default, Forte will not freeze core/virtual orbitals in MCSCF,\n unless the option MCSCF_IGNORE_FROZEN_ORBS is set to False.\n"
198+
msg += f"\n Your input file specifies the FROZEN_DOCC ({data.mo_space_info.size('FROZEN_DOCC')} MOs) / FROZEN_UOCC ({data.mo_space_info.size('FROZEN_UOCC')} MOs) arrays in the\n MO_SPACE_INFO block, but the option MCSCF_IGNORE_FROZEN_ORBS is set to True.\n"
199+
msg += "\n If you want to freeze the core/virtual orbitals in MCSCF, set MCSCF_IGNORE_FROZEN_ORBS to False,\n otherwise change the FROZEN_DOCC/FROZEN_UOCC arrays to zero(s) and update the RESTRICTED_DOCC/RESTRICTED_UOCC arrays.\n"
200200
print(msg)
201201
psi4.core.print_out(msg)
202202

@@ -294,7 +294,7 @@ def gradient_forte(name, **kwargs):
294294
OrbitalTransformation(orb_type, job_type != "NONE").run(data)
295295

296296
active_space_solver_type = data.options.get_str("ACTIVE_SPACE_SOLVER")
297-
mcscf_freeze_core = data.options.get_bool("MCSCF_FREEZE_CORE")
297+
mcscf_ignore_frozen = data.options.get_bool("MCSCF_IGNORE_FROZEN_ORBS")
298298
data = MCSCF(active_space_solver_type).run(data)
299299
energy = data.results.value("energy")
300300

‎tests/manual/mcscf-5/input.dat

+2-2
Original file line numberDiff line numberDiff line change
@@ -26,8 +26,8 @@ set forte {
2626
mcscf_e_convergence 8 # energy convergence of the MCSCF iterations
2727
mcscf_g_convergence 6 # gradient convergence of the MCSCF iterations
2828
mcscf_micro_maxiter 4 # do at most 4 micro iterations per macro iteration
29-
mcscf_freeze_core true # enables freezing the MCSCF core orbitals
29+
mcscf_ignore_frozen_orbs false # enables freezing the MCSCF core orbitals
3030
}
3131

3232
energy('forte')
33-
compare_values(-112.871834862958, variable("CURRENT ENERGY"), 11, "MCSCF energy")
33+
compare_values(-112.871834862958, variable("CURRENT ENERGY"), 8, "MCSCF energy")

‎tests/methods/aci_scf-1/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ set forte{
4141
int_type conventional
4242
sigma 0.0
4343
gamma 0.0
44-
mcscf_freeze_core true
44+
mcscf_ignore_frozen_orbs false
4545
mcscf_do_diis false
4646
}
4747

‎tests/methods/casscf-3/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,7 @@ set forte {
3030
active [3, 0, 1, 2]
3131
mcscf_e_convergence 1e-10
3232
mcscf_g_convergence 1e-6
33-
mcscf_freeze_core true
33+
mcscf_ignore_frozen_orbs false
3434
}
3535

3636
mcscf_forte = energy('forte')

‎tests/methods/casscf-5/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,7 @@ set forte {
3131
root_sym 0
3232
nroot 1
3333
multiplicity 1
34-
mcscf_freeze_core true
34+
mcscf_ignore_frozen_orbs false
3535
}
3636

3737
energy('forte')

‎tests/methods/casscf-6/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ set forte{
2020
frozen_docc [1,0,0,0]
2121
restricted_docc [1,0,1,1]
2222
active [2,0,0,0]
23-
mcscf_freeze_core true
23+
mcscf_ignore_frozen_orbs false
2424
}
2525

2626
energy('forte')

‎tests/methods/casscf-fcidump-1/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ set forte {
2424
frozen_uocc [0,0,0,0]
2525
restricted_docc [2,0,0,0]
2626
active [2,0,2,2]
27-
mcscf_freeze_core true
27+
mcscf_ignore_frozen_orbs false
2828
e_convergence 8
2929
r_convergence 8
3030
}

‎tests/methods/casscf-gradient-2/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@ set forte{
3333
mcscf_g_convergence 1e-12
3434
mcscf_e_convergence 1e-12
3535
cpscf_convergence 1e-10
36-
mcscf_freeze_core true
36+
mcscf_ignore_frozen_orbs false
3737
}
3838

3939
grad = gradient('forte')

‎tests/methods/casscf-gradient-3/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,7 @@ set forte{
4848
mcscf_maxiter 100
4949
mcscf_g_convergence 1e-10
5050
mcscf_e_convergence 1e-12
51-
mcscf_freeze_core true
51+
mcscf_ignore_frozen_orbs false
5252
cpscf_convergence 1e-10
5353
}
5454

‎tests/methods/casscf-opt-3/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ set globals{
2222
set forte{
2323
active_space_solver fci
2424
frozen_docc [1,0,0,0]
25-
mcscf_freeze_core true
25+
mcscf_ignore_frozen_orbs false
2626
}
2727

2828
Eopt = optimize('forte')

‎tests/methods/df-casscf-1/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@ set forte{
4040
mcscf_maxiter 25
4141
mcscf_e_convergence 7
4242
mcscf_g_convergence 7
43-
mcscf_freeze_core true
43+
mcscf_ignore_frozen_orbs false
4444
print 0
4545
}
4646
Ecas_1 = energy('forte')

‎tests/methods/df-casscf-gradient-1/input.dat

+54
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
# A test of analytic DF-CASSCF gradients on N2
2+
3+
import forte
4+
5+
ref_grad = psi4.Matrix.from_list([
6+
[ 0.0000000000, 0.0000000000, 0.4930152443],
7+
[-0.0000000000, 0.0000000000, -0.4930152443]
8+
])
9+
ref_grad_fc = psi4.Matrix.from_list([
10+
[ 0.0000000000, -0.0000000000, 0.4931039074],
11+
[-0.0000000000, 0.0000000000, -0.4931039074]
12+
])
13+
14+
molecule N2{
15+
N
16+
N 1 1.0
17+
}
18+
19+
set globals {
20+
scf_type df
21+
reference rhf
22+
e_convergence 10
23+
d_convergence 8
24+
maxiter 100
25+
basis cc-pvdz
26+
df_basis_mp2 cc-pvdz-jkfit
27+
df_basis_scf cc-pvdz-jkfit
28+
docc [3,0,0,0,0,2,1,1]
29+
}
30+
31+
set forte{
32+
int_type df
33+
active_space_solver fci
34+
restricted_docc [2,0,0,0,0,2,0,0]
35+
active [1,0,1,1,0,1,1,1]
36+
e_convergence 12
37+
mcscf_maxiter 100
38+
mcscf_g_convergence 1e-12
39+
mcscf_e_convergence 1e-12
40+
cpscf_convergence 1e-12
41+
mcscf_ignore_frozen_orbs false
42+
}
43+
44+
#set gradient_write true
45+
#set findif points 5
46+
grad = gradient('forte')
47+
compare_matrices(ref_grad, grad, 8, "DF-CASSCF(6,6)/cc-pVDZ/AE gradient on N2")
48+
49+
set forte{
50+
frozen_docc [1,0,0,0,0,1,0,0]
51+
restricted_docc [1,0,0,0,0,1,0,0]
52+
}
53+
grad = gradient('forte')
54+
compare_matrices(ref_grad_fc, grad, 8, "DF-CASSCF(6,6)/cc-pVDZ/FC gradient on N2")

‎tests/methods/df-casscf-gradient-1/output.ref

+1,722
Large diffs are not rendered by default.

‎tests/methods/df-casscf-gradient-2/input.dat

+62
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
# Test DF-MCSCF analytic gradients with frozen orbitals
2+
3+
import forte
4+
5+
# ref_grad from 5-point finite difference
6+
ref_grad = psi4.Matrix.from_list([
7+
[-0.0009775460, -0.0082154776, 0.0000000000],
8+
[ 0.0009775460, 0.0082154776, 0.0000000000],
9+
[ 0.0068717247, 0.0043753883, 0.0000000000],
10+
[-0.0068717247, -0.0043753883, 0.0000000000],
11+
[-0.0085117226, -0.0044009360, 0.0000000000],
12+
[ 0.0085117226, 0.0044009360, 0.0000000000],
13+
[-0.0186375215, -0.0052285808, 0.0000000000],
14+
[ 0.0186375215, 0.0052285808, 0.0000000000],
15+
[-0.0190885742, 0.0141712291, 0.0000000000],
16+
[ 0.0190885742, -0.0141712291, 0.0000000000]
17+
])
18+
19+
molecule butadiene{
20+
H 1.080977 -2.558832 0.000000
21+
H -1.080977 2.558832 0.000000
22+
H 2.103773 -1.017723 0.000000
23+
H -2.103773 1.017723 0.000000
24+
H -0.973565 -1.219040 0.000000
25+
H 0.973565 1.219040 0.000000
26+
C 0.000000 0.728881 0.000000
27+
C 0.000000 -0.728881 0.000000
28+
C 1.117962 -1.474815 0.000000
29+
C -1.117962 1.474815 0.000000
30+
}
31+
32+
set globals {
33+
reference rhf
34+
scf_type df
35+
e_convergence 12
36+
d_convergence 10
37+
maxiter 100
38+
basis dz
39+
docc [7,1,1,6]
40+
df_basis_scf def2-universal-jkfit
41+
df_basis_mp2 def2-universal-jkfit
42+
}
43+
44+
set forte{
45+
int_type df
46+
active_space_solver fci
47+
frozen_docc [2,0,0,2]
48+
restricted_docc [5,0,0,4]
49+
active [0,2,2,0]
50+
frozen_uocc [2,0,0,2]
51+
e_convergence 1e-12
52+
mcscf_maxiter 100
53+
mcscf_g_convergence 1e-10
54+
mcscf_e_convergence 1e-12
55+
cpscf_convergence 1e-10
56+
mcscf_ignore_frozen_orbs false
57+
}
58+
59+
#set gradient_write true
60+
#set findif points 5
61+
grad = gradient('forte')
62+
compare_matrices(ref_grad, grad, 7, "CASSCF(4,4)/DZ gradient on butadiene with frozen orbitals")

‎tests/methods/df-casscf-gradient-2/output.ref

+960
Large diffs are not rendered by default.

‎tests/methods/df-dsrg-mrpt2-7-localized-actv/input.dat

+1-1
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ set forte{
4949
mcscf_g_convergence 1.0e-6
5050
mcscf_e_convergence 1.0e-12
5151
ci_spin_adapt true
52-
mcscf_freeze_core true
52+
mcscf_ignore_frozen_orbs false
5353
}
5454

5555
Ecas, wfn = energy('forte', ref_wfn=wfn, return_wfn=True)

‎tests/methods/tests.yaml

+3-1
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,7 @@ casscf:
7777
- casscf-gradient-1
7878
- casscf-gradient-3
7979
- casscf-opt-3
80+
- df-casscf-gradient-1
8081
medium:
8182
- casscf-1
8283
- casscf-2
@@ -86,6 +87,7 @@ casscf:
8687
- casscf-6
8788
- casscf-9
8889
- casscf-gradient-2
90+
- df-casscf-gradient-2
8991
- casscf-opt-1
9092
# - casscf-opt-2
9193
- sa-casscf-1
@@ -430,4 +432,4 @@ x2c:
430432
- x2c-1
431433
pyscf:
432434
short:
433-
- pyscf_interface
435+
- pyscf_interface

0 commit comments

Comments
 (0)
Please sign in to comment.