Permalink
Browse files

Experimental hard-coded scf routine

might not work properly
  • Loading branch information...
mfherbst committed Jun 25, 2016
1 parent 2147793 commit e264896d28c651ab6035d616385c4f728011c6ee
@@ -9,3 +9,9 @@

# our build directory
build

# gtags files:
GPATH
GRTAGS
GTAGS

@@ -0,0 +1,35 @@
#pragma once
#include <linalgwrap/ArmadilloMatrix.hh>
namespace hard_coded_hf {
using namespace linalgwrap;

ArmadilloMatrix<double> hack_guess(const ArmadilloMatrix<double>& s_bb) {
// apply Löwdin normalisation to the basis functions
// - Diagonalise the overlap
// - Take 1/\sqrt{evals} at the diagonal
// - results in orthonormalised basis functions

const arma::mat& s_bb_data = s_bb.data();

// Diagonalize s_bb
arma::vec eval;
arma::mat evec;
arma::eig_sym(eval, evec, s_bb_data);

// take 1/sqrt( . ) for each eigenvalue:
std::transform(std::begin(eval), std::end(eval), std::begin(eval),
[](double elem) { return 1. / std::sqrt(elem); });

// construct new Armadillo matrix of guess vectors:
arma::mat guess = evec * arma::diagmat(eval);
// If number of orbitals is different from number of basis functions,
// guess should be rectangular.
// int b = s_bb.n_rows;
// guess = slice(evec * diagmat(eval), b, 5,0,0);

// Return it properly enwrapped:
// Note, that our matrices are row-major, but since armadillo
// is column-major, we need to transpose it first
return ArmadilloMatrix<double>(guess.t());
}
}
@@ -1,12 +1,188 @@
#include "hack_guess.hh"
#include <gint/IntegralLookup.hh>
#include <gint/version.hh>
#include <gscf/version.hh>
#include <iostream>
#include <linalgwrap/io.hh>
#include <linalgwrap/version.hh>
#include <molsturm/IopPlainScfSolver.hh>
#include <molsturm/RestrictedClosedIntegralOperator.hh>
#include <molsturm/version.hh>

namespace hard_coded_hf {} // namespace hard_coded_hf
namespace hard_coded_hf {
using namespace molsturm;
using namespace linalgwrap;

int main() {
/** Run an (atomic) SCF based on the integral data Sturmian14.
*
* \param Z Number of nuclei
* \param k_exp Sturmian exponent
* \param n_alpha Number of alpha electrons
* \param n_beta Number of beta electrons
*/
void run_rhf_sturmian_debug(double Z, double k_exp, size_t n_alpha,
size_t n_beta) {
//
// Types and settings
//
// Types of scalar and matrix
typedef double scalar_type;
typedef SmallMatrix<scalar_type> stored_matrix_type;

// Types of the integrals we use:
const gint::OrbitalType otype = gint::REAL_ATOMIC;
const std::string basis_type("Sturmians");

// The lookup class type to get the actual integrals
typedef gint::IntegralLookup<stored_matrix_type, otype> int_lookup_type;

// The type of the integral terms:
typedef typename int_lookup_type::integral_matrix_type int_term_type;

//
// Integral terms
//
// Generate integral lookup object
int_lookup_type integrals(basis_type);

// Get the integral as actual objects.
int_term_type S_bb = integrals("overlap");
int_term_type T_bb = integrals("kinetic");
int_term_type V0_bb = integrals("nuclear_attraction");
int_term_type J_bb = integrals("coulomb");
int_term_type K_bb = integrals("exchange");

// Combine 1e terms:
std::vector<int_term_type> terms_1e{std::move(T_bb), std::move(V0_bb)};

//
// Debug output
//
std::ofstream mathematicafile("/tmp/debug_molsturm_rhf_sturmian.m");
auto debugout = linalgwrap::io::make_formatted_stream_writer<
linalgwrap::io::Mathematica, scalar_type>(mathematicafile, 1e-5);

//
// Problem setup
//
// TODO hack a guess:
stored_matrix_type guess_bf =
hack_guess(static_cast<stored_matrix_type>(S_bb));
debugout.write("sbb", S_bb);
debugout.write("guess", guess_bf);

// The term container for the fock operator matrix
IntegralTermContainer<stored_matrix_type> integral_container(
std::move(terms_1e), std::move(J_bb), std::move(K_bb));

RestrictedClosedIntegralOperator<stored_matrix_type> m_fock(
integral_container, guess_bf, n_alpha, n_beta);

IopPlainScfSolver<decltype(m_fock)> solver(debugout);
solver.solve(m_fock, static_cast<stored_matrix_type>(S_bb));
}

struct args_type {
double Z = 4.0;
double k_exp = 1.0;
size_t n_alpha = 2;
size_t n_beta = n_alpha;
};

/** Quick and dirty function to parse a string to a different type.
* Return false if not possible */
template <typename T>
bool str_to_type(const std::string& in, T& out) {
std::stringstream ss(in);
if (!(ss >> out)) {
return false;
}
return true;
}

/** \brief Quick and dirty function to parse the commandline arguments
*
* \returns true if all is fine, else false
*/
bool parse_args(int argc, char** argv, args_type& parsed) {
bool had_Z = false;
bool had_alpha = false;
bool had_beta = false;
bool had_k_exp = false;

// Parsing
for (int i = 1; i < argc; ++i) {
std::string flag(argv[i]);

// Check that an argument exists:
if (++i >= argc) {
std::cerr << "Flag " << flag << " needs an argument." << std::endl;
return false;
}

// Store the argument:
std::string argument(argv[i]);

// Interpret argument
if (flag == std::string("--Z")) {
had_Z = true;
if (!str_to_type<double>(argument, parsed.Z)) {
std::cerr << "Invalid double provided to --Z: " << argument
<< std::endl;
return false;
}
} else if (flag == std::string("--alpha")) {
had_alpha = true;
if (!str_to_type<size_t>(argument, parsed.n_alpha)) {
std::cerr << "Invalid int provided to --alpha: " << argument
<< std::endl;
return false;
}
} else if (flag == std::string("--beta")) {
had_beta = true;
if (!str_to_type<size_t>(argument, parsed.n_beta)) {
std::cerr << "Invalid int provided to --beta: " << argument
<< std::endl;
return false;
}
} else if (flag == std::string("--kexp")) {
had_k_exp = true;
if (!str_to_type<double>(argument, parsed.k_exp)) {
std::cerr << "Invalid double provided to --kexp: " << argument
<< std::endl;
return false;
}
} else {
std::cerr << "Unknown flag: " << flag << std::endl;
std::cerr << "Valid are: --Z, --alpha, --beta, --kexp" << std::endl;
return false;
}
}

// Error handling
if (!had_Z) {
std::cerr << "Need flag --Z <double> to supply nuclear charge."
<< std::endl;
}
if (!had_alpha) {
std::cerr << "Need flag --alpha <int> to supply number of alpha electrons."
<< std::endl;
}
if (!had_beta) {
std::cerr << "Need flag --beta <int> to supply number of beta electrons."
<< std::endl;
}
if (!had_k_exp) {
std::cerr << "Need flag --kexp <double> to supply k exponent." << std::endl;
}

if (had_Z && had_alpha && had_beta && had_k_exp) {
return true;
}
return false;
}

int main(int argc, char** argv) {
std::cout << "molsturm version: " << molsturm::version::version_string()
<< std::endl
<< "gscf version: " << gscf::version::version_string()
@@ -16,5 +192,22 @@ int main() {
<< "linalgwrap version: " << linalgwrap::version::version_string()
<< std::endl;

args_type args;
if (!parse_args(argc, argv, args)) {
return 1;
}

std::cout << "The following configuration was read:" << std::endl
<< "Z: " << args.Z << std::endl
<< "k_exp: " << args.k_exp << std::endl
<< "n_alpha: " << args.n_alpha << std::endl
<< "n_beta: " << args.n_beta << std::endl
<< std::endl;

run_rhf_sturmian_debug(args.Z, args.k_exp, args.n_alpha, args.n_beta);
return 0;
}

} // namespace hard_coded_hf

int main(int argc, char** argv) { return hard_coded_hf::main(argc, argv); }
@@ -0,0 +1,35 @@
#pragma once
#include <linalgwrap/LazyMatrix_i.hh>

namespace molsturm {

template <typename StoredMatrix>
class IntegralOperatorBase : public linalgwrap::LazyMatrix_i<StoredMatrix> {
public:
typedef linalgwrap::LazyMatrix_i<StoredMatrix> base_type;
typedef typename base_type::scalar_type scalar_type;
typedef typename base_type::stored_matrix_type stored_matrix_type;
typedef typename base_type::size_type size_type;
typedef typename base_type::lazy_matrix_expression_ptr_type
lazy_matrix_expression_ptr_type;

// TODO implement general integral operator stuff in here
};

//@{
/** \brief struct representing a type (std::true_type, std::false_type) which
* indicates whether T is an IntegralOperator
*
* The definition is done using SFINAE, such that even for types not having a
* typedef scalar_type this expression is valid.
* */
template <typename T, typename = void>
struct IsIntegralOperator : public std::false_type {};

template <typename T>
struct IsIntegralOperator<T, linalgwrap::void_t<typename T::stored_matrix_type>>
: public std::is_base_of<
IntegralOperatorBase<typename T::stored_matrix_type>, T> {};
//@}

} // namespace molsturm
@@ -0,0 +1,52 @@
#pragma once
#include <gint/Integral.hh>
#include <linalgwrap/Constants.hh>

namespace molsturm {

template <typename StoredMatrix>
struct IntegralTermContainer {
typedef StoredMatrix stored_matrix_type;
typedef typename stored_matrix_type::scalar_type scalar_type;

//! The type to use for the integral terms.
typedef typename gint::Integral<stored_matrix_type> int_term_type;

/** The one electron terms */
std::vector<int_term_type> integral_terms_1e;

/** Set of coefficients for the one electron integral terms */
std::vector<scalar_type> coefficients_1e;

/** The coulomb term to use */
int_term_type coulomb_term;

/** The coefficient for the coulomb term */
scalar_type coefficient_coulomb;

/** The exchange term to use */
int_term_type exchange_term;

/** The coefficient for the HF exchange term */
scalar_type coefficient_exchange;

/** \brief Setup the class with the given integral terms.
* \note By default all coefficients, but the exchange coefficient, are set
* to one. The exchange term coefficient is set to -1.
* This is done in order to make sure that by default this data structure is
* ready for a Hartree-Fock calculation.
*/
IntegralTermContainer(std::vector<int_term_type> integral_terms_1e_,
int_term_type coulomb_term_,
int_term_type exchange_term_)
: integral_terms_1e{std::move(integral_terms_1e_)},
coefficients_1e{std::vector<scalar_type>(
integral_terms_1e.size(),
linalgwrap::Constants<scalar_type>::one)},
coulomb_term{std::move(coulomb_term_)},
coefficient_coulomb{linalgwrap::Constants<scalar_type>::one},
exchange_term{std::move(exchange_term_)},
coefficient_exchange{-linalgwrap::Constants<scalar_type>::one} {}
};

} // namespace molsturm
Oops, something went wrong.

0 comments on commit e264896

Please sign in to comment.