Last active
August 29, 2015 14:05
-
-
Save wpoely86/99cfad8adcecd7390606 to your computer and use it in GitHub Desktop.
mointegrals for PSI4
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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