Frozen core CI
#848
Replies: 5 comments
-
Frozen core DOCI density matricesTest 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 "FrozenCoreDOCIRDM_test"
#include <boost/test/unit_test.hpp>
#include "DensityMatrix/CIDMCalculators/FrozenCoreDOCIRDMBuilder.hpp"
#include "DensityMatrix/CIDMCalculators/SeniorityZeroDMCalculator.hpp"
#include "DensityMatrix/CIDMCalculators/SpinResolvedSelectedDMCalculator.hpp"
#include "Operator/SecondQuantized/SQHamiltonian.hpp"
#include "QCMethod/CI/HamiltonianBuilder/FrozenCoreDOCI.hpp"
#include "Utilities/linalg.hpp"
BOOST_AUTO_TEST_CASE(FrozenCoreDOCI_one_DMs) {
size_t K = 5;
GQCP::Molecule H5 = GQCP::Molecule::HChain(K, 1.1);
GQCP::RSpinOrbitalBasis<double, GQCP::GTOShell> spinor_basis {H5, "STO-3G"};
auto sq_hamiltonian = GQCP::RSQHamiltonian<double>::Molecular(spinor_basis, H5); // in an AO basis
GQCP::SpinUnresolvedFrozenONVBasis fock_space {K, 3, 2};
GQCP::SpinResolvedSelectedONVBasis selected_fock_space {fock_space};
GQCP::FrozenCoreDOCI doci {fock_space};
// Specify solver options and solve the eigenvalue problem
// Solve the dense DOCI eigenvalue problem
GQCP::CISolver ci_solver {doci, sq_hamiltonian};
GQCP::DenseSolverOptions solver_options;
ci_solver.solve(solver_options);
GQCP::VectorX<double> coef = ci_solver.eigenpair().eigenvector();
// Get the frozen core DOCI and SelectedCI 1-DMs
GQCP::SpinResolvedSelectedDMCalculator sci_rdm {selected_fock_space};
GQCP::FrozenCoreDOCIRDMBuilder doci_rdm {fock_space};
GQCP::SpinResolved1DM<double> one_DMs_s = sci_rdm.calculateSpinResolved1DM(coef);
GQCP::SpinResolved1DM<double> one_DMs = doci_rdm.calculateSpinResolved1DM(coef);
BOOST_CHECK(one_DMs_s.orbitalDensity().isApprox(one_DMs.orbitalDensity()));
BOOST_CHECK(one_DMs_s.alpha().isApprox(one_DMs.alpha()));
BOOST_CHECK(one_DMs_sbeta().isApprox(one_DMs.beta()));
}
BOOST_AUTO_TEST_CASE(FrozenCoreDOCI_two_DMs) {
size_t K = 5;
GQCP::Molecule H5 = GQCP::Molecule::HChain(K, 1.1);
GQCP::RSpinOrbitalBasis<double, GQCP::GTOShell> spinor_basis {H5, "STO-3G"};
auto sq_hamiltonian = GQCP::RSQHamiltonian<double>::Molecular(spinor_basis, H5); // in an AO basis
GQCP::SpinUnresolvedFrozenONVBasis fock_space {K, 3, 2};
GQCP::SpinResolvedSelectedONVBasis selected_fock_space {fock_space};
GQCP::FrozenCoreDOCI doci {fock_space};
// Specify solver options and solve the eigenvalue problem
// Solve the dense DOCI eigenvalue problem
GQCP::CISolver ci_solver {doci, sq_hamiltonian};
GQCP::DenseSolverOptions solver_options;
ci_solver.solve(solver_options);
GQCP::VectorX<double> coef = ci_solver.eigenpair().eigenvector();
// Get the frozen core DOCI and SelectedCI 2-DMs
GQCP::SpinResolvedSelectedDMCalculator sci_rdm {selected_fock_space};
GQCP::FrozenCoreDOCIRDMBuilder doci_rdm {fock_space};
GQCP::SpinResolved2DM<double> two_DMs_s = sci_rdm.calculateSpinResolved2DM(coef);
GQCP::SpinResolved2DM<double> two_DMs = doci_rdm.calculateSpinResolved2DM(coef);
BOOST_CHECK(two_DMs_s.alphaAlpha().isApprox(two_DMs.alphaAlpha(), 1.0e-06));
BOOST_CHECK(two_DMs_s.alphaBeta().isApprox(two_DMs.alphaBeta(), 1.0e-06));
BOOST_CHECK(two_DMs_s.betaAlpha().isApprox(two_DMs.betaAlpha(), 1.0e-06));
BOOST_CHECK(two_DMs_s.betaBeta().isApprox(two_DMs.betaBeta(), 1.0e-06));
BOOST_CHECK(two_DMs_s.orbitalDensity().isApprox(two_DMs.orbitalDensity(), 1.0e-06));
}
|
Beta Was this translation helpful? Give feedback.
0 replies
-
Frozen core DOCI
|
Beta Was this translation helpful? Give feedback.
0 replies
-
Frozen core 'generalized' evaluationsHere's the code that we used to support: Header file// 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/>.
#pragma once
#include "ONVBasis/BaseFrozenCoreONVBasis.hpp"
#include "ONVBasis/BaseONVBasis.hpp"
#include "ONVBasis/ONVManipulator.hpp"
#include "ONVBasis/SpinUnresolvedONVBasis.hpp"
namespace GQCP {
/**
* A spin-unresolved frozen ONV basis: this is a subset of the N-electron spin-unresolved ONV basis in which the first X orbitals are always occupied.
*/
class SpinUnresolvedFrozenONVBasis: public BaseFrozenCoreONVBasis, public ONVManipulator<SpinUnresolvedFrozenONVBasis> {
protected:
size_t X; // number of frozen orbitals/electrons
SpinUnresolvedONVBasis active_onv_basis; // active (non-frozen) spin-unresolved ONV basis containing only the active electrons (N-X) and orbitals (K-X)
public:
// CONSTRUCTORS
/**
* @param K the total number of orbitals
* @param N the total number of electrons
* @param X the number of frozen orbitals and electrons
*/
SpinUnresolvedFrozenONVBasis(const size_t K, const size_t N, const size_t X);
/**
* @param onv_basis (to be frozen) full spin-unresolved ONV basis
* @param X the number of frozen orbitals and electrons
*/
SpinUnresolvedFrozenONVBasis(const SpinUnresolvedONVBasis& onv_basis, const size_t X);
// DESTRUCTOR
/**
* The default destructor.
*/
~SpinUnresolvedFrozenONVBasis() override = default;
// OVERRIDEN PUBLIC METHODS
/**
* @param representation a representation of a spin-unresolved ONV
*
* @return the address (i.e. the ordering number) of the given spin-unresolved ONV
*/
size_t addressOf(const size_t representation) const override;
/**
* @param onv the spin-unresolved ONV
*
* @return the number of ONVs (with a larger address) this spin-unresolved ONV would couple with given a one electron operator
*/
size_t countOneElectronCouplings(const SpinUnresolvedONV& onv) const override;
/**
* @return the number of non-zero (non-diagonal) couplings of a one electron coupling scheme in the spin-unresolved ONV basis
*/
size_t countTotalOneElectronCouplings() const override { return this->active_onv_basis.countTotalOneElectronCouplings(); }
/**
* @return the number of non-zero (non-diagonal) couplings of a two electron coupling scheme in the spin-unresolved ONV basis
*/
size_t countTotalTwoElectronCouplings() const override { return this->active_onv_basis.countTotalTwoElectronCouplings(); }
/**
* @param onv the spin-unresolved ONV
*
* @return the number of ONVs (with a larger address) this spin-unresolved ONV would couple with given a two electron operator
*/
size_t countTwoElectronCouplings(const SpinUnresolvedONV& onv) const override;
/**
* @param representation a representation of a spin-unresolved ONV
*
* @return the next bitstring permutation in the frozen spin-unresolved ONV basis
*
* Examples (first orbital is frozen):
* 0101 -> 1001
* 01101 -> 10011
*/
size_t nextPermutationOf(const size_t representation) const override;
/**
* Calculate unsigned representation for a given address
*
* @param address the address of the representation is calculated
*
* @return unsigned representation of the address
*/
size_t representationOf(const size_t address) const override;
// PUBLIC METHODS
/**
* @return this' ONV basis that is considered active
*/
const SpinUnresolvedONVBasis& activeONVBasis() const { return this->active_onv_basis; }
using ONVManipulator<SpinUnresolvedFrozenONVBasis>::addressOf;
/**
* @return the number of orbitals that are considered frozen in this ONV basis
*/
size_t numberOfFrozenOrbitals() const { return this->X; }
};
} // namespace GQCPSource file// 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 "ONVBasis/SpinUnresolvedFrozenONVBasis.hpp"
namespace GQCP {
/*
* CONSTRUCTORS
*/
/**
* @param M the total number of orbitals
* @param N the total number of electrons
* @param X the number of frozen orbitals and electrons
*/
SpinUnresolvedFrozenONVBasis::SpinUnresolvedFrozenONVBasis(const size_t M, const size_t N, const size_t X) :
BaseFrozenCoreONVBasis(std::make_shared<SpinUnresolvedONVBasis>(SpinUnresolvedONVBasis(M - X, N - X)), X),
ONVManipulator(N),
active_onv_basis {M - X, N - X},
X {X} {}
/**
* @param fock_space (to be frozen) full spin-unresolved ONV basis
* @param X the number of frozen orbitals and electrons
*/
SpinUnresolvedFrozenONVBasis::SpinUnresolvedFrozenONVBasis(const SpinUnresolvedONVBasis& fock_space, const size_t X) :
SpinUnresolvedFrozenONVBasis(fock_space.numberOfOrbitals(), fock_space.numberOfElectrons(), X) {}
/*
* OVERRIDDEN PUBLIC METHODS
*/
/**
* @param representation a representation of an ONV
*
* @return the address (i.e. the ordering number) of the given ONV
*/
size_t SpinUnresolvedFrozenONVBasis::addressOf(const size_t representation) const {
// Transform the representation to the sub space, by bitshifting left X number of times to remove the frozen orbital indices
// Address of the total ONV in this frozen ONV basis is identical to that of the sub ONV in the sub ONV basis.
return this->active_onv_basis.addressOf(representation >> this->X);
};
/**
* @param onv the ONV
*
* @return the number of ONVs (with a larger address) the given ONV would couple with given a one electron operator
*/
size_t SpinUnresolvedFrozenONVBasis::countOneElectronCouplings(const SpinUnresolvedONV& onv) const {
size_t V = this->M - N; // number of virtual orbitals
size_t coupling_count = 0;
for (size_t e1 = this->X; e1 < this->N; e1++) { // start from X to ignore the frozen electrons
size_t p = onv.occupationIndexOf(e1);
coupling_count += (V + e1 - p); // number of virtuals with an index larger than p
}
return coupling_count;
}
/**
* @param onv the ONV
*
* @return the number of ONVs (with a larger address) the given ONV would couple with given a two electron operator
*/
size_t SpinUnresolvedFrozenONVBasis::countTwoElectronCouplings(const SpinUnresolvedONV& onv) const {
size_t V = this->M - N; // number of virtual orbitals
size_t coupling_count = 0;
for (size_t e1 = this->X; e1 < this->N; e1++) { // start from X to ignore the frozen electrons
size_t p = onv.occupationIndexOf(e1);
coupling_count += (V + e1 - p); // one electron part
for (size_t e2 = e1 + 1; e2 < this->N; e2++) {
size_t q = onv.occupationIndexOf(e2);
size_t coupling_count2 = (V + e2 - q);
coupling_count += (V - coupling_count2) * coupling_count2;
if (coupling_count2 > 1) {
coupling_count += SpinUnresolvedONVBasis::calculateDimension(coupling_count2, 2);
}
}
}
return coupling_count;
}
/**
* @param representation a representation of an ONV
*
* @return the next bitstring permutation in this ONV basis
*
* Examples (first orbital is frozen):
* 0101 -> 1001
* 01101 -> 10011
*/
size_t SpinUnresolvedFrozenONVBasis::nextPermutationOf(size_t representation) const {
// Generate the permutation from the active space, bitshift left X number of times to remove the frozen orbital indices before passing it to the active space
size_t sub_permutation = this->active_onv_basis.nextPermutationOf(representation >> this->X);
// Transform the permutation to the frozen core space, by bitshifting right X number of times and filling the new 0 bits with 1's (by adding 2^X - 1)
sub_permutation <<= X;
sub_permutation += static_cast<size_t>(pow(2, X) - 1);
return sub_permutation;
};
/**
* Calculate unsigned representation for a given address
*
* @param address the address of the representation is calculated
*
* @return unsigned representation of the address
*/
size_t SpinUnresolvedFrozenONVBasis::representationOf(size_t address) const {
// generate the representation in the active space
size_t representation = this->active_onv_basis.representationOf(address);
// transform the permutation to the frozen core space, by bitshifting right X number of times and filling the new 0 bits with 1's (by adding 2^X - 1)
representation <<= X;
representation += static_cast<size_t>(pow(2, X) - 1);
return representation;
};
} // namespace GQCP
|
Beta Was this translation helpful? Give feedback.
0 replies
-
Frozen core 'restricted' evaluations
|
Beta Was this translation helpful? Give feedback.
0 replies
-
|
Beta Was this translation helpful? Give feedback.
0 replies
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
-
In #688, we made a judgment to drop support for frozen core CI calculations. The reason why is that, at the time of its first implementation, there was no difference between a seniority-zero ONV basis and a (spin-unresolved) FCI ONV basis. The formulas we had implemented did not distinguish between these two cases, so they should be re-checked with caution.
Frozen core evaluations
Frozen core density matrices
Beta Was this translation helpful? Give feedback.
All reactions