Skip to content

Instantly share code, notes, and snippets.

@wpoely86

wpoely86/input.dat

Last active August 29, 2015 14:05
Show Gist options
  • Save wpoely86/99cfad8adcecd7390606 to your computer and use it in GitHub Desktop.
Save wpoely86/99cfad8adcecd7390606 to your computer and use it in GitHub Desktop.
Download ZIP
mointegrals for PSI4
Raw
input.dat
molecule BH {
B 0 0 0.0
H 0 0 2.328898308325
# symmetry c1
units au
}
set {
# exact ERI
SCF_TYPE PK
# strict convergence
E_CONVERGENCE 1e-12
D_CONVERGENCE 1e-12
ints_tolerance 1e-16
docc [3, 0, 0, 0]
basis sto-3g
}
plugin_load("./mointegrals.so")
set mointegrals {
# print 1
do_tei True
savehdf5 true
HDF5_FILENAME "mo-integrals.h5"
SAVE_WITH_SYMMETRY false
}
scf()
plugin("mointegrals.so")
Raw
mointegrals.cc
#include <libplugin/plugin.h>
#include <psi4-dec.h>
#include <libdpd/dpd.h>
#include <psifiles.h>
#include <libpsio/psio.hpp>
#include <libiwl/iwl.hpp>
#include <libtrans/integraltransform.h>
#include <libmints/wavefunction.h>
#include <libmints/mints.h>
#include <libciomr/libciomr.h>
#include <liboptions/liboptions.h>
#include <libchkpt/chkpt.h>
#include <hdf5.h>
#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
// macro to help check return status of HDF5 functions
#define HDF5_STATUS_CHECK(status) { \
if(status < 0) \
fprintf(outfile, "%s:%d: Error with HDF5. status code=%d\n", __FILE__, __LINE__, status); \
}
// This allows us to be lazy in getting the spaces in DPD calls
#define ID(x) ints.DPD_ID(x)
INIT_PLUGIN
namespace psi{ namespace mointegrals{
extern "C" int
read_options(std::string name, Options &options)
{
if (name == "MOINTEGRALS"|| options.read_globals()) {
/*- The amount of information printed
to the output file -*/
options.add_bool("PRINT_INTEGRALS", true);
/*- Whether to compute two-electron integrals -*/
options.add_bool("DO_TEI", true);
// save to a HDF5 file
options.add_bool("SAVEHDF5", false);
options.add_str("HDF5_FILENAME", "integrals.h5");
// save the matrices as a block matrix (if true),
// else save them in one big matrix (with lots of zeros)
options.add_bool("SAVE_WITH_SYMMETRY", false);
}
return true;
}
extern "C" PsiReturnType
mointegrals(Options &options)
{
/*
* This plugin shows a simple way of obtaining MO basis integrals, directly from a DPD buffer. It is also
* possible to generate integrals with labels (IWL) formatted files, but that's not shown here.
*/
bool print = options.get_bool("PRINT_INTEGRALS");
bool doTei = options.get_bool("DO_TEI");
bool savehdf5 = options.get_bool("SAVEHDF5");
std::string filename = options.get_str("HDF5_FILENAME");
boost::algorithm::to_lower(filename);
bool savewithsym = options.get_bool("SAVE_WITH_SYMMETRY");
if(savewithsym)
{
fprintf(outfile, "ERROR: SAVE_WITH_SYMMETRY does not yet work for TEI");
return Failure;
}
if(savewithsym && !savehdf5)
{
fprintf(outfile, "ERROR: SAVE_WITH_SYMMETRY is set but SAVEHDF5 is not. This does make any sense?");
return Failure;
}
// Grab the global (default) PSIO object, for file I/O
boost::shared_ptr<PSIO> psio(_default_psio_lib_);
// Now we want the reference (SCF) wavefunction
boost::shared_ptr<Wavefunction> wfn = Process::environment.wavefunction();
if(!wfn) throw PSIEXCEPTION("SCF has not been run yet!");
// Quickly check that there are no open shell orbitals here...
int nirrep = wfn->nirrep();
// For now, we'll just transform for closed shells and generate all integrals. For more elaborate use of the
// LibTrans object, check out the plugin_mp2 example.
std::vector<boost::shared_ptr<MOSpace> > spaces;
spaces.push_back(MOSpace::all);
IntegralTransform ints(wfn, spaces, IntegralTransform::Restricted);
ints.transform_tei(MOSpace::all, MOSpace::all, MOSpace::all, MOSpace::all);
// Use the IntegralTransform object's DPD instance, for convenience
dpd_set_default(ints.get_dpd_id());
//Readin the MO OEI in moOei & print everything
int nmo = Process::environment.wavefunction()->nmo();
int nIrreps = Process::environment.wavefunction()->nirrep();
int *orbspi = Process::environment.wavefunction()->nmopi();
int *docc = Process::environment.wavefunction()->doccpi();
int *socc = Process::environment.wavefunction()->soccpi();
int nTriMo = nmo * (nmo + 1) / 2;
double *temp = new double[nTriMo];
Matrix moOei("MO OEI", nIrreps, orbspi, orbspi);
IWL::read_one(psio.get(), PSIF_OEI, PSIF_MO_OEI, temp, nTriMo, 0, 0, outfile);
moOei.set(temp);
moOei.print();
fprintf(outfile, "**** Molecular Integrals Start Here \n");
if(savehdf5)
{
fprintf(outfile, "**** Saving integrals to %s \n", filename.c_str());
if(savewithsym)
fprintf(outfile, "**** Saving to HDF5 as block matrices \n");
}
std::string SymmLabel = Process::environment.molecule()->sym_label();
fprintf(outfile, "Symmetry Label = ");
fprintf(outfile, SymmLabel.c_str());
fprintf(outfile, " \n");
fprintf(outfile, "Nirreps = %1d \n", nirrep);
double NuclRepulsion = Process::environment.molecule()->nuclear_repulsion_energy();
fprintf(outfile, "Nuclear Repulsion Energy = %16.48f \n", NuclRepulsion);
fprintf(outfile, "Number Of Molecular Orbitals = %2d \n", nmo);
fprintf(outfile, "Irreps Of Molecular Orbitals = \n");
for (int h=0; h<nirrep; ++h){
for (int cnt=0; cnt<moOei.rowspi(h); ++cnt){
fprintf(outfile, "%1d ",h);
}
}
fprintf(outfile, "\n");
fprintf(outfile, "DOCC = ");
for (int h=0; h<nirrep; h++){
fprintf(outfile, "%2d ",docc[h]);
}
fprintf(outfile, "\n");
fprintf(outfile, "SOCC = ");
for (int h=0; h<nirrep; h++){
fprintf(outfile, "%2d ",socc[h]);
}
fprintf(outfile, "\n");
int nelectrons = 0;
for(int i=0;i<Process::environment.molecule()->natom();i++)
nelectrons += Process::environment.molecule()->true_atomic_number(i);
nelectrons -= Process::environment.molecule()->molecular_charge();
fprintf(outfile, "**** MO OEI \n");
boost::shared_ptr<Matrix> hMat;
if(savehdf5 && !savewithsym)
{
hMat.reset(new Matrix(nmo, nmo));
hMat->set(0.0);
}
int nTot = 0;
for(int h = 0; h < nirrep; ++h){
for(int cnt = 0; cnt < moOei.rowspi(h); ++cnt){
for (int cnt2 = 0; cnt2 < moOei.colspi(h); ++cnt2)
{
fprintf(outfile, "%1d\t%1d\t%16.48f \n",
nTot+cnt, nTot+cnt2, moOei[h][cnt][cnt2]);
if(savehdf5 && !savewithsym)
{
hMat->set(nTot+cnt, nTot+cnt2, moOei[h][cnt][cnt2]);
hMat->set(nTot+cnt2, nTot+cnt, moOei[h][cnt][cnt2]);
}
}
}
nTot += moOei.rowspi(h);
}
hid_t file_id, group_id, dataset_id, dataspace_id, attribute_id;
herr_t status;
if(savehdf5)
{
file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
HDF5_STATUS_CHECK(file_id);
group_id = H5Gcreate(file_id, "/integrals", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
HDF5_STATUS_CHECK(group_id);
dataspace_id = H5Screate(H5S_SCALAR);
attribute_id = H5Acreate (group_id, "nelectrons", H5T_STD_I32LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
status = H5Awrite (attribute_id, H5T_NATIVE_INT, &nelectrons );
HDF5_STATUS_CHECK(status);
status = H5Aclose(attribute_id);
HDF5_STATUS_CHECK(status);
attribute_id = H5Acreate (group_id, "sp_dim", H5T_STD_I32LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
status = H5Awrite (attribute_id, H5T_NATIVE_INT, &nmo );
HDF5_STATUS_CHECK(status);
status = H5Aclose(attribute_id);
HDF5_STATUS_CHECK(status);
attribute_id = H5Acreate (group_id, "nuclear_repulsion_energy", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
status = H5Awrite (attribute_id, H5T_NATIVE_DOUBLE, &NuclRepulsion );
HDF5_STATUS_CHECK(status);
status = H5Aclose(attribute_id);
HDF5_STATUS_CHECK(status);
status = H5Sclose(dataspace_id);
HDF5_STATUS_CHECK(status);
if(savewithsym)
{
hid_t oei_group_id;
oei_group_id = H5Gcreate(group_id, "OEI", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
HDF5_STATUS_CHECK(oei_group_id);
for(int b=0;b<nirrep;b++)
{
hsize_t dimarray = moOei.rowspi(b) * moOei.colspi(b);
dataspace_id = H5Screate_simple(1, &dimarray, NULL);
std::string block_name = boost::lexical_cast<std::string>(b);
dataset_id = H5Dcreate(oei_group_id, block_name.c_str(), H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if(dimarray)
status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, moOei.get_pointer(b) );
HDF5_STATUS_CHECK(status);
status = H5Dclose(dataset_id);
HDF5_STATUS_CHECK(status);
status = H5Sclose(dataspace_id);
HDF5_STATUS_CHECK(status);
}
status = H5Gclose(oei_group_id);
HDF5_STATUS_CHECK(status);
} else {
hsize_t dimarray = nmo * nmo;
dataspace_id = H5Screate_simple(1, &dimarray, NULL);
dataset_id = H5Dcreate(group_id, "OEI", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, hMat->get_pointer(0) );
HDF5_STATUS_CHECK(status);
status = H5Dclose(dataset_id);
HDF5_STATUS_CHECK(status);
status = H5Sclose(dataspace_id);
HDF5_STATUS_CHECK(status);
}
}
fprintf(outfile, "**** MO TEI \n");
boost::shared_ptr<Matrix> VMat;
if(savehdf5 && !savewithsym)
{
VMat.reset(new Matrix(nmo * nmo, nmo * nmo));
VMat->set(0.0);
}
hid_t tei_group_id;
if(savehdf5 && savewithsym)
{
tei_group_id = H5Gcreate(group_id, "TEI", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
HDF5_STATUS_CHECK(tei_group_id);
}
/*
* Now, loop over the DPD buffer, printing the integrals
*/
dpdbuf4 K;
psio->open(PSIF_LIBTRANS_DPD, PSIO_OPEN_OLD);
// To only process the permutationally unique integrals, change the ID("[A,A]") to ID("[A>=A]+")
global_dpd_->buf4_init(&K, PSIF_LIBTRANS_DPD, 0, ID("[A,A]"), ID("[A,A]"),
ID("[A>=A]+"), ID("[A>=A]+"), 0, "MO Ints (AA|AA)");
for(int h = 0; h < nirrep; ++h){
global_dpd_->buf4_mat_irrep_init(&K, h);
global_dpd_->buf4_mat_irrep_rd(&K, h);
boost::shared_ptr<Matrix> block_vMat(new Matrix(K.params->rowtot[h], K.params->coltot[h]));
assert(K.params->rowtot[h] == K.params->coltot[h]);
for(int pq = 0; pq < K.params->rowtot[h]; ++pq){
int p = K.params->roworb[h][pq][0];
int q = K.params->roworb[h][pq][1];
int psym = K.params->psym[p];
int qsym = K.params->qsym[q];
int prel = p - K.params->poff[psym];
int qrel = q - K.params->qoff[qsym];
for(int rs = 0; rs < K.params->coltot[h]; ++rs){
int r = K.params->colorb[h][rs][0];
int s = K.params->colorb[h][rs][1];
int rsym = K.params->rsym[r];
int ssym = K.params->ssym[s];
int rrel = r - K.params->roff[rsym];
int srel = s - K.params->soff[ssym];
// Print out the absolute orbital numbers, the relative (within irrep)
// numbers, the symmetries, and the integral itself
fprintf(outfile, "(%2d %2d | %2d %2d) = %16.48f, "
"symmetries = (%1d %1d | %1d %1d), "
"relative indices = (%2d %2d | %2d %2d)\n",
p, q, r, s, K.matrix[h][pq][rs],
psym, qsym, rsym, ssym,
prel, qrel, rrel, srel);
// PSI is chemical notation, we get: (pr|V|qs)
// so: V_abcd = V(a*nmo+b,c*nmo+d)
if(savehdf5 && !savewithsym)
{
int idx1 = p * nmo + r;
int idx2 = q * nmo + s;
VMat->set(idx1, idx2, K.matrix[h][pq][rs]);
VMat->set(idx2, idx1, K.matrix[h][pq][rs]);
} else if(savehdf5)
{
int idx1 = prel * orbspi[h] + rrel;
int idx2 = qrel * orbspi[h] + srel;
block_vMat->set(idx1, idx2, K.matrix[h][pq][rs]);
block_vMat->set(idx2, idx1, K.matrix[h][pq][rs]);
}
}
}
if(savehdf5 && savewithsym)
{
hsize_t dimarray = block_vMat->rowspi(0) * block_vMat->colspi(0);
dataspace_id = H5Screate_simple(1, &dimarray, NULL);
std::string block_name = boost::lexical_cast<std::string>(h);
dataset_id = H5Dcreate(tei_group_id, block_name.c_str(), H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if(dimarray)
status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, block_vMat->get_pointer(0) );
HDF5_STATUS_CHECK(status);
status = H5Dclose(dataset_id);
HDF5_STATUS_CHECK(status);
status = H5Sclose(dataspace_id);
HDF5_STATUS_CHECK(status);
}
global_dpd_->buf4_mat_irrep_close(&K, h);
}
global_dpd_->buf4_close(&K);
psio->close(PSIF_LIBTRANS_DPD, PSIO_OPEN_OLD);
if(savehdf5 && !savewithsym)
{
hsize_t dimarray = nmo * nmo * nmo * nmo;
dataspace_id = H5Screate_simple(1, &dimarray, NULL);
dataset_id = H5Dcreate(group_id, "TEI", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, VMat->get_pointer(0) );
HDF5_STATUS_CHECK(status);
status = H5Dclose(dataset_id);
HDF5_STATUS_CHECK(status);
status = H5Sclose(dataspace_id);
HDF5_STATUS_CHECK(status);
} else if(savehdf5)
{
status = H5Gclose(tei_group_id);
HDF5_STATUS_CHECK(status);
}
if(savehdf5)
{
status = H5Gclose(group_id);
HDF5_STATUS_CHECK(status);
status = H5Fclose(file_id);
HDF5_STATUS_CHECK(status);
}
return Success;
}
}} // End Namespaces
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment