Skip to content

Adding explicit input occupancies and stabilizing virtual states in moldft #307

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
lratcliff wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Adding explicit input occupancies and stabilizing virtual states in moldft #307

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
51a4fad
Adding option for explicit input occupancies in moldft, e.g. for core…
lratcliff Jan 8, 2019
34e1d34
More stable approach to gradual decreasing virtual states in moldft
lratcliff Jan 8, 2019
a28606f
bug fix for explicit occupancy (especially for mixed AE/psp cases)
lratcliff Feb 6, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 59 additions & 5 deletions src/apps/chem/CalculationParameters.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,7 @@ struct CalculationParameters {
bool conv_only_dens; ///< If true remove bsh_residual from convergence criteria how ugly name is...
bool psp_calc; ///< pseudopotential calculation for all atoms
bool print_dipole_matels; ///< If true output dipole matrix elements
bool explicit_occ; ///< If true there should be an accompanying "occup" block which sets some input occupancies
// Next list inferred parameters
int nalpha; ///< Number of alpha spin electrons
int nbeta; ///< Number of beta spin electrons
Expand All @@ -114,11 +115,16 @@ struct CalculationParameters {
bool tdksprop; ///< time-dependent Kohn-Sham equation propagate
std::string nuclear_corrfac; ///< nuclear correlation factor
bool pure_ae; ///< pure all electron calculation with no pseudo-atoms
int nv_factor; ///< factor to multiply number of virtual orbitals with when automatically decreasing nvirt
// variables for calculations with automatically decreasing number of virtual states:
// do nv_its iters for nvalpha, then nv_its iters for nvalpha-nv_step etc. until reach nvalpha-nv_extra states
int nv_extra; ///< number of extra virtual states when automatically decreasing nvirt
int nv_step; ///< number of virtual states to decrease by
int nv_its; ///< maximum number of iterations for each number of virtual states
int vnucextra; // load balance parameter for nuclear pot.
int loadbalparts = 2; // was 6
std::string pcm_data; ///< do a PCM (solvent) calculation
std::string ac_data; ///< do a calculation with asymptotic correction (see ACParameters class in chem/AC.h for details)
std::vector<long> input_aocc, input_bocc; // input occupancies in case of explicit occupancies

// Next list for response code from a4v4
bool response; ///< response function calculation
Expand Down Expand Up @@ -150,7 +156,8 @@ struct CalculationParameters {
ar & xc_data & protocol_data;
ar & gopt & gtol & gtest & gval & gprec & gmaxiter & ginitial_hessian & algopt & tdksprop
& nuclear_corrfac & psp_calc & print_dipole_matels & pure_ae & hessian & read_cphf & restart_cphf
& purify_hessian & vnucextra & loadbalparts & pcm_data & ac_data;
& purify_hessian & vnucextra & loadbalparts & pcm_data & ac_data & explicit_occ & input_aocc & input_bocc
& nv_extra & nv_step & nv_its;
}

CalculationParameters()
Expand Down Expand Up @@ -191,6 +198,7 @@ struct CalculationParameters {
, conv_only_dens(false)
, psp_calc(false)
, print_dipole_matels(false)
, explicit_occ(false)
, nalpha(0)
, nbeta(0)
, nmo_alpha(0)
Expand All @@ -213,7 +221,9 @@ struct CalculationParameters {
, tdksprop(false)
, nuclear_corrfac("none")
, pure_ae(true)
, nv_factor(1)
, nv_extra(0)
, nv_step(0)
, nv_its(5)
, vnucextra(12)
, loadbalparts(2)
, pcm_data("none")
Expand Down Expand Up @@ -470,9 +480,18 @@ struct CalculationParameters {
else if (s == "print_dipole_matels") {
print_dipole_matels = true;
}
else if (s == "nv_factor") {
f >> nv_factor;
else if (s == "explicit_occ") {
explicit_occ = true;
}
else if (s == "nv_extra") {
f >> nv_extra;
}
else if (s == "nv_step") {
f >> nv_step;
}
else if (s == "nv_its") {
f >> nv_its;
}
else if (s == "response") {
response = true;
}
Expand Down Expand Up @@ -517,6 +536,30 @@ struct CalculationParameters {
}
if (nopen != 0) spin_restricted = false;
if (nvalpha || nvbeta) localize = false; // must use canonical orbitals if computing virtuals

// if and only if the variable is present, read the input block corresponding to the explicit occupancies
// note that everything beyond the actual number of occupied states is ignored
// i.e. virtual state occupancies are ignored
// similarly, if not enough states are specified, everything beyond that point remains unmodified
// also, if spin unrestricted, boccs will be ignored
if (explicit_occ) {
position_stream(f, "occup");
double a, b;

while (std::getline(f,s)) {
std::istringstream ss(s);
std::string tag;
if (s == "end") {
break;
}
else {
ss >> a >> b;
input_aocc.push_back(a);
input_bocc.push_back(b);
}
}
}

}

void set_molecular_info(const Molecule& molecule, const AtomicBasisSet& aobasis, unsigned int n_core) {
Expand Down Expand Up @@ -647,6 +690,17 @@ struct CalculationParameters {
madness::print(" psp or all electron ", "all electron");
else
madness::print(" psp or all electron ", "mixed psp/AE");
if (explicit_occ){
madness::print(" explicit occupancy ");
}
if (nvalpha > 0){
madness::print(" number of virtuals ", nvalpha, nvbeta);
if (nv_extra > 0 && nv_its > 0){
madness::print(" num. extra virtuals ", nv_extra);
madness::print("decrease virtuals by ", nv_step);
madness::print(" max. virtual its ", nv_its);
}
}
}

void gprint(World& world) const {
Expand Down
51 changes: 51 additions & 0 deletions src/apps/chem/SCF.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -277,6 +277,14 @@ namespace madness {
aeps = copy(aeps(Slice(n_core, n_core + param.nmo_alpha - 1)));
aocc = copy(aocc(Slice(n_core, n_core + param.nmo_alpha - 1)));
}

// override the occupancies in the case of explicit input occupancies
if (param.explicit_occ) {
int nocc = std::min<unsigned int>(param.input_aocc.size(), param.nalpha);
for (int i = 0; i < nocc; ++i){
aocc[i] = param.input_aocc[i];
}
}

if (amo[0].k() != k) {
reconstruct(world, amo);
Expand Down Expand Up @@ -332,6 +340,24 @@ namespace madness {
// normalize(world, bmo);

}

// override the occupancies in the case of explicit input occupancies
if (param.explicit_occ) {
int nocc = std::min<unsigned int>(param.input_bocc.size(), param.nbeta);
for (int i = 0; i < nocc; ++i){
bocc[i] = param.input_bocc[i];
}
}

}

// output occupancies in the case of explicit input occupancies
if (param.explicit_occ) {
if (world.rank()==0){
printf("\n\tUsing occupancies from input file:\n");
print_occs();
}

}
}

Expand Down Expand Up @@ -1257,6 +1283,14 @@ namespace madness {
aocc = tensorT(param.nmo_alpha);
for (int i = 0; i < param.nalpha; ++i)
aocc[i] = 1.0;

// override the occupancies in the case of explicit input occupancies
if (param.explicit_occ) {
int nocc = std::min<unsigned int>(param.input_aocc.size(), param.nalpha);
for (int i = 0; i < nocc; ++i){
aocc[i] = param.input_aocc[i];
}
}

if (world.rank()==0) print("grouping alpha orbitals into sets");
aset=group_orbital_sets(world,aeps,aocc,param.nmo_alpha);
Expand All @@ -1270,10 +1304,27 @@ namespace madness {
for (int i = 0; i < param.nbeta; ++i)
bocc[i] = 1.0;

// override the occupancies in the case of explicit input occupancies
if (param.explicit_occ) {
int nocc = std::min<unsigned int>(param.input_bocc.size(), param.nbeta);
for (int i = 0; i < nocc; ++i){
bocc[i] = param.input_bocc[i];
}
}

if (world.rank()==0) print("grouping beta orbitals into sets");
bset=group_orbital_sets(world,beps,bocc,param.nmo_beta);

}

// output occupancies in the case of explicit input occupancies
if (param.explicit_occ) {
if (world.rank()==0){
printf("\n\tUsing occupancies from input file:\n");
print_occs();
}
}

END_TIMER(world, "guess orbital grouping");
}
}
Expand Down
62 changes: 50 additions & 12 deletions src/apps/chem/SCF.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -407,6 +407,28 @@ namespace madness {

/// getter for the occupation numbers, alpha spin
const tensorT& get_bocc() const {return bocc;}

/// function for printing occupancies, useful in the case where explicit occupancies were given
void print_occs()
{
long naocc = aocc.size();

long nbocc = 0;
if (param.nbeta != 0 && !param.spin_restricted) {
nbocc = bocc.size();
}

// assuming aocc must always be equal to or longer than bocc
printf("\tOccupancies:\n");
for (long i = 0; i < naocc; ++i) {
printf("\t%10.6f\t", aocc(i));
if (i < nbocc){
printf("%10.6f\t", bocc(i));
}
printf("\n");
}
printf("\n");
}

bool is_spin_restricted() const {return param.spin_restricted;}

Expand Down Expand Up @@ -741,7 +763,8 @@ namespace madness {

int nvalpha = calc.param.nmo_alpha - calc.param.nalpha;
int nvbeta = calc.param.nmo_beta - calc.param.nbeta;
int nvalpha_start, nv_old;
int nvalpha_start, nvalpha_end, nv_old;
int maxiter_final = calc.param.maxiter;

// initialize the PCM solver for this geometry
if (calc.param.pcm_data != "none") {
Expand All @@ -755,19 +778,34 @@ namespace madness {

//repeat with gradually decreasing nvirt, only for first protocol
if (proto == 0 && nvalpha > 0){
nvalpha_start = nvalpha * calc.param.nv_factor;}
nvalpha_start = nvalpha;
nvalpha_end = nvalpha - calc.param.nv_extra;}
else{
nvalpha_start = nvalpha;}
nvalpha_start = nvalpha - calc.param.nv_extra;
nvalpha_end = nvalpha - calc.param.nv_extra;}

nv_old = nvalpha_start;

for (int nv=nvalpha_start;nv>=nvalpha;nv-=nvalpha){
int step;
if (calc.param.nv_extra > 0 && calc.param.nv_step == 0)
step=calc.param.nv_extra;
else
step=calc.param.nv_step;
// to make sure we eventually exit loop
if (step==0) step=1;

for (int nv=nvalpha_start;nv>=nvalpha_end;nv-=step){

if (nv > 0 && world.rank() == 0) std::cout << "Running with " << nv << " virtual states" << std::endl;

if (nv > nvalpha_end){
calc.param.maxiter = calc.param.nv_its;}
else if (nvalpha > 0){
calc.param.maxiter = maxiter_final;}

calc.param.nmo_alpha = calc.param.nalpha + nv;
// check whether this is sensible for spin restricted case
if (calc.param.nbeta && !calc.param.spin_restricted){
if (calc.param.nbeta){
if (nvbeta == nvalpha){
calc.param.nmo_beta = calc.param.nbeta + nv;}
else{
Expand Down Expand Up @@ -799,17 +837,17 @@ namespace madness {
else {
if (nv != nv_old){
calc.amo.resize(calc.param.nmo_alpha);
calc.bmo.resize(calc.param.nmo_beta);

calc.aocc = tensorT(calc.param.nmo_alpha);
for (int i = 0; i < calc.param.nalpha; ++i)
calc.aocc[i] = 1.0;

calc.bocc = tensorT(calc.param.nmo_beta);
for (int i = 0; i < calc.param.nbeta; ++i)
calc.bocc[i] = 1.0;

// might need to resize aset, bset, but for the moment this doesn't seem to be necessary

if (calc.param.nbeta){
calc.bmo.resize(calc.param.nmo_beta);

calc.bocc = tensorT(calc.param.nmo_beta);
for (int i = 0; i < calc.param.nbeta; ++i)
calc.bocc[i] = 1.0;}

}
calc.project(world);
Expand Down
⚡