search

GQCG / GQCP

The Ghent Quantum Chemistry Package for electronic structure calculations
https://gqcg.github.io/GQCP/
GNU Lesser General Public License v3.0
34 stars 10 forks source link

Rewrite frozen core (spin-resolved) CI DMs #741

Closed lelemmen closed 3 years ago

lelemmen commented 3 years ago

Describe a change that would improve the code base We should make sure that this implementation is as clear as possible.

lelemmen commented 3 years ago

BaseSpinResolvedFrozenCoreCIDMCalculator.hpp

// This file is part of GQCG-GQCP.
//
// Copyright (C) 2017-2020  the GQCG developers
//
// GQCG-GQCP is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// GQCG-GQCP is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with GQCG-GQCP.  If not, see <http://www.gnu.org/licenses/>.

#include "DensityMatrix/CIDMCalculators/BaseSpinResolvedFrozenDMCalculator.hpp"

#include "Utilities/linalg.hpp"

namespace GQCP {

/*
 *  CONSTRUCTORS
 */

/**
 *  @param dm_calculator                shared pointer to active (non-frozen core) DM builder
 *  @param X                            the number of frozen orbitals
 */
BaseSpinResolvedFrozenDMCalculator::BaseSpinResolvedFrozenDMCalculator(const std::shared_ptr<BaseSpinResolvedDMCalculator> dm_calculator, const size_t X) :
    BaseSpinResolvedDMCalculator(),
    active_dm_calculator {std::move(dm_calculator)},
    X {X} {}

/*
 * PUBLIC OVERRIDDEN METHODS
 */

/**
 *  @param x        the coefficient vector representing the wave function
 *
 *  @return all 1-DMs given a coefficient vector
 */
SpinResolved1DM<double> BaseSpinResolvedFrozenDMCalculator::calculateSpinResolved1DM(const VectorX<double>& x) const {

    auto K = this->onvBasis()->numberOfOrbitals();

    SpinResolved1DMComponent<double> D_aa = SpinResolved1DMComponent<double>::Zero(K);
    SpinResolved1DMComponent<double> D_bb = SpinResolved1DMComponent<double>::Zero(K);

    auto K_active = K - this->X;

    // Set the values for the frozen orbital
    for (size_t i = 0; i < this->X; i++) {
        D_aa(i, i) = 1;
        D_bb(i, i) = 1;
    }

    SpinResolved1DM<double> sub_1DMs = this->active_dm_calculator->calculateSpinResolved1DM(x);

    // Incorporate the submatrices from the active space
    D_aa.block(this->X, this->X, K_active, K_active) += sub_1DMs.alpha();
    D_bb.block(this->X, this->X, K_active, K_active) += sub_1DMs.beta();

    return SpinResolved1DM<double>(D_aa, D_bb);
};

/**
 *  @param x        the coefficient vector representing the wave function
 *
 *  @return all 2-DMs given a coefficient vector
 */
SpinResolved2DM<double> BaseSpinResolvedFrozenDMCalculator::calculateSpinResolved2DM(const VectorX<double>& x) const {

    auto K = this->onvBasis()->numberOfOrbitals();

    SpinResolved2DMComponent<double> d_aaaa = SpinResolved2DMComponent<double>::Zero(K);
    SpinResolved2DMComponent<double> d_aabb = SpinResolved2DMComponent<double>::Zero(K);
    SpinResolved2DMComponent<double> d_bbaa = SpinResolved2DMComponent<double>::Zero(K);
    SpinResolved2DMComponent<double> d_bbbb = SpinResolved2DMComponent<double>::Zero(K);

    SpinResolved1DM<double> one_DMs = this->active_dm_calculator->calculateSpinResolved1DM(x);
    auto D_aa = one_DMs.alpha();
    auto D_bb = one_DMs.beta();

    // Implement frozen DM formulas
    for (size_t p = 0; p < this->X; p++) {  // iterate over frozen orbitals

        // DM Overlap between frozen and active space:
        //      frozen orbital indices (p) must always have one annihilation and one creation index (always occupied)
        //      values are dictated by the 'active' orbital indices and correspond to that of the active 1DMs
        //      Hence we start adding the 1DMs starting from index 'X' the number frozen orbitals
        d_aaaa.addBlock<0, 1>(D_aa, this->X, this->X, p, p);
        d_aaaa.addBlock<2, 3>(D_aa, p, p, this->X, this->X);
        d_aaaa.addBlock<2, 1>(-D_aa, p, this->X, this->X, p);
        d_aaaa.addBlock<3, 0>(-D_aa, this->X, p, p, this->X);

        d_bbbb.addBlock<0, 1>(D_bb, this->X, this->X, p, p);
        d_bbbb.addBlock<2, 3>(D_bb, p, p, this->X, this->X);
        d_bbbb.addBlock<2, 1>(-D_bb, p, this->X, this->X, p);
        d_bbbb.addBlock<3, 0>(-D_bb, this->X, p, p, this->X);

        d_aabb.addBlock<2, 3>(D_bb, p, p, this->X, this->X);
        d_aabb.addBlock<0, 1>(D_aa, this->X, this->X, p, p);

        d_bbaa.addBlock<2, 3>(D_aa, p, p, this->X, this->X);
        d_bbaa.addBlock<0, 1>(D_bb, this->X, this->X, p, p);

        // Set the values for the frozen orbital
        d_bbaa(p, p, p, p) = 1;
        d_aabb(p, p, p, p) = 1;

        for (size_t q = p + 1; q < this->X; q++) {  // iterate over frozen orbitals

            d_aaaa(p, p, q, q) = 1;
            d_aaaa(q, q, p, p) = 1;
            d_aaaa(p, q, q, p) = -1;
            d_aaaa(q, p, p, q) = -1;

            d_bbbb(p, p, q, q) = 1;
            d_bbbb(q, q, p, p) = 1;
            d_bbbb(p, q, q, p) = -1;
            d_bbbb(q, p, p, q) = -1;

            d_aabb(p, p, q, q) = 1;
            d_bbaa(p, p, q, q) = 1;

            d_aabb(q, q, p, p) = 1;
            d_bbaa(q, q, p, p) = 1;
        }
    }

    // Incorporate the 2-DM subblocks into the total 2DMs
    SpinResolved2DM<double> sub_2DMs = this->active_dm_calculator->calculateSpinResolved2DM(x);

    d_aaaa.addBlock(sub_2DMs.alphaAlpha(), this->X, this->X, this->X, this->X);
    d_bbbb.addBlock(sub_2DMs.betaBeta(), this->X, this->X, this->X, this->X);
    d_aabb.addBlock(sub_2DMs.alphaBeta(), this->X, this->X, this->X, this->X);
    d_bbaa.addBlock(sub_2DMs.betaAlpha(), this->X, this->X, this->X, this->X);

    return SpinResolved2DM<double>(d_aaaa, d_aabb, d_bbaa, d_bbbb);
};

}  // namespace GQCP

SpinResolvedFrozenDMCalculator.cpp

// This file is part of GQCG-GQCP.
//
// Copyright (C) 2017-2020  the GQCG developers
//
// GQCG-GQCP is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// GQCG-GQCP is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with GQCG-GQCP.  If not, see <http://www.gnu.org/licenses/>.

#include "DensityMatrix/CIDMCalculators/SpinResolvedFrozenDMCalculator.hpp"

#include "DensityMatrix/CIDMCalculators/SpinResolvedDMCalculator.hpp"

namespace GQCP {

/*
 *  CONSTRUCTORS
 */

/**
 *  @param onv_basis       the frozen spin-resolved ONV basis
 */
SpinResolvedFrozenDMCalculator::SpinResolvedFrozenDMCalculator(const SpinResolvedFrozenONVBasis& onv_basis) :
    BaseSpinResolvedFrozenDMCalculator(std::make_shared<SpinResolvedDMCalculator>(onv_basis.activeONVBasis()), onv_basis.numberOfFrozenOrbitals()),
    onv_basis {onv_basis} {}

}  // namespace GQCP
lelemmen commented 3 years ago

Test case

// // This file is part of GQCG-GQCP.
// //
// // Copyright (C) 2017-2020  the GQCG developers
// //
// // GQCG-GQCP is free software: you can redistribute it and/or modify
// // it under the terms of the GNU Lesser General Public License as published by
// // the Free Software Foundation, either version 3 of the License, or
// // (at your option) any later version.
// //
// // GQCG-GQCP is distributed in the hope that it will be useful,
// // but WITHOUT ANY WARRANTY; without even the implied warranty of
// // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// // GNU Lesser General Public License for more details.
// //
// // You should have received a copy of the GNU Lesser General Public License
// // along with GQCG-GQCP.  If not, see <http://www.gnu.org/licenses/>.

// #define BOOST_TEST_MODULE "SpinResolvedFrozenDMCalculator_test"

// #include <boost/test/unit_test.hpp>

// #include "DensityMatrix/CIDMCalculators/SpinResolvedFrozenDMCalculator.hpp"
// #include "DensityMatrix/CIDMCalculators/SpinResolvedSelectedDMCalculator.hpp"
// #include "Mathematical/Optimization/Eigenproblem/EigenproblemSolver.hpp"
// #include "Operator/SecondQuantized/SQHamiltonian.hpp"
// #include "QCMethod/CI/CI.hpp"
// #include "QCMethod/CI/CIEnvironment.hpp"

// /**
//  *  Check if the 1- and 2-DMs for a frozen core spin-resolved ONV basis are equal to the 'selected' case.
//  *  The system of interest is a linear chain of 5 H atoms, 1.1 bohr apart, using an STO-3G basisset.
//  */
// BOOST_AUTO_TEST_CASE(FrozenCoreFCI_one_DMs) {

//     // Set up the molecular Hamiltonian for H5//STO-3G in the Löwdin basis.
//     const GQCP::Molecule molecule = GQCP::Molecule::HChain(5, 1.1);
//     const auto N_alpha = 3;
//     const auto N_beta = 2;

//     GQCP::RSpinOrbitalBasis<double, GQCP::GTOShell> spinor_basis {molecule, "STO-3G"};
//     const auto K = spinor_basis.numberOfSpatialOrbitals();
//     spinor_basis.lowdinOrthonormalize();

//     const auto sq_hamiltonian = GQCP::RSQHamiltonian<double>::Molecular(spinor_basis, molecule);

//     // Do a dense frozen core CI calculation for 2 frozen orbitals.
//     const GQCP::SpinResolvedFrozenONVBasis onv_basis {K, N_alpha, N_beta, 2};

//     auto environment = GQCP::CIEnvironment::Dense(sq_hamiltonian, onv_basis);
//     auto solver = GQCP::EigenproblemSolver::Dense();

//     const auto linear_expansion = GQCP::QCMethod::CI<GQCP::SpinResolvedFrozenONVBasis>(onv_basis).optimize(solver, environment).groundStateParameters();

//     // Calculate the 1-DMs using specialized spin-resolved and 'selected' routines, and check if they are equal.
//     const GQCP::SpinResolvedFrozenDMCalculator spin_resolved_dm_calculator {onv_basis};
//     const auto one_DMs_specialized = spin_resolved_dm_calculator.calculateSpinResolved1DM(linear_expansion.coefficients());

//     const GQCP::SpinResolvedSelectedONVBasis selected_onv_basis {onv_basis};
//     const GQCP::SpinResolvedSelectedDMCalculator selected_dm_calculator {selected_onv_basis};
//     const auto one_DMs_selected = selected_dm_calculator.calculateSpinResolved1DM(linear_expansion.coefficients());

//     BOOST_CHECK(one_DMs_specialized.orbitalDensity().isApprox(one_DMs_selected.orbitalDensity(), 1.0e-12));
//     BOOST_CHECK(one_DMs_specialized.alpha().isApprox(one_DMs_selected.alpha(), 1.0e-12));
//     BOOST_CHECK(one_DMs_specialized.beta().isApprox(one_DMs_selected.beta(), 1.0e-12));

//     // Calculate the 2-DMs using specialized spin-resolved and 'selected' routines, and check if they are equal.
//     const auto two_DMs_specialized = spin_resolved_dm_calculator.calculateSpinResolved2DM(linear_expansion.coefficients());
//     const auto two_DMs_selected = selected_dm_calculator.calculateSpinResolved2DM(linear_expansion.coefficients());

//     BOOST_CHECK(two_DMs_specialized.alphaAlpha().isApprox(two_DMs_selected.alphaAlpha(), 1.0e-12));
//     BOOST_CHECK(two_DMs_specialized.alphaBeta().isApprox(two_DMs_selected.alphaBeta(), 1.0e-12));
//     BOOST_CHECK(two_DMs_specialized.betaAlpha().isApprox(two_DMs_selected.betaAlpha(), 1.0e-12));
//     BOOST_CHECK(two_DMs_specialized.betaBeta().isApprox(two_DMs_selected.betaBeta(), 1.0e-12));
//     BOOST_CHECK(two_DMs_specialized.orbitalDensity().isApprox(two_DMs_selected.orbitalDensity(), 1.0e-12));
// }