diff --git a/share/python/driver.py b/share/python/driver.py index 6881152ad26..a12f9bc98bc 100644 --- a/share/python/driver.py +++ b/share/python/driver.py @@ -98,7 +98,7 @@ 'cc3' : run_ccenergy, 'mrcc' : run_mrcc, # interface to Kallay's MRCC program 'bccd' : run_bccd, - 'bccd(t)' : run_bccd_t, + 'bccd(t)' : run_bccd, 'eom-ccsd' : run_eom_cc, 'eom-cc2' : run_eom_cc, 'eom-cc3' : run_eom_cc, diff --git a/share/python/proc.py b/share/python/proc.py index 020e43885a3..22134ce49cb 100644 --- a/share/python/proc.py +++ b/share/python/proc.py @@ -1837,11 +1837,25 @@ def run_bccd(name, **kwargs): psi4.set_local_option('CCTRANSORT', 'WFN', 'BCCD') psi4.set_local_option('CCENERGY', 'WFN', 'BCCD') + elif (name.lower() == 'bccd(t)'): + psi4.set_local_option('TRANSQT2', 'WFN', 'BCCD_T') + psi4.set_local_option('CCSORT', 'WFN', 'BCCD_T') + psi4.set_local_option('CCENERGY', 'WFN', 'BCCD_T') + psi4.set_local_option('CCTRANSORT', 'WFN', 'BCCD_T') + psi4.set_local_option('CCTRIPLES', 'WFN', 'BCCD_T') + else: + raise ValidationError("proc.py:run_bccd name %s not recognized" % name) + + # Bypass routine scf if user did something special to get it to converge ref_wfn = kwargs.get('ref_wfn', None) if ref_wfn is None: ref_wfn = scf_helper(name, **kwargs) # C1 certified + # Needed for (T). + if (psi4.get_option('SCF', 'REFERENCE') == 'ROHF'): + ref_wfn.semicanonicalize() + if psi4.get_option('SCF', 'SCF_TYPE') in ['CD', 'DF']: mints = psi4.MintsHelper(ref_wfn.molecule().basisset()) mints.integrals() @@ -1856,41 +1870,23 @@ def run_bccd(name, **kwargs): else: psi4.transqt2(ref_wfn) psi4.ccsort() + ref_wfn = psi4.ccenergy(ref_wfn) - psi4.print_out('Brueckner convergence check: %d\n' % psi4.get_variable('BRUECKNER CONVERGED')) + psi4.print_out('Brueckner convergence check: %s\n' % bool(psi4.get_variable('BRUECKNER CONVERGED'))) if (psi4.get_variable('BRUECKNER CONVERGED') == True): break - # This is quite arbitrary, blame DGAS - if bcc_iter_cnt > 10: - raise ValidationError("Stopped at 50th BCCSD iteration, something is going wrong!") + if bcc_iter_cnt >= psi4.get_option('CCENERGY', 'BCCD_MAXITER'): + psi4.print_out("\n\nWarning! BCCD did not converge within the maximum number of iterations.") + psi4.print_out("You can increase the number of BCCD iterations by changing BCCD_MAXITER.\n\n") break bcc_iter_cnt += 1 - optstash.restore() - return ref_wfn - - -def run_bccd_t(name, **kwargs): - """Function encoding sequence of PSI module calls for - a Brueckner CCD(T) calculation. - - """ - optstash = p4util.OptionsState( - ['TRANSQT2', 'WFN'], - ['CCSORT', 'WFN'], - ['CCENERGY', 'WFN'], - ['CCTRIPLES', 'WFN']) - - psi4.set_local_option('TRANSQT2', 'WFN', 'BCCD_T') - psi4.set_local_option('CCSORT', 'WFN', 'BCCD_T') - psi4.set_local_option('CCENERGY', 'WFN', 'BCCD_T') - psi4.set_local_option('CCTRIPLES', 'WFN', 'BCCD_T') - bccd_wfn = run_bccd(name, **kwargs) - psi4.cctriples(bccd_wfn) + if (name.lower() == 'bccd(t)'): + psi4.cctriples(ref_wfn) optstash.restore() - return bccd_wfn + return ref_wfn def run_scf_property(name, **kwargs): diff --git a/src/bin/ccenergy/get_params.cc b/src/bin/ccenergy/get_params.cc index 8338ba53d8c..963baf56dda 100644 --- a/src/bin/ccenergy/get_params.cc +++ b/src/bin/ccenergy/get_params.cc @@ -64,7 +64,8 @@ void CCEnergyWavefunction::get_params(Options &options) else if(junk == "ROHF" && (params_.wfn == "MP2" || params_.wfn == "CCSD_T" || params_.wfn == "CCSD_AT" || params_.wfn == "CC3" || params_.wfn == "EOM_CC3" || - params_.wfn == "CC2" || params_.wfn == "EOM_CC2")) { + params_.wfn == "CC2" || params_.wfn == "EOM_CC2" || + params_.wfn == "BCCD" || params_.wfn == "BCCD_T")) { params_.ref = 2; params_.semicanonical = 1; } diff --git a/src/bin/cctransort/cctransort.cc b/src/bin/cctransort/cctransort.cc index b23d9bb7df3..e528cec8c48 100644 --- a/src/bin/cctransort/cctransort.cc +++ b/src/bin/cctransort/cctransort.cc @@ -83,7 +83,8 @@ PsiReturnType cctransort(SharedWavefunction ref, Options& options) else if(options.get_str("REFERENCE") =="ROHF" && (options.get_str("WFN")=="MP2" || options.get_str("WFN")=="CCSD_T" || options.get_str("WFN")=="CCSD_AT" || options.get_str("WFN")=="CC3" || options.get_str("WFN")=="EOM_CC3" || - options.get_str("WFN")=="CC2" || options.get_str("WFN")=="EOM_CC2")) { + options.get_str("WFN")=="CC2" || options.get_str("WFN")=="EOM_CC2" || + options.get_str("WFN")=="BCCD" || options.get_str("WFN")=="BCCD_T")) { reference = 2; semicanonical = true; } @@ -97,7 +98,7 @@ PsiReturnType cctransort(SharedWavefunction ref, Options& options) // Allow user to force semicanonical if(options["SEMICANONICAL"].has_changed()) { semicanonical = options.get_bool("SEMICANONICAL"); - reference = 2; + if (semicanonical) reference = 2; } int nirreps = ref->nirrep(); diff --git a/src/bin/cctriples/get_moinfo.cc b/src/bin/cctriples/get_moinfo.cc index b5cdaf64d7f..f28cc46139f 100644 --- a/src/bin/cctriples/get_moinfo.cc +++ b/src/bin/cctriples/get_moinfo.cc @@ -84,7 +84,7 @@ void get_moinfo(boost::shared_ptr wfn, Options &options) junk = options.get_str("REFERENCE"); /* if no reference is given, assume rhf */ if(junk == "RHF") params.ref = 0; - else if(junk == "ROHF" && params.wfn == "CCSD_T") { + else if(junk == "ROHF" && (params.wfn == "CCSD_T" || params.wfn == "BCCD_T")) { params.ref = 2; params.semicanonical = 1; } diff --git a/src/bin/psi4/read_options.cc b/src/bin/psi4/read_options.cc index 309f6b9ee8e..e44f45ec43a 100644 --- a/src/bin/psi4/read_options.cc +++ b/src/bin/psi4/read_options.cc @@ -2194,6 +2194,8 @@ int read_options(const std::string &name, Options & options, bool suppress_print options.add_double("CC_SS_SCALE",1.13); /*- Convert ROHF MOs to semicanonical MOs -*/ options.add_bool("SEMICANONICAL", true); + /*- Convert ROHF MOs to semicanonical MOs -*/ + options.add_int("BCCD_MAXITER", 50); } if(name == "DFMP2"|| options.read_globals()) { /*- MODULEDESCRIPTION Performs density-fitted MP2 computations for RHF/UHF/ROHF reference wavefunctions. -*/ diff --git a/src/bin/transqt2/get_params.cc b/src/bin/transqt2/get_params.cc index ba0281902f0..7e3526fe8e9 100644 --- a/src/bin/transqt2/get_params.cc +++ b/src/bin/transqt2/get_params.cc @@ -58,7 +58,8 @@ void get_params(Options & options) else if((reference == "ROHF") and ((params.wfn == "MP2") or (params.wfn == "CCSD_T") or (params.wfn == "CCSD_AT") or (params.wfn == "CC3") or (params.wfn == "EOM_CC3") or - (params.wfn == "CC2") or (params.wfn == "EOM_CC2"))) { + (params.wfn == "CC2") or (params.wfn == "EOM_CC2") or + (params.wfn == "BCCD") or (params.wfn == "BCCD_T"))) { params.ref = 2; params.semicanonical = 1; } diff --git a/tests/cc16/input.dat b/tests/cc16/input.dat index 8e8b9e38473..8a7177ae762 100644 --- a/tests/cc16/input.dat +++ b/tests/cc16/input.dat @@ -1,4 +1,4 @@ -#! UHF-B-CCD(T)/cc-pVDZ $^{3}B@@1$ CH2 single-point energy (fzc, MO-basis $\langle ab|cd \rangle$ ) +#! ROHF and UHF-B-CCD(T)/cc-pVDZ $^{3}B@@1$ CH2 single-point energy (fzc, MO-basis $\langle ab|cd \rangle$ ) memory 250 mb @@ -12,9 +12,9 @@ molecule ch2 { } set { - reference = uhf - basis cc-pVDZ - freeze_core true + reference uhf + basis cc-pVDZ + freeze_core true } energy("bccd(t)") @@ -25,3 +25,18 @@ ebccd_t = -39.032691827829460 #TEST compare_values(escf, get_variable("SCF TOTAL ENERGY"), 7, "SCF energy") #TEST compare_values(ebccd, get_variable("CCSD TOTAL ENERGY"), 7, "B-CCD energy") #TEST compare_values(ebccd_t, get_variable("CCSD(T) TOTAL ENERGY"), 7, "B-CCD(T) energy") #TEST + + +# We should obtain the same result as above, but start with different orbitals +# Energy will be slightly different as the core is frozen +set { + reference rohf +} + +energy("bccd(t)") +escf = -38.91341670976116 #TEST +ebccd = -39.030807046983838 #TEST +ebccd_t = -39.032665163463861 #TEST +compare_values(escf, get_variable("SCF TOTAL ENERGY"), 7, "SCF energy") #TEST +compare_values(ebccd, get_variable("CCSD TOTAL ENERGY"), 7, "B-CCD energy") #TEST +compare_values(ebccd_t, get_variable("CCSD(T) TOTAL ENERGY"), 7, "B-CCD(T) energy") #TEST