use std::process::ExitCode;
use hartree::core::Molecule;
use hartree::core::units::{ANGSTROM_TO_BOHR, AU_DIPOLE_TO_DEBYE};
use hartree::dft::FunctionalSpec;
use hartree::disp::Dispersion;
use hartree::opt::OptResult;
use hartree::opt::ts::{TsResult, TsStatus};
use hartree::props::frequencies::FrequencyResult;
use hartree::props::population::PopulationAnalysis;
use hartree::props::thermo::ThermoResult;
use hartree::scf::{Reference, ScfResult, Smearing};
use hartree::{DftDiagnostics, Job, JobOptions, Method, PostHfResult};
mod periodic;
mod ts_flags;
fn main() -> ExitCode {
match run() {
Ok(true) => ExitCode::SUCCESS,
Ok(false) => ExitCode::FAILURE,
Err(message) => {
eprintln!("error: {message}");
ExitCode::FAILURE
}
}
}
fn run() -> Result<bool, String> {
let args: Vec<String> = std::env::args().skip(1).collect();
if args.is_empty() || args.iter().any(|a| a == "-h" || a == "--help") {
print_usage();
return Ok(true);
}
if periodic::is_periodic(&args) {
return periodic::run_periodic_cli(&args);
}
let mut xyz_path: Option<String> = None;
let mut basis = String::from("sto-3g");
let mut basis_explicit = false;
let mut method_str = String::from("rhf");
let mut charge: i32 = 0;
let mut multiplicity: u32 = 1;
let mut do_opt = false;
let mut do_ts = false;
let mut ts_options = hartree::opt::ts::TsOptions::default();
let mut ts_flag_seen = false;
let mut all_electron = false;
let mut direct = false;
let mut ri = false;
let mut ri_mp2 = false;
let mut cosx = false;
let mut x2c = false;
let mut do_properties = false;
let mut do_freq = false;
let mut do_sph = false;
let mut conformers: Option<Option<String>> = None;
let mut conformers_out: Option<String> = None;
let mut rerank = false;
let mut symmetry_number: u32 = 1;
let mut qrrho_w0: f64 = hartree::props::thermo::QRRHO_W0_DEFAULT_CM1;
let mut grid: Option<usize> = None;
let mut solvent: Option<String> = None;
let mut eps: Option<f64> = None;
let mut smd: Option<String> = None;
let mut alpb: Option<String> = None;
let mut gbsa: Option<String> = None;
let mut cosmo_file: Option<String> = None;
let mut smear: Option<f64> = None;
let mut method_explicit = false;
let mut fod = false;
let mut fod_cube: Option<String> = None;
let mut cp_na: Option<usize> = None;
let mut cp_charges: Option<(i32, i32)> = None;
let mut cp_mults: Option<(u32, u32)> = None;
let mut gcp_keyword: Option<String> = None;
let mut no_method_warnings = false;
let mut recommend_task: Option<String> = None;
let mut protocol: Option<String> = None;
let mut i = 0;
while i < args.len() {
match args[i].as_str() {
"--basis" => {
basis = take(&args, &mut i, "--basis")?;
basis_explicit = true;
}
"--method" => {
method_str = take(&args, &mut i, "--method")?.to_ascii_lowercase();
method_explicit = true;
}
"--charge" => {
charge = take(&args, &mut i, "--charge")?
.parse()
.map_err(|_| "--charge must be an integer".to_string())?
}
"--spin" | "--mult" => {
multiplicity = take(&args, &mut i, "--spin")?
.parse()
.map_err(|_| "--spin must be a positive integer (2S+1)".to_string())?
}
"--opt" => do_opt = true,
"--ts" => do_ts = true,
"--all-electron" => all_electron = true,
"--direct" => direct = true,
"--ri" => ri = true,
"--ri-mp2" => ri_mp2 = true,
"--cosx" => cosx = true,
"--x2c" => x2c = true,
"--properties" | "--props" => do_properties = true,
"--freq" => do_freq = true,
"--sph" => do_sph = true,
"--conformers" => {
if args.get(i + 1).map(|s| s.as_str()) == Some("crest") {
i += 1;
conformers = Some(Some("crest".to_string()));
} else {
conformers = Some(None);
}
}
"--conformers-out" => conformers_out = Some(take(&args, &mut i, "--conformers-out")?),
"--rerank" => rerank = true,
"--qrrho-w0" => {
qrrho_w0 = take(&args, &mut i, "--qrrho-w0")?
.parse()
.map_err(|_| "--qrrho-w0 must be a frequency in cm^-1".to_string())?;
if !qrrho_w0.is_finite() || qrrho_w0 <= 0.0 {
return Err("--qrrho-w0 must be a positive frequency in cm^-1".into());
}
}
"--symmetry-number" => {
symmetry_number = take(&args, &mut i, "--symmetry-number")?
.parse()
.map_err(|_| "--symmetry-number must be a positive integer".to_string())?
}
"--grid" => {
grid = Some(
take(&args, &mut i, "--grid")?
.parse()
.map_err(|_| "--grid must be an integer 0..=4".to_string())?,
)
}
"--solvent" => solvent = Some(take(&args, &mut i, "--solvent")?.to_ascii_lowercase()),
"--smd" => smd = Some(take(&args, &mut i, "--smd")?.to_ascii_lowercase()),
"--alpb" => alpb = Some(take(&args, &mut i, "--alpb")?.to_ascii_lowercase()),
"--gbsa" => gbsa = Some(take(&args, &mut i, "--gbsa")?.to_ascii_lowercase()),
"--cosmo-file" => cosmo_file = Some(take(&args, &mut i, "--cosmo-file")?),
"--smear" => {
let t: f64 = take(&args, &mut i, "--smear")?
.parse()
.map_err(|_| "--smear must be a temperature in kelvin".to_string())?;
if !t.is_finite() || t <= 0.0 {
return Err("--smear temperature must be positive".into());
}
smear = Some(t);
}
"--cp" => {
let n: usize = take(&args, &mut i, "--cp")?.parse().map_err(|_| {
"--cp must be the number of atoms in fragment A (the first n atoms \
of the XYZ, in input order)"
.to_string()
})?;
if n == 0 {
return Err("--cp needs at least one atom in fragment A".into());
}
cp_na = Some(n);
}
"--cp-charges" => {
let v = take(&args, &mut i, "--cp-charges")?;
cp_charges = Some(parse_pair(&v, "--cp-charges")?);
}
"--cp-mults" => {
let v = take(&args, &mut i, "--cp-mults")?;
let (a, b): (u32, u32) = parse_pair(&v, "--cp-mults")?;
if a == 0 || b == 0 {
return Err("--cp-mults multiplicities must be >= 1".into());
}
cp_mults = Some((a, b));
}
"--gcp" => gcp_keyword = Some(take(&args, &mut i, "--gcp")?.to_ascii_lowercase()),
"--protocol" => {
protocol = Some(take(&args, &mut i, "--protocol")?.to_ascii_lowercase())
}
"--no-method-warnings" => no_method_warnings = true,
"--recommend" => {
recommend_task = Some(take(&args, &mut i, "--recommend")?.to_ascii_lowercase())
}
"--fod" => fod = true,
"--fod-cube" => fod_cube = Some(take(&args, &mut i, "--fod-cube")?),
"--eps" => {
eps = Some(
take(&args, &mut i, "--eps")?
.parse()
.map_err(|_| "--eps must be a number (dielectric constant)".to_string())?,
)
}
"--ts-irc" => {
ts_options.confirm_irc = true;
ts_flag_seen = true;
}
f if ts_flags::TS_VALUE_FLAGS.contains(&f) => {
let flag = args[i].clone();
let value = take(&args, &mut i, &flag)?;
ts_flags::apply_ts_option(&mut ts_options, &flag, &value)?;
ts_flag_seen = true;
}
other if other.starts_with("--") => return Err(format!("unknown option {other}")),
path => xyz_path = Some(path.to_string()),
}
i += 1;
}
if ts_flag_seen && !do_ts {
eprintln!("warning: --ts-* flags are ignored without --ts");
}
if let Some(task) = &recommend_task {
return report_recommendation(task);
}
let xyz_path = xyz_path.ok_or("missing XYZ file argument (try --help)")?;
let xyz = std::fs::read_to_string(&xyz_path).map_err(|e| format!("reading {xyz_path}: {e}"))?;
let molecule = Molecule::from_xyz(&xyz)
.map_err(|e| e.to_string())?
.with_charge(charge)
.with_multiplicity(multiplicity);
molecule.validate().map_err(|e| e.to_string())?;
if do_sph && !do_freq {
return Err("--sph modifies --freq; pass --freq as well".into());
}
let solvent_eps = match (&solvent, eps) {
(Some(_), Some(_)) => {
return Err(
"--solvent and --eps are mutually exclusive (a named solvent already fixes ε)"
.into(),
);
}
(Some(name), None) => Some(hartree::solv::solvent_epsilon(name).ok_or_else(|| {
let names: Vec<&str> = hartree::solv::SOLVENTS.iter().map(|(n, _)| *n).collect();
format!(
"unknown solvent {name:?} (available: {}; or give --eps <value> directly)",
names.join(", ")
)
})?),
(None, Some(e)) => Some(e),
(None, None) => None,
};
if smd.is_some() && (solvent.is_some() || eps.is_some()) {
return Err(
"--smd and --solvent/--eps are mutually exclusive (SMD fixes both the dielectric constant and the non-electrostatic parameterization)"
.into(),
);
}
if let Some(name) = &smd
&& hartree::solv::smd_solvent(name).is_none()
{
let names: Vec<&str> = hartree::solv::SMD_SOLVENTS.iter().map(|s| s.name).collect();
return Err(format!(
"unknown SMD solvent {name:?} (available: {})",
names.join(", ")
));
}
let n_solv = [
solvent.is_some() || eps.is_some(),
smd.is_some(),
alpb.is_some(),
gbsa.is_some(),
cosmo_file.is_some(),
]
.iter()
.filter(|&&x| x)
.count();
if n_solv > 1 {
return Err(
"the solvation options are mutually exclusive: choose at most one of \
--solvent/--eps (C-PCM), --smd, --alpb, --gbsa, or --cosmo-file"
.into(),
);
}
if let Some(name) = &alpb
&& hartree::solv::alpb_solvent(name).is_none()
{
return Err(format!(
"unknown ALPB solvent {name:?} (available: {})",
hartree::solv::alpb_solvent_names().join(", ")
));
}
if let Some(name) = &gbsa
&& hartree::solv::gbsa_solvent(name).is_none()
{
return Err(format!(
"unknown GBSA solvent {name:?} (available: {})",
hartree::solv::gbsa_solvent_names().join(", ")
));
}
if let Some(name) = &protocol {
if !matches!(name.as_str(), "w1" | "hartree-w1") {
return Err(format!(
"unknown protocol {name:?} (available: w1 — the hartree-W1 composite \
thermochemistry protocol)"
));
}
if method_explicit || basis_explicit {
return Err(
"--protocol w1 defines its own methods and bases (B3LYP/cc-pVTZ geometry \
+ frequencies; CCSD(T)/cc-pVTZ; CCSD/cc-pVQZ; CBS extrapolation); drop \
--method/--basis"
.into(),
);
}
if do_opt || do_freq {
return Err(
"--protocol w1 already optimizes the geometry and computes frequencies as \
protocol stages; drop --opt/--freq"
.into(),
);
}
for (on, what) in [
(direct, "--direct"),
(ri, "--ri"),
(ri_mp2, "--ri-mp2"),
(cosx, "--cosx"),
(x2c, "--x2c"),
(fod, "--fod"),
(smear.is_some(), "--smear"),
(do_properties, "--properties"),
(do_sph, "--sph"),
(conformers.is_some(), "--conformers"),
(rerank, "--rerank"),
(cosmo_file.is_some(), "--cosmo-file"),
(gcp_keyword.is_some(), "--gcp"),
(cp_na.is_some(), "--cp"),
(grid.is_some(), "--grid"),
(solvent_eps.is_some(), "--solvent/--eps"),
(smd.is_some(), "--smd"),
(alpb.is_some(), "--alpb"),
(gbsa.is_some(), "--gbsa"),
] {
if on {
return Err(format!(
"{what} is not supported with --protocol w1 (the protocol is a fixed \
gas-phase staged workflow); run the stages as separate jobs instead"
));
}
}
let w1_opts = hartree::w1::W1Options {
symmetry_number,
qrrho_w0_cm1: qrrho_w0,
all_electron,
..hartree::w1::W1Options::default()
};
let res = hartree::w1::run_w1(&molecule, &w1_opts)?;
report_w1(&molecule, &w1_opts, &res);
return Ok(true);
}
let ml_spec = hartree::multilevel::parse_spec(&method_str, multiplicity)?;
if rerank && ml_spec.is_none() {
return Err(
"--rerank requires --conformers and a multi-level '//' --method spec \
(e.g. --method \"r2scan-3c // hf/sto-3g\")"
.into(),
);
}
if let Some(spec) = ml_spec {
if basis_explicit {
return Err(
"--basis conflicts with a multi-level '//' spec: each level carries its \
basis inside the spec itself (method/basis, or a composite keyword)"
.into(),
);
}
if do_opt {
return Err(
"--opt is redundant with a multi-level '//' spec: the OPT level (right of \
'//') already optimizes the geometry; drop --opt"
.into(),
);
}
for (on, what) in [
(direct, "--direct"),
(ri, "--ri"),
(cosx, "--cosx"),
(x2c, "--x2c"),
(fod, "--fod"),
(smear.is_some(), "--smear"),
(do_properties, "--properties"),
(do_sph, "--sph"),
(cosmo_file.is_some(), "--cosmo-file"),
(gcp_keyword.is_some(), "--gcp"),
(cp_na.is_some(), "--cp"),
] {
if on {
return Err(format!(
"{what} is not supported with a multi-level '//' spec; \
run the stages as separate single-level jobs instead"
));
}
}
if let Some(g) = grid
&& g > 4
{
return Err("--grid must be an integer 0..=4".into());
}
let ml_opts = hartree::multilevel::MultiLevelOptions {
compute_frequencies: do_freq,
symmetry_number,
qrrho_w0_cm1: qrrho_w0,
grid_override: grid,
all_electron,
ri_mp2,
solvent_eps,
smd: smd.clone(),
alpb: alpb.clone(),
gbsa: gbsa.clone(),
};
if let Some(engine) = &conformers {
if !rerank {
return Err(
"--conformers with a multi-level '//' spec needs --rerank (the CENSO-lite \
re-ranking hook); without it the spec is ambiguous for the screening step"
.into(),
);
}
return run_multilevel_conformers(
&molecule,
&spec,
&ml_opts,
engine.as_deref(),
alpb.clone(),
conformers_out.as_deref(),
);
}
let res = hartree::multilevel::run_multilevel(&molecule, &spec, &ml_opts)?;
report_multilevel(&molecule, &spec, &res);
if !no_method_warnings {
let mut warnings = res.sp.method_warnings.clone();
for w in &res.opt.method_warnings {
if !warnings.contains(w) {
warnings.push(w.clone());
}
}
print_method_warnings(&warnings);
}
return Ok(true);
}
if let Some(engine) = &conformers {
return run_conformers(
&molecule,
&method_str,
&basis,
engine.as_deref(),
alpb.clone(),
conformers_out.as_deref(),
);
}
if let Some(xtb_method) = hartree::ext::xtb::XtbMethod::from_keyword(&method_str) {
return run_xtb(&molecule, &method_str, xtb_method, do_opt, alpb.clone());
}
if fod_cube.is_some() && !fod {
return Err("--fod-cube requires --fod".into());
}
if fod {
if method_explicit {
if method_str != "tpss" {
eprintln!(
"warning: FOD analysis is calibrated at TPSS/def2-TZVP (T_el = 5000 K); \
using --method {method_str} — interpret N_FOD thresholds with care"
);
}
} else {
method_str = String::from("tpss");
}
if basis_explicit {
if !basis.eq_ignore_ascii_case("def2-tzvp") {
eprintln!(
"warning: FOD analysis is calibrated at TPSS/def2-TZVP; using --basis {basis}"
);
}
} else {
basis = String::from("def2-tzvp");
}
}
{
let gcp_is_dftc = gcp_keyword
.as_deref()
.is_some_and(|k| k.replace('-', "") == "dftc");
if method_str.replace('-', "").contains("dftc") || gcp_is_dftc {
return Err(
"DFT-C is not implemented: its parameter tables (the 1296 pairwise \
c_AB/alpha_AB/beta_AB values for H\u{2013}Kr plus the fitted atomic and \
gCP/def2-SVPD parameters) are published only in the AIP supplementary \
material of J. Chem. Phys. 146, 234105 (2017), which hartree cannot \
vendor. Obtain the SI zip in a browser (https://doi.org/10.1063/1.4986962) \
or email the corresponding author. The B97M-V/def2-SVPD NCI protocol \
runs without DFT-C: --method b97m-v --basis def2-svpd."
.into(),
);
}
}
let (base_method_str, want_disp) = if let Some(base) = method_str.strip_suffix("-d3") {
(base.to_string(), Some(false))
} else if let Some(base) = method_str.strip_suffix("-d4") {
(base.to_string(), Some(true))
} else {
(method_str.clone(), None)
};
let composite = hartree::composite::composite(&base_method_str);
if let Some(c) = composite {
if want_disp.is_some() {
return Err(format!(
"{} defines its own {} and short-range corrections; a -d3/-d4 suffix is not allowed",
c.keyword,
c.dispersion.label()
));
}
if basis_explicit && !basis.eq_ignore_ascii_case(c.basis) {
return Err(format!(
"{} implies the {} basis; --basis {basis} conflicts \
(omit --basis or pass {})",
c.keyword, c.basis_label, c.basis
));
}
basis = String::from(c.basis);
}
let method = if let Some(c) = composite {
Method::Dft(FunctionalSpec::parse(c.functional).map_err(|e| e.to_string())?)
} else {
match base_method_str.as_str() {
"rhf" => Method::Rhf,
"uhf" => Method::Uhf,
"hf" => {
if multiplicity > 1 {
Method::Uhf
} else {
Method::Rhf
}
}
"rohf" => Method::Rohf,
"mp2" => Method::Mp2,
"ccsd" => Method::Ccsd,
"ccsd(t)" | "ccsdt" => Method::CcsdT,
other => match FunctionalSpec::parse(other) {
Ok(spec) => Method::Dft(spec),
Err(hartree::dft::DftError::UnknownFunctional(_)) => {
return Err(format!(
"unknown method {other:?} (expected hf, rhf, uhf, rohf, mp2, ccsd, ccsd(t), \
or a functional like svwn/pbe/blyp/b3lyp/pbe0/tpss/r2scan, \
optionally with a -d3 or -d4 suffix)"
));
}
Err(err) => return Err(err.to_string()),
},
}
};
let dispersion = if let Some(c) = composite {
Some(c.dispersion)
} else if let Some(d4) = want_disp {
let (suffix, model) = if d4 { ("-d4", "D4") } else { ("-d3", "D3(BJ)") };
if matches!(method, Method::Mp2 | Method::Ccsd | Method::CcsdT) {
return Err(format!(
"{suffix} is not supported for post-HF methods ({base_method_str}); \
the {model} correction applies to HF and DFT functionals"
));
}
let param_key = match &method {
Method::Rhf | Method::Uhf | Method::Rohf => "hf".to_string(),
Method::Dft(spec) => spec
.d4_param_set()
.map(str::to_string)
.unwrap_or_else(|| spec.name().to_string()),
_ => unreachable!(),
};
match Dispersion::for_method(d4, ¶m_key) {
Some(d) => Some(d),
None => {
return Err(format!(
"no {model} parametrization exists for {param_key} \
(supported: pbe, blyp, b3lyp, b3lyp5, pbe0, tpss, r2scan, hf; \n D4 additionally: b2plyp, revdsd-pbep86, pwpb95)"
));
}
}
} else {
None
};
let standalone_gcp = match &gcp_keyword {
Some(k) => {
if composite.is_some() {
return Err(format!(
"{} defines its own gCP treatment; --gcp is not allowed with a composite \
method",
base_method_str
));
}
Some(hartree::disp::GcpParams::by_keyword(k).ok_or_else(|| {
format!(
"unknown gCP parameter set {k:?} (available: r2scan-3c / def2-mtzvpp, \
dft/sv(p) / b3lyp-3c, pbeh-3c / def2-msvp; only the sets vendored from \
the reference implementation mctc-gcp are bundled)"
)
})?)
}
None => None,
};
if grid.is_some() && !matches!(method, Method::Dft(_)) {
return Err("--grid requires a functional method (e.g. --method pbe)".into());
}
let metadata_grid = match &method {
Method::Dft(spec) if grid.is_none() && spec.grid_sensitive() => {
let level = spec.recommended_grid_level();
eprintln!(
"note: {} is grid-sensitive; defaulting to grid level {level} \
(override with --grid)",
spec.name()
);
Some(level)
}
_ => None,
};
let grid_level = grid
.or(metadata_grid)
.unwrap_or(composite.map_or(3, |c| c.grid_level));
if grid_level > 4 {
return Err("--grid must be an integer 0..=4".into());
}
if do_opt && (solvent_eps.is_some() || smd.is_some()) {
eprintln!("note: solvated geometry optimization uses finite-difference gradients (slower)");
}
if do_opt && do_freq {
eprintln!(
"note: --freq is not computed in the same run as --opt for single-level jobs; \
re-run --freq at the printed optimized geometry (multi-level '//' specs and \
--protocol w1 combine the two)"
);
}
if !method_explicit && !fod && !no_method_warnings {
eprintln!(
"note: no --method given; defaulting to RHF (the historic default). For \
production work r2scan-3c is recommended — see `hartree --recommend general` \
(suppress this note with --no-method-warnings)"
);
}
let job = Job {
molecule: molecule.clone(),
basis: basis.clone(),
method,
options: JobOptions {
all_electron,
direct,
ri,
compute_properties: do_properties,
compute_frequencies: do_freq,
single_point_hessian: do_sph,
optimize_geometry: do_opt,
transition_state: do_ts,
symmetry_number,
qrrho_w0_cm1: qrrho_w0,
grid_level,
dispersion,
solvent_eps,
smd: smd.clone(),
alpb: alpb.clone(),
gbsa: gbsa.clone(),
cosmo_file: cosmo_file.as_ref().map(Into::into),
gcp: composite.and_then(|c| c.gcp).or(standalone_gcp),
srb: composite.and_then(|c| c.srb),
smearing: smear.map(|temperature_k| Smearing::Fermi { temperature_k }),
fod,
fod_cube: fod_cube.as_ref().map(Into::into),
ri_mp2,
cosx,
x2c,
ts_options,
n_threads: None,
mem_budget_bytes: None,
},
};
if let Some(na) = cp_na {
if na >= molecule.len() {
return Err(format!(
"--cp {na}: fragment A must leave at least one atom for fragment B \
(the complex has {} atoms)",
molecule.len()
));
}
let (qa, qb) = cp_charges.unwrap_or((0, 0));
let (ma, mb) = cp_mults.unwrap_or((1, 1));
let frags = hartree::CpFragments {
fragment_a: (0..na).collect(),
charge_a: qa,
multiplicity_a: ma,
charge_b: qb,
multiplicity_b: mb,
};
let cp = hartree::counterpoise(&job, &frags)?;
println!(
"hartree -- counterpoise (Boys-Bernardi) -- {} / {}",
method_str.to_ascii_uppercase(),
basis
);
println!(
"complex: {} atoms (charge {}, mult {}); fragment A: atoms 1-{na} \
(charge {qa}, mult {ma}); fragment B: atoms {}-{} (charge {qb}, mult {mb})",
molecule.len(),
molecule.charge,
molecule.multiplicity,
na + 1,
molecule.len()
);
println!();
println!(" E_AB^(AB) (complex) {:>20.12} Eh", cp.e_complex);
println!(
" E_A^(AB) (A + ghost-B) {:>20.12} Eh",
cp.e_a_in_dimer_basis
);
println!(
" E_B^(AB) (ghost-A + B) {:>20.12} Eh",
cp.e_b_in_dimer_basis
);
println!(" E_A^(A) (A alone) {:>20.12} Eh", cp.e_a);
println!(" E_B^(B) (B alone) {:>20.12} Eh", cp.e_b);
println!();
let to_kcal = 627.509474063;
println!(
" dE_int (uncorrected) {:>20.12} Eh ({:>10.4} kcal/mol)",
cp.interaction_uncorrected(),
cp.interaction_uncorrected() * to_kcal
);
println!(
" dE_int^CP (corrected) {:>20.12} Eh ({:>10.4} kcal/mol)",
cp.interaction_cp(),
cp.interaction_cp() * to_kcal
);
println!(
" delta_BSSE {:>20.12} Eh ({:>10.4} kcal/mol)",
cp.bsse(),
cp.bsse() * to_kcal
);
return Ok(true);
}
if cp_charges.is_some() || cp_mults.is_some() {
return Err("--cp-charges/--cp-mults require --cp <nA>".into());
}
let result = job.run()?;
let ecp_atoms = hartree::BasisSet::load(&basis)
.map(|set| hartree::ecp_summary(&molecule, &set))
.unwrap_or_default();
if let Some(opt) = &result.optimized_geometry {
report_optimization(&molecule, &basis, &method_str, &ecp_atoms, opt);
if !no_method_warnings {
print_method_warnings(&result.method_warnings);
}
return Ok(opt.converged);
}
if let Some(ts) = &result.transition_state {
report_transition_state(&molecule, &basis, &method_str, &ecp_atoms, ts);
if !no_method_warnings {
print_method_warnings(&result.method_warnings);
}
return Ok(ts.converged());
}
report(
&molecule,
&basis,
&method_str,
&ecp_atoms,
&result.scf,
result.ri.as_ref(),
result.cosx.as_ref(),
);
if x2c {
println!();
println!("X2C-1e scalar-relativistic Hamiltonian: active");
println!(
" spin-free one-electron X2C (c = 137.035999084 a.u., CODATA 2018); \
two-electron integrals nonrelativistic"
);
if do_properties {
println!(
" note: property operators are NOT picture-change corrected (X2C-1e caveat); \
dipole moments and populations use the nonrelativistic operators"
);
}
}
if let Some(dft) = &result.dft {
report_dft(dft, &result.scf);
}
if let (Some(e_nl), None) = (result.vv10_energy, &result.double_hybrid) {
println!(
" E_nl (VV10, non-SC) {:>20.12} Eh (post-SCF; included in total below)",
e_nl
);
println!(
" total energy + E_nl {:>20.12} Eh",
result.scf.energy + e_nl
);
}
if let Some((occ_a, occ_b)) = &result.scf.occupations {
let frac = |occ: &[f64]| occ.iter().filter(|&&f| f > 1e-6 && f < 1.0 - 1e-6).count();
let ts = result.scf.electronic_entropy.unwrap_or(0.0);
let free = result.scf.free_energy.unwrap_or(result.scf.energy);
let t_el = smear
.or(result.fod.as_ref().map(|f| f.temperature_k))
.expect("occupations imply a smearing temperature");
println!();
println!("Fermi smearing (T = {:.1} K):", t_el);
println!(
" fractionally occupied {:>10} α / {:>5} β (of {} orbitals each)",
frac(occ_a),
frac(occ_b),
occ_a.len()
);
println!(" T*S_el {:>20.12} Eh", ts);
println!(" free energy F=E-T*S {:>20.12} Eh", free);
}
if let Some(f) = &result.fod {
println!();
println!("FOD analysis (Grimme):");
println!(
" T_el {:>20.1} K (T = 5000 K + 20000 K * a_x)",
f.temperature_k
);
println!(
" N_FOD {:>20.6} ({:.6} alpha, {:.6} beta)",
f.n_fod, f.n_fod_alpha, f.n_fod_beta
);
if let Some(path) = &fod_cube {
println!(" rho_FOD cube written to {path}");
}
if f.n_fod >= 1.0 {
println!(
" WARNING: N_FOD = {:.3} >= 1.0 — strong static correlation \
(multireference character); single-reference results (DFT, MP2, CCSD(T)) \
may be unreliable for this system",
f.n_fod
);
} else if f.n_fod >= 0.5 {
println!(
" note: N_FOD = {:.3} is in the borderline 0.5–1.0 range — mild static \
correlation; inspect rho_FOD and frontier occupations",
f.n_fod
);
} else {
println!(
" N_FOD < 0.5: no significant static correlation indicated \
(single-reference methods are appropriate)"
);
}
}
if let Some(smd) = &result.smd {
const KCAL: f64 = 627.509_451;
println!();
println!("SMD solvation (Marenich, Cramer, Truhlar 2009):");
println!(" solvent {:>20}", smd.solvent);
println!(" dielectric ε {:>20.4}", smd.epsilon);
println!(" E (gas, this geom) {:>20.12} Eh", smd.e_gas);
println!(
" E (solution) {:>20.12} Eh (SCF total in solvent)",
smd.e_solution
);
println!(
" ΔG_EP {:>20.12} Eh ({:>9.3} kcal/mol)",
smd.g_ep,
smd.g_ep * KCAL
);
println!(
" G_CDS {:>20.12} Eh ({:>9.3} kcal/mol)",
smd.g_cds,
smd.g_cds * KCAL
);
println!(
" ΔG_solv {:>20.12} Eh ({:>9.3} kcal/mol)",
smd.dg_solv,
smd.dg_solv * KCAL
);
println!(
" standard state: 298 K, fixed 1 mol/L in gas and solution (no concentration term)"
);
} else if let Some(g) = &result.gbsa {
const KCAL: f64 = 627.509_451;
println!();
println!(
"{} solvation (xtb GFN2 parameters; post-SCF on Mulliken charges):",
g.model
);
println!(" solvent {:>20}", g.solvent);
println!(" dielectric ε {:>20.4}", g.epsilon);
println!(
" G_born {:>20.12} Eh ({:>9.3} kcal/mol)",
g.g_born,
g.g_born * KCAL
);
if g.g_hb != 0.0 {
println!(
" G_hb {:>20.12} Eh ({:>9.3} kcal/mol)",
g.g_hb,
g.g_hb * KCAL
);
}
println!(
" G_sasa {:>20.12} Eh ({:>9.3} kcal/mol)",
g.g_sasa,
g.g_sasa * KCAL
);
println!(
" G_shift {:>20.12} Eh ({:>9.3} kcal/mol)",
g.g_shift,
g.g_shift * KCAL
);
println!(
" ΔG_solv {:>20.12} Eh ({:>9.3} kcal/mol)",
g.g_solv,
g.g_solv * KCAL
);
println!(
" caveat: GFN2-fit parameters applied to ab-initio Mulliken charges (provenance \
documented)"
);
} else if let Some(e_solv) = result.scf.solvation_energy {
println!();
if let Some(path) = &cosmo_file {
println!("C-PCM (ideal conductor, ε = ∞) for COSMO-RS export:");
println!(
" E_diel {:>20.12} Eh (included in total energy)",
e_solv
);
println!(" .cosmo file written to {path}");
} else {
println!("C-PCM solvation (electrostatics only):");
if let Some(name) = &solvent {
println!(" solvent {:>20}", name);
}
println!(" dielectric ε {:>20.4}", solvent_eps.unwrap());
println!(
" E_solv {:>20.12} Eh (included in total energy)",
e_solv
);
}
}
if let Some(c) = composite {
let e_disp = result
.dispersion_energy
.expect("every composite carries a dispersion correction");
println!();
println!("{} composite:", c.keyword);
println!(
" {:<20}{:>20.12} Eh",
format!("E_SCF ({})", c.functional),
result.scf.energy
);
println!(
" {:<20}{:>20.12} Eh",
format!("E_{}", c.disp_label),
e_disp
);
if let (Some(e_gcp), Some(gcp)) = (result.gcp_energy, c.gcp) {
println!(
" {:<20}{:>20.12} Eh",
format!("E_gCP ({})", gcp.label),
e_gcp
);
}
if let Some(e_srb) = result.srb_energy {
println!(" E_SRB {:>20.12} Eh", e_srb);
}
println!(
" composite total {:>20.12} Eh",
result.scf.energy
+ e_disp
+ result.gcp_energy.unwrap_or(0.0)
+ result.srb_energy.unwrap_or(0.0)
);
} else if let (Some(e_disp), None) = (result.dispersion_energy, &result.double_hybrid) {
let label = dispersion
.expect("dispersion energy implies a model")
.label();
println!();
println!("dispersion {label:<9}{:>20.12} Eh", e_disp);
println!(
"total energy + disp {:>20.12} Eh",
result.scf.energy + e_disp
);
}
if composite.is_none()
&& let (Some(e_gcp), Some(p)) = (result.gcp_energy, standalone_gcp)
{
println!();
println!("gCP ({}) {:>20.12} Eh", p.label, e_gcp);
println!(
"total energy incl. corrections {:>20.12} Eh",
result.best_energy()
);
}
if let Some(dh) = &result.double_hybrid {
println!();
println!(
"double hybrid {} (PT2 on {} orbitals):",
dh.functional_name, dh.scf_functional_name
);
println!(" E_SCF (DH, no PT2) {:>20.12} Eh", dh.e_scf);
println!(
" E_PT2 os {:>20.12} Eh (c_os = {:.5} x {:.12})",
dh.c_os * dh.e_os,
dh.c_os,
dh.e_os
);
println!(
" E_PT2 ss {:>20.12} Eh (c_ss = {:.5} x {:.12})",
dh.c_ss * dh.e_ss,
dh.c_ss,
dh.e_ss
);
match &dh.pt2_aux_basis {
Some(aux) => println!(" PT2 backend {:>20}", format!("RI-MP2 ({aux})")),
None => println!(" PT2 backend {:>20}", "conventional MP2"),
}
println!(" frozen core {:>20} orbitals", dh.n_frozen);
if let Some(e_nl) = result.vv10_energy {
println!(" E_nl (VV10 x {:.5}){:>20.12} Eh", dh.vv10_scale, e_nl);
}
if let Some(e_disp) = result.dispersion_energy {
println!(" E_disp (D4) {:>20.12} Eh", e_disp);
}
println!(" total energy {:>20.12} Eh", result.best_energy());
}
if let Some(post) = &result.post_hf {
match post {
PostHfResult::Mp2 {
result: r,
n_frozen,
} => report_mp2(*n_frozen, r, result.scf.reference),
PostHfResult::RiMp2 {
result: r,
n_frozen,
aux_basis,
} => report_ri_mp2(*n_frozen, r, aux_basis, result.scf.reference),
PostHfResult::Ccsd {
result: r,
n_frozen,
} => report_ccsd(*n_frozen, r),
PostHfResult::CcsdT {
result: r,
n_frozen,
} => report_ccsdt(*n_frozen, r),
}
}
if let Some(props) = &result.properties {
report_properties(&molecule, props.dipole_au, &props.population);
}
if let Some(freq) = &result.frequencies {
if freq.is_sph {
println!();
println!(
"NOTE: single-point Hessian (SPH) — geometry taken as-is, gradient direction \
projected out (Spicher & Grimme 2021). SPH frequencies are APPROXIMATE."
);
}
report_frequencies(&freq.frequencies);
report_thermo(&freq.thermochemistry);
}
if !no_method_warnings {
print_method_warnings(&result.method_warnings);
}
Ok(result.converged())
}
fn print_method_warnings(warnings: &[String]) {
if warnings.is_empty() {
return;
}
println!();
println!("method-quality assessment (informational; suppress with --no-method-warnings):");
for w in warnings {
println!(" {w}");
}
}
fn report_recommendation(task: &str) -> Result<bool, String> {
let rec = hartree::guardrails::recommend(task).ok_or_else(|| {
format!(
"unknown --recommend task {task:?} (available: {})",
hartree::guardrails::recommendation_tasks().join(", ")
)
})?;
println!("hartree -- recommended level of theory: {}", rec.task);
println!();
println!(" level: {}", rec.level);
println!(" rationale: {}", rec.rationale);
println!();
println!(" run:");
for inv in rec.invocation {
println!(" {inv}");
}
if !rec.notes.is_empty() {
println!();
for note in rec.notes {
println!(" note: {note}");
}
}
Ok(true)
}
fn parse_pair<T: std::str::FromStr>(value: &str, flag: &str) -> Result<(T, T), String> {
let mut parts = value.split(',');
let err = || format!("{flag} expects two comma-separated values, e.g. {flag} 0,1");
let a = parts.next().ok_or_else(err)?.trim();
let b = parts.next().ok_or_else(err)?.trim();
if parts.next().is_some() {
return Err(err());
}
Ok((a.parse().map_err(|_| err())?, b.parse().map_err(|_| err())?))
}
pub(crate) fn take(args: &[String], i: &mut usize, flag: &str) -> Result<String, String> {
*i += 1;
args.get(*i)
.cloned()
.ok_or_else(|| format!("{flag} needs a value"))
}
const GAP_WARN_HF: f64 = 0.1;
const GAP_WARN_KS: f64 = 0.05;
const T1_WARN: f64 = 0.02;
fn report(
molecule: &Molecule,
basis: &str,
method: &str,
ecp_atoms: &[(String, u32, u32)],
result: &ScfResult,
ri: Option<&hartree::RiDiagnostics>,
cosx: Option<&hartree::CosxDiagnostics>,
) {
println!("hartree -- {} / {}", method.to_ascii_uppercase(), basis);
println!(
"atoms: {} charge: {} multiplicity: {} electrons: {} ({} alpha, {} beta) basis fns: {}",
molecule.len(),
molecule.charge,
molecule.multiplicity,
result.n_alpha + result.n_beta,
result.n_alpha,
result.n_beta,
result.n_basis,
);
print_ecp_line(ecp_atoms);
if let Some(ri) = ri {
println!(
"RI-JK density fitting: aux basis {} ({} aux fns)",
ri.aux_basis, ri.naux
);
}
if let Some(cosx) = cosx {
let rs_note = match cosx.rs_omega {
Some(omega) => format!(", RS Coulomb + erf(omega={omega}) kernels"),
None => String::new(),
};
println!(
"COSX semi-numerical exchange: grid {} ({} points{}{})",
cosx.grid,
cosx.n_points,
if cosx.overlap_fitted {
", overlap-fitted"
} else {
""
},
rs_note
);
}
println!("nuclear repulsion: {:.10} Eh", result.nuclear_repulsion);
println!();
println!("{:>5} {:>22} {:>11}", "iter", "energy / Eh", "|error|");
for step in &result.history {
println!(
"{:>5} {:>22.12} {:>11.2e}",
step.iteration, step.energy, step.error_norm
);
}
println!();
if result.converged {
println!("converged in {} iterations", result.iterations);
} else {
println!("NOT CONVERGED after {} iterations", result.iterations);
}
println!();
println!("total energy {:>20.12} Eh", result.energy);
println!("electronic {:>20.12} Eh", result.electronic_energy);
println!("nuclear repulsion {:>20.12} Eh", result.nuclear_repulsion);
if result.reference != Reference::Rhf {
println!("<S^2> {:>20.6}", result.spin_squared);
}
let (gap_alpha, gap_beta) = result.homo_lumo_gap();
let is_ks = result.xc_energy.is_some();
let threshold = if is_ks { GAP_WARN_KS } else { GAP_WARN_HF };
if result.reference == Reference::Uhf {
if let Some(g) = gap_alpha {
println!("HOMO-LUMO gap (α) {:>20.6} Eh", g);
}
if let Some(g) = gap_beta {
println!("HOMO-LUMO gap (β) {:>20.6} Eh", g);
}
} else if let Some(g) = gap_alpha {
println!("HOMO-LUMO gap {:>20.6} Eh", g);
}
for (label, gap) in [("α", gap_alpha), ("β", gap_beta)] {
if let Some(g) = gap {
if g < threshold {
println!(
"warning: small HOMO-LUMO gap ({g:.4} Eh < {threshold} Eh{}) — possible \
multi-reference character or SCF instability",
if result.reference == Reference::Uhf {
format!(", {label} channel")
} else {
String::new()
}
);
if result.reference != Reference::Uhf {
break; }
}
}
}
}
fn report_dft(dft: &DftDiagnostics, scf: &ScfResult) {
let ks = if scf.reference == Reference::Uhf {
"UKS"
} else {
"RKS"
};
println!();
println!("Kohn-Sham DFT: {} ({})", dft.functional_name, ks);
println!(
" grid level {} ({} points)",
dft.grid_level, dft.n_grid_points
);
if let Some(exc) = scf.xc_energy {
println!(" E_xc {:>20.12} Eh", exc);
}
if let Some(ne) = scf.n_elec_grid {
println!(
" grid electrons ∫ρ {:>20.10} (N = {})",
ne,
scf.n_alpha + scf.n_beta
);
}
if dft.exx_fraction > 0.0 {
println!(" exact exchange c_x {:>20.4}", dft.exx_fraction);
}
}
fn report_mp2(n_frozen: usize, mp2: &hartree::cc::Mp2Result, reference: hartree::scf::Reference) {
let kind = if reference == hartree::scf::Reference::Uhf {
"UHF"
} else {
"RHF"
};
println!();
println!(
"{kind}-MP2 (frozen core: {} orbital{})",
n_frozen,
if n_frozen == 1 { "" } else { "s" }
);
println!(" opposite-spin {:>20.12} Eh", mp2.opposite_spin);
println!(" same-spin {:>20.12} Eh", mp2.same_spin);
println!(" correlation energy {:>20.12} Eh", mp2.correlation_energy);
println!(" total energy {:>20.12} Eh", mp2.total_energy);
}
fn report_ri_mp2(
n_frozen: usize,
mp2: &hartree::cc::RiMp2Result,
aux_basis: &str,
reference: hartree::scf::Reference,
) {
let kind = if reference == hartree::scf::Reference::Uhf {
"UHF"
} else {
"RHF"
};
println!();
println!(
"{kind}-RI-MP2 (frozen core: {} orbital{})",
n_frozen,
if n_frozen == 1 { "" } else { "s" }
);
println!(
" aux basis (MP2-fit) {:>20} ({} aux fns)",
aux_basis, mp2.naux
);
println!(" opposite-spin {:>20.12} Eh", mp2.opposite_spin);
println!(" same-spin {:>20.12} Eh", mp2.same_spin);
println!(" correlation energy {:>20.12} Eh", mp2.correlation_energy);
println!(" total energy {:>20.12} Eh", mp2.total_energy);
}
fn report_ccsd(n_frozen: usize, cc: &hartree::cc::CcsdResult) {
println!();
println!(
"RHF-CCSD (spin-adapted; frozen core: {} orbital{})",
n_frozen,
if n_frozen == 1 { "" } else { "s" }
);
if cc.converged {
println!(" converged in {} iterations", cc.iterations);
} else {
println!(" NOT CONVERGED after {} iterations", cc.iterations);
}
println!(" MP2 (iter-0) check {:>20.12} Eh", cc.mp2_correlation);
println!(" correlation energy {:>20.12} Eh", cc.correlation_energy);
println!(" total energy {:>20.12} Eh", cc.total_energy);
println!(" T1 diagnostic {:>20.6}", cc.t1_diagnostic);
if cc.t1_diagnostic > T1_WARN {
println!(
"warning: T1 diagnostic {:.4} > {T1_WARN} — significant singles amplitudes; \
a single-reference CC treatment may be unreliable",
cc.t1_diagnostic
);
}
}
fn report_ccsdt(n_frozen: usize, r: &hartree::cc::CcsdTResult) {
report_ccsd(n_frozen, &r.ccsd);
println!();
println!(
"RHF-CCSD(T) (frozen core: {} orbital{})",
n_frozen,
if n_frozen == 1 { "" } else { "s" }
);
println!(" triples correction {:>20.12} Eh", r.triples_energy);
println!(" total energy {:>20.12} Eh", r.total_energy);
}
fn print_ecp_line(ecp_atoms: &[(String, u32, u32)]) {
if ecp_atoms.is_empty() {
return;
}
let list = ecp_atoms
.iter()
.map(|(sym, z, n_core)| format!("{sym} (Z={z}, {n_core} core electrons replaced)"))
.collect::<Vec<_>>()
.join(", ");
println!("effective core potentials (def2-ECP): {list}");
}
fn report_optimization(
molecule: &Molecule,
basis: &str,
method: &str,
ecp_atoms: &[(String, u32, u32)],
result: &OptResult,
) {
println!(
"hartree -- {} / {} [geometry optimization]",
method.to_ascii_uppercase(),
basis
);
println!(
"atoms: {} charge: {} multiplicity: {}",
molecule.len(),
molecule.charge,
molecule.multiplicity,
);
print_ecp_line(ecp_atoms);
println!();
println!(
"{:>5} {:>22} {:>11} {:>11} {:>11}",
"step", "energy / Eh", "max force", "rms force", "max disp"
);
for step in &result.history {
println!(
"{:>5} {:>22.12} {:>11.2e} {:>11.2e} {:>11.2e}",
step.iteration, step.energy, step.max_force, step.rms_force, step.max_disp
);
}
println!();
if result.converged {
println!("optimization converged in {} steps", result.iterations);
} else {
println!(
"optimization NOT CONVERGED after {} steps",
result.iterations
);
}
println!();
println!("optimized geometry (angstrom):");
for (atom, pos) in molecule.atoms.iter().zip(&result.positions) {
println!(
"{:<2} {:>14.8} {:>14.8} {:>14.8}",
atom.element.symbol(),
pos[0] / ANGSTROM_TO_BOHR,
pos[1] / ANGSTROM_TO_BOHR,
pos[2] / ANGSTROM_TO_BOHR,
);
}
println!();
println!("total energy {:>20.12} Eh", result.energy);
}
fn report_transition_state(
molecule: &Molecule,
basis: &str,
method: &str,
ecp_atoms: &[(String, u32, u32)],
result: &TsResult,
) {
println!(
"hartree -- {} / {} [transition-state search]",
method.to_ascii_uppercase(),
basis
);
println!(
"atoms: {} charge: {} multiplicity: {}",
molecule.len(),
molecule.charge,
molecule.multiplicity,
);
print_ecp_line(ecp_atoms);
println!();
println!(
"{:>5} {:>22} {:>11} {:>11} {:>11}",
"step", "energy / Eh", "max force", "rms force", "max disp"
);
for step in &result.history {
println!(
"{:>5} {:>22.12} {:>11.2e} {:>11.2e} {:>11.2e}",
step.iteration, step.energy, step.max_force, step.rms_force, step.max_disp
);
}
println!();
match result.status {
TsStatus::Converged => {
println!("transition state converged in {} steps", result.iterations)
}
TsStatus::NotConverged => println!(
"transition-state search NOT CONVERGED after {} steps",
result.iterations
),
TsStatus::WrongImaginaryModeCount => println!(
"geometry converged in {} steps, but it is NOT a first-order saddle \
(wrong number of imaginary modes)",
result.iterations
),
TsStatus::StoppedEarly => {
println!("transition-state search stopped early by an observer")
}
_ => println!("transition-state search finished with an unrecognized status"),
}
if let Some(v) = &result.verification {
let n_neg = v.negative_eigenvalues.len();
println!(
"saddle check: {} imaginary mode{}",
n_neg,
if n_neg == 1 { "" } else { "s" }
);
if let Some(freq) = v.imaginary_frequency_cm1 {
println!(" imaginary frequency {:>13.2} cm^-1", freq);
}
}
println!();
println!("transition-state geometry (angstrom):");
for (atom, pos) in molecule.atoms.iter().zip(&result.positions) {
println!(
"{:<2} {:>14.8} {:>14.8} {:>14.8}",
atom.element.symbol(),
pos[0] / ANGSTROM_TO_BOHR,
pos[1] / ANGSTROM_TO_BOHR,
pos[2] / ANGSTROM_TO_BOHR,
);
}
println!();
println!("total energy {:>20.12} Eh", result.energy);
if let Some(irc) = &result.irc {
println!();
println!("IRC endpoint confirmation (short downhill trace along the reaction mode):");
let df = irc.forward_energy - result.energy;
let dr = irc.reverse_energy - result.energy;
println!(
" forward endpoint energy {:>20.12} Eh ΔE = {:>9.2e} Eh",
irc.forward_energy, df
);
println!(
" reverse endpoint energy {:>20.12} Eh ΔE = {:>9.2e} Eh",
irc.reverse_energy, dr
);
if df <= 0.0 && dr <= 0.0 {
println!(" both endpoints relaxed below the saddle (ΔE <= 0).");
} else {
println!(
" note: an endpoint did not relax below the saddle (ΔE > 0); \
inspect the geometry."
);
}
}
}
fn report_properties(molecule: &Molecule, mu_au: [f64; 3], pop: &PopulationAnalysis) {
println!();
println!("one-electron properties:");
let mu_debye: [f64; 3] = mu_au.map(|x| x * AU_DIPOLE_TO_DEBYE);
let mu_mag = mu_debye.iter().map(|x| x * x).sum::<f64>().sqrt();
println!(
" dipole moment: x = {:>10.6} D y = {:>10.6} D z = {:>10.6} D",
mu_debye[0], mu_debye[1], mu_debye[2]
);
println!(" |μ| {:>10.6} Debye", mu_mag);
println!();
println!(" atom element Mulliken Löwdin");
for (i, atom) in molecule.atoms.iter().enumerate() {
println!(
" {:>4} {:>7} {:>10.6} {:>10.6}",
i + 1,
atom.element.symbol(),
pop.mulliken_charges[i],
pop.lowdin_charges[i],
);
}
println!();
println!(" Mayer bond orders:");
let n = molecule.len();
for i in 0..n {
for j in i + 1..n {
let b = pop.mayer_bond_orders[i][j];
if b > 0.1 {
println!(
" {:>2}({}) – {:>2}({}) {:>8.4}",
i + 1,
molecule.atoms[i].element.symbol(),
j + 1,
molecule.atoms[j].element.symbol(),
b
);
}
}
}
}
fn run_xtb(
molecule: &Molecule,
method_str: &str,
method: hartree::ext::xtb::XtbMethod,
do_opt: bool,
alpb: Option<String>,
) -> Result<bool, String> {
use hartree::ext::xtb::{XtbInput, XtbRun, run};
let run_kind = if do_opt { XtbRun::Opt } else { XtbRun::Energy };
let input = XtbInput::from_molecule(method, molecule, alpb);
let workdir = std::env::temp_dir().join(format!("hartree_xtb_{}", std::process::id()));
let result = run(molecule, &input, run_kind, &workdir).map_err(|e| e.to_string())?;
println!(
"hartree -- {} (external xtb)",
method_str.to_ascii_uppercase()
);
println!(
"atoms: {} charge: {} multiplicity: {}",
molecule.len(),
molecule.charge,
molecule.multiplicity
);
println!();
println!("total energy {:>20.12} Eh", result.energy);
if let Some(opt) = &result.optimized {
println!();
println!("optimized geometry (angstrom):");
for atom in &opt.atoms {
println!(
"{:<2} {:>14.8} {:>14.8} {:>14.8}",
atom.element.symbol(),
atom.position[0] / ANGSTROM_TO_BOHR,
atom.position[1] / ANGSTROM_TO_BOHR,
atom.position[2] / ANGSTROM_TO_BOHR,
);
}
}
Ok(true)
}
fn run_conformers(
molecule: &Molecule,
method_str: &str,
basis: &str,
engine: Option<&str>,
alpb: Option<String>,
out_path: Option<&str>,
) -> Result<bool, String> {
use hartree::ext::ensemble::Ensemble;
let ensemble: Ensemble = if engine == Some("crest") {
use hartree::ext::crest::{CrestInput, run};
let method = hartree::ext::xtb::XtbMethod::from_keyword(method_str)
.unwrap_or(hartree::ext::xtb::XtbMethod::Gfn2Xtb);
let input = CrestInput::from_molecule(method, molecule, alpb);
let workdir = std::env::temp_dir().join(format!("hartree_crest_{}", std::process::id()));
println!(
"hartree -- conformer ensemble (external CREST, {})",
method.keyword()
);
run(molecule, &input, &workdir).map_err(|e| e.to_string())?
} else {
use hartree::ext::confgen::{ConfGenOptions, generate_conformers};
let method = simple_method(method_str)?;
let opts = ConfGenOptions::default();
println!(
"hartree -- conformer ensemble (fallback torsion-driving generator, {} / {})",
method_str.to_ascii_uppercase(),
basis
);
let res = generate_conformers(molecule, &opts, |m| {
let job = Job {
molecule: m.clone(),
basis: basis.to_string(),
method: method.clone(),
options: JobOptions::default(),
};
job.run().ok().and_then(|r| {
if r.converged() {
Some(r.best_energy())
} else {
None
}
})
})
.map_err(|e| e.to_string())?;
println!(
"rotatable bonds: {} driven: {} candidates: {} screened: {}",
res.rotatable_bonds.len(),
res.driven_bonds.len(),
res.n_candidates,
res.n_screened
);
res.ensemble
};
if ensemble.is_empty() {
return Err("conformer search produced no valid conformers".into());
}
let rel_kcal = ensemble.relative_energies_kcal();
let weights = ensemble.boltzmann_weights(298.15);
println!();
println!(
"ensemble ({} conformers, Boltzmann weights at 298.15 K):",
ensemble.len()
);
println!(
" {:>4} {:>20} {:>12} {:>10}",
"rank", "energy / Eh", "ΔE / kcal", "weight"
);
for (i, c) in ensemble.conformers.iter().enumerate() {
println!(
" {:>4} {:>20.12} {:>12.4} {:>10.4}",
i + 1,
c.energy,
rel_kcal[i],
weights[i]
);
}
if let Some(path) = out_path {
let mut text = String::new();
for c in &ensemble.conformers {
text.push_str(
&hartree::ext::xyz::write_xyz(&c.molecule, &format!("{:.12}", c.energy))
.map_err(|e| e.to_string())?,
);
}
std::fs::write(path, text).map_err(|e| format!("writing {path}: {e}"))?;
println!();
println!("ensemble written to {path}");
}
Ok(true)
}
fn report_w1(molecule: &Molecule, opts: &hartree::w1::W1Options, res: &hartree::w1::W1Result) {
println!("hartree -- hartree-W1 composite thermochemistry protocol");
println!("(a W1-STYLE protocol from hartree's native methods -- NOT literal W1/W1-F12)");
println!(
"atoms: {} charge: {} multiplicity: {}",
molecule.len(),
molecule.charge,
molecule.multiplicity
);
println!();
println!("protocol stages:");
println!(" 1. geometry + frequencies {}", res.opt_label);
println!(
" 2. HF/CBS {} / {} (Karton-Martin 2-pt, E(L) = E_CBS + A(L+1)e^(-9*sqrt(L)))",
res.basis_small, res.basis_large
);
println!(
" 3. CCSD corr/CBS {} / {} (Halkier 2-pt n^-3), frozen core: {}",
res.basis_small, res.basis_large, res.n_frozen
);
println!(
" 4. (T) additivity {} (unextrapolated)",
res.basis_small
);
println!(" omitted (documented): core-valence correlation, scalar relativity, spin-orbit");
println!();
println!("optimized geometry (angstrom):");
for atom in &res.geometry.atoms {
println!(
"{:<2} {:>14.8} {:>14.8} {:>14.8}",
atom.element.symbol(),
atom.position[0] / ANGSTROM_TO_BOHR,
atom.position[1] / ANGSTROM_TO_BOHR,
atom.position[2] / ANGSTROM_TO_BOHR,
);
}
println!();
println!("hartree-W1 energy breakdown:");
println!(
" E_opt ({:<}) {:>20.12} Eh (at its optimized geometry)",
res.opt_label, res.e_opt
);
println!(
" E_HF ({:<}) {:>20.12} Eh",
res.basis_small, res.e_hf_small
);
println!(
" E_HF ({:<}) {:>20.12} Eh",
res.basis_large, res.e_hf_large
);
println!(
" E_HF/CBS {:>20.12} Eh",
res.e_hf_cbs
);
println!(
" E_corr(CCSD, {:<}) {:>20.12} Eh",
res.basis_small, res.e_ccsd_corr_small
);
println!(
" E_corr(CCSD, {:<}) {:>20.12} Eh",
res.basis_large, res.e_ccsd_corr_large
);
println!(
" E_corr(CCSD)/CBS {:>20.12} Eh",
res.e_ccsd_corr_cbs
);
println!(
" E_(T) ({:<}) {:>20.12} Eh",
res.basis_small, res.e_t_small
);
println!(
" E(hartree-W1) {:>20.12} Eh = E_HF/CBS + E_corr/CBS + E_(T)",
res.electronic_energy()
);
if let Some(t) = &res.thermo {
report_frequencies(&t.freq.frequencies);
report_thermo(&t.freq.thermochemistry);
println!();
println!(
"hartree-W1 thermochemistry (thermal corrections at {}, sigma = {}):",
res.opt_label, opts.symmetry_number
);
println!(" H_corr(low) = H-E {:>20.12} Eh", t.h_corr);
println!(" G_corr(low) = G-E {:>20.12} Eh (RRHO)", t.g_corr);
println!(" G_corr(low, mRRHO) {:>20.12} Eh", t.g_corr_qrrho);
println!(
" H(hartree-W1) {:>20.12} Eh = E + H_corr(low)",
t.enthalpy
);
println!(
" G(hartree-W1, RRHO) {:>20.12} Eh = E + G_corr(low)",
t.gibbs
);
println!(
" G(hartree-W1, mRRHO) {:>20.12} Eh (recommended)",
t.gibbs_qrrho
);
}
}
fn report_multilevel(
molecule: &Molecule,
spec: &hartree::multilevel::MultiLevelSpec,
res: &hartree::multilevel::MultiLevelResult,
) {
println!("hartree -- multi-level: {}", spec.label());
println!(
"atoms: {} charge: {} multiplicity: {}",
molecule.len(),
molecule.charge,
molecule.multiplicity
);
println!(" SP level (energy) {}", spec.sp.label);
println!(" OPT level (geometry) {}", spec.opt.label);
println!();
let opt = res
.opt
.optimized_geometry
.as_ref()
.expect("multi-level result carries the optimization");
println!(
"geometry: optimized at {} ({} steps, converged)",
spec.opt.label, opt.iterations
);
println!("optimized geometry (angstrom):");
for atom in &res.geometry.atoms {
println!(
"{:<2} {:>14.8} {:>14.8} {:>14.8}",
atom.element.symbol(),
atom.position[0] / ANGSTROM_TO_BOHR,
atom.position[1] / ANGSTROM_TO_BOHR,
atom.position[2] / ANGSTROM_TO_BOHR,
);
}
println!();
println!("multi-level energy breakdown:");
println!(
" E_low ({:<}) {:>20.12} Eh (at its optimized geometry)",
spec.opt.label, res.e_low
);
println!(
" E(SCF) ({:<}) {:>20.12} Eh",
spec.sp.label, res.sp.scf.energy
);
for (label, value) in [
("E_disp", res.sp.dispersion_energy),
("E_gCP", res.sp.gcp_energy),
("E_SRB", res.sp.srb_energy),
("E_nl (VV10)", res.sp.vv10_energy),
] {
if let Some(v) = value {
println!(" {label:<21} {v:>20.12} Eh (SP level)");
}
}
if let Some(post) = &res.sp.post_hf {
println!(
" E_corr (post-HF) {:>20.12} Eh (SP level)",
post.total_energy() - res.sp.scf.energy
);
}
if let Some(dh) = &res.sp.double_hybrid {
println!(
" E_PT2 ({:<}) {:>20.12} Eh (SP level)",
dh.functional_name,
dh.pt2_energy()
);
}
println!(
" E_high//low {:>20.12} Eh (final electronic energy)",
res.e_high
);
if let Some(t) = &res.thermo {
report_frequencies(&t.freq.frequencies);
report_thermo(&t.freq.thermochemistry);
println!();
println!(
"composite thermochemistry (thermal corrections at {}, electronic energy at {}):",
spec.opt.label, spec.sp.label
);
println!(" E_high//low {:>20.12} Eh", res.e_high);
println!(" H_corr(low) = H-E {:>20.12} Eh", t.h_corr);
println!(" G_corr(low) = G-E {:>20.12} Eh (RRHO)", t.g_corr);
println!(" G_corr(low, mRRHO) {:>20.12} Eh", t.g_corr_qrrho);
println!(
" H(composite) {:>20.12} Eh = E_high + H_corr(low)",
t.enthalpy
);
println!(
" G(composite, RRHO) {:>20.12} Eh = E_high + G_corr(low)",
t.gibbs
);
println!(
" G(composite, mRRHO) {:>20.12} Eh (recommended)",
t.gibbs_qrrho
);
}
}
fn run_multilevel_conformers(
molecule: &Molecule,
spec: &hartree::multilevel::MultiLevelSpec,
ml_opts: &hartree::multilevel::MultiLevelOptions,
engine: Option<&str>,
alpb: Option<String>,
out_path: Option<&str>,
) -> Result<bool, String> {
use hartree::ext::ensemble::Ensemble;
use hartree::multilevel::ENSEMBLE_RERANK_CAP;
let ensemble: Ensemble = if engine == Some("crest") {
use hartree::ext::crest::{CrestInput, run};
let method = hartree::ext::xtb::XtbMethod::Gfn2Xtb;
let input = CrestInput::from_molecule(method, molecule, alpb);
let workdir = std::env::temp_dir().join(format!("hartree_crest_{}", std::process::id()));
println!(
"hartree -- multi-level conformer re-ranking (external CREST, {})",
method.keyword()
);
run(molecule, &input, &workdir).map_err(|e| e.to_string())?
} else {
use hartree::ext::confgen::{ConfGenOptions, generate_conformers};
println!(
"hartree -- multi-level conformer re-ranking (fallback generator, \
screened at RHF/STO-3G)"
);
let res = generate_conformers(molecule, &ConfGenOptions::default(), |m| {
let job = Job {
molecule: m.clone(),
basis: "sto-3g".into(),
method: Method::Rhf,
options: JobOptions::default(),
};
job.run().ok().and_then(|r| {
if r.converged() {
Some(r.best_energy())
} else {
None
}
})
})
.map_err(|e| e.to_string())?;
res.ensemble
};
if ensemble.is_empty() {
return Err("conformer search produced no valid conformers".into());
}
if ensemble.len() > ENSEMBLE_RERANK_CAP {
println!(
"note: re-ranking the {ENSEMBLE_RERANK_CAP} lowest of {} conformers \
(the documented cap keeps the hook cheap)",
ensemble.len()
);
}
let ranked =
hartree::multilevel::rerank_ensemble(&ensemble, spec, ml_opts, ENSEMBLE_RERANK_CAP)?;
println!();
println!(
"re-ranked ensemble ({} conformers, {}; Boltzmann weights at 298.15 K \
from the composite energies):",
ranked.len(),
spec.label()
);
println!(
" {:>4} {:>20} {:>20} {:>12} {:>10}",
"rank", "E_low / Eh", "E_high//low / Eh", "ΔE / kcal", "weight"
);
const KCAL: f64 = 627.509_474_063;
let e0 = ranked.first().map(|r| r.e_high).unwrap_or(0.0);
for (i, r) in ranked.iter().enumerate() {
println!(
" {:>4} {:>20.12} {:>20.12} {:>12.4} {:>10.4}",
i + 1,
r.e_low,
r.e_high,
(r.e_high - e0) * KCAL,
r.weight
);
}
if let Some(path) = out_path {
let mut text = String::new();
for r in &ranked {
text.push_str(
&hartree::ext::xyz::write_xyz(&r.molecule, &format!("{:.12}", r.e_high))
.map_err(|e| e.to_string())?,
);
}
std::fs::write(path, text).map_err(|e| format!("writing {path}: {e}"))?;
println!();
println!("re-ranked ensemble written to {path}");
}
Ok(true)
}
fn simple_method(method_str: &str) -> Result<Method, String> {
match method_str {
"rhf" | "hf" => Ok(Method::Rhf),
"uhf" => Ok(Method::Uhf),
"rohf" => Ok(Method::Rohf),
other => FunctionalSpec::parse(other).map(Method::Dft).map_err(|_| {
format!("unsupported conformer-screen method {other:?} (use hf or a functional)")
}),
}
}
fn report_frequencies(freq: &FrequencyResult) {
println!();
println!("harmonic vibrational frequencies (numerical Hessian of the gradient):");
for (i, &f) in freq.frequencies_cm1.iter().enumerate() {
let tag = if f < -1.0 {
"i"
} else if f.abs() < 10.0 {
"(trans/rot)"
} else {
""
};
let f_display = if f < 0.0 {
format!("{:.2}i", (-f))
} else {
format!("{:.2}", f)
};
println!(" {:>3}: {:>12} cm⁻¹ {}", i + 1, f_display, tag);
}
println!();
if freq.n_imaginary > 0 {
println!(
" WARNING: {} imaginary mode(s) — not at a minimum",
freq.n_imaginary
);
} else {
let n_vib = freq.frequencies_cm1.iter().filter(|&&f| f >= 10.0).count();
println!(" {n_vib} real vibrational modes (0 imaginary)");
}
}
fn report_thermo(thermo: &ThermoResult) {
println!();
println!(
"RRHO thermochemistry at {:.2} K (σ = {}, linear = {}):",
thermo.temperature, thermo.symmetry_number, thermo.is_linear
);
println!(
" zero-point energy (ZPE) {:>20.12} Eh",
thermo.zpe
);
println!(
" thermal correction to E {:>20.12} Eh",
thermo.thermal_energy_corr
);
println!(
" thermal correction to H {:>20.12} Eh",
thermo.enthalpy_corr
);
println!(
" total enthalpy H({:.0} K) {:>20.12} Eh",
thermo.temperature, thermo.enthalpy
);
println!(
" entropy S {:>20.12} Eh/K",
thermo.entropy
);
println!(
" T·S {:>20.12} Eh",
thermo.temperature * thermo.entropy
);
println!(
" Gibbs free energy G({:.0} K) {:>20.12} Eh",
thermo.temperature, thermo.gibbs
);
println!();
println!("quasi-RRHO (mRRHO) thermochemistry, Grimme Chem. Eur. J. 18, 9955 (2012),");
println!(
"w0 = {:.1} cm⁻¹ (entropy-only convention: H is the harmonic RRHO enthalpy):",
thermo.qrrho_w0_cm1
);
println!(
" entropy S(mRRHO) {:>20.12} Eh/K",
thermo.entropy_qrrho
);
println!(
" T·S(mRRHO) {:>20.12} Eh",
thermo.temperature * thermo.entropy_qrrho
);
println!(
" Gibbs free energy G(mRRHO) {:>20.12} Eh (recommended)",
thermo.gibbs_qrrho
);
}
fn print_usage() {
println!(
"hartree -- pure-Rust quantum chemistry\n\n\
USAGE:\n hartree <molecule.xyz> [options]\n\n\
OPTIONS:\n\
\x20 --basis <name> basis set [default: sto-3g]. Bundled: sto-3g, 6-31g;\n\
\x20 6-311g, 6-311g(d,p), 6-311+g(d,p), 6-311++g(d,p)\n\
\x20 (omits He); cc-pvdz, cc-pvtz, cc-pvqz; aug-cc-pvtz\n\
\x20 (all H–Ar); the Karlsruhe family, all-electron H–Kr:\n\
\x20 def2-svp, def2-svpd, def2-tzvp, def2-tzvpp, def2-tzvpd,\n\
\x20 def2-tzvppd, def2-qzvp, def2-qzvpp, def2-mtzvpp,\n\
\x20 def2-mtzvp (mTZVP), def2-msvp,\n\
\x20 ma-def2-svp, ma-def2-tzvp (minimally augmented);\n\
\x20 def2-svp and def2-tzvp also cover Ag, Sn, I, Au via\n\
\x20 the def2-ECP effective core potentials (auto-attached)\n\
\x20 --method <name> hf | rhf | uhf | rohf | mp2 | ccsd | ccsd(t) | <functional>\n\
\x20 (hf auto-selects rhf/uhf by multiplicity)\n\
\x20 [default: rhf] (mp2 auto-selects an RHF or UHF\n\
\x20 reference by multiplicity; ccsd/ccsd(t) are\n\
\x20 closed-shell RHF single points; a functional name — svwn, pbe, blyp,\n\
\x20 b3lyp, pbe0, tpss, r2scan, m06-2x, pw6b95, b97m-v,\n\
\x20 wb97x-v, wb97m-v, … — runs Kohn–Sham DFT, auto RKS/UKS.\n\
\x20 Double hybrids b2plyp, revdsd-pbep86, pwpb95,
\n \x20 wb97m(2)/wb97m2: closed-shell in-core single points
\n \x20 (E = E_SCF + c_os*E_os + c_ss*E_ss, frozen-core PT2
\n \x20 on the converged KS orbitals; wb97m(2) is evaluated
\n \x20 non-SC on wb97m-v orbitals, VV10 scaled by 0.65904).
\n \x20 Range-separated (wb97x-v/wb97m-v) and -V (VV10)\n\
\x20 functionals: conventional in-core single points, or\n\
\x20 --ri with --cosx for the RS hybrids (RI-J Coulomb +\n\
\x20 semi-numerical K; --ri alone and --direct are\n\
\x20 rejected for RS; no --opt/--freq); E_nl is evaluated\n\
\x20 non-self-consistently after the SCF; grid-sensitive\n\
\x20 functionals default to grid level 4)\n\
\x20 A -d3 suffix (pbe-d3, b3lyp-d3, hf-d3, …) adds the\n\
\x20 D3(BJ) dispersion correction; a -d4 suffix (pbe-d4,\n\
\x20 r2scan-d4, …) adds DFT-D4 (EEQ charges + ATM term)\n\
\x20 instead (neither for svwn or post-HF, not both at once)\n\
\x20 Composite (\"3c\") methods bundle functional + basis +\n\
\x20 corrections under one keyword; the basis is implied\n\
\x20 and -d3/-d4 suffixes are rejected (each defines its\n\
\x20 own corrections):\n\
\x20 r2scan-3c: r2scan/def2-mTZVPP + D4 (method-specific\n\
\x20 parameters) + gCP, grid defaults to 4\n\
\x20 b3lyp-3c: b3lyp5/def2-mSVP + D3(BJ)-ATM + gCP\n\
\x20 (DFT/SV(P) parameters)\n\
\x20 b97-3c: refitted B97 GGA/mTZVP + D3(BJ)-ATM + SRB\n\
\x20 (short-range bond-length correction; no gCP)\n\
\x20 pbeh-3c: modified PBE hybrid (42% EXX)/def2-mSVP\n\
\x20 + D3(BJ)-ATM + damped gCP\n\
\x20 Multi-level (\"//\") specs: --method \"SP // OPT\" runs\n\
\x20 the multi-level workflow — optimize at the OPT (low)\n\
\x20 level, final single point at the SP (high) level on\n\
\x20 that geometry (ORCA notation), e.g.\n\
\x20 --method \"wb97m-v/def2-tzvpp // r2scan-3c\". Each\n\
\x20 side is method[/basis]: composites carry their own\n\
\x20 basis (an explicit one is rejected); plain methods\n\
\x20 require one; -d3/-d4 suffixes apply per level. With\n\
\x20 --freq the frequencies/thermal corrections run at the\n\
\x20 OPT level and G = E_high + (G_low - E_low) is\n\
\x20 reported (composite free energy). Solvation flags\n\
\x20 apply to every stage by that stage's own rules;\n\
\x20 --opt/--basis and the single-point-only backends\n\
\x20 (--direct/--ri/--cosx/--x2c/--fod/--smear/--cp/\n\
\x20 --props/--sph/--gcp) are rejected with the spec\n\
\x20 --protocol <name> composite thermochemistry protocol. w1 (hartree-W1):\n\
\x20 a W1-STYLE staged workflow (NOT literal W1/W1-F12):\n\
\x20 B3LYP/cc-pVTZ geometry +\n\
\x20 frequencies/RRHO; HF/CBS from cc-pVTZ/cc-pVQZ\n\
\x20 (2-pt Karton-Martin); frozen-core CCSD-corr/CBS from\n\
\x20 the same pair (2-pt Halkier n^-3); (T) at cc-pVTZ\n\
\x20 (additivity). Closed-shell molecules only. Prints the\n\
\x20 stage breakdown and E/H/G. Core-valence correlation,\n\
\x20 scalar relativity, and spin-orbit are omitted\n\
\x20 (documented). Honours --symmetry-number, --qrrho-w0,\n\
\x20 --all-electron; rejects --method/--basis/--opt/--freq\n\
\x20 and the backend/solvation flags\n\
\x20 --rerank with --conformers and a multi-level \"//\" spec:\n\
\x20 CENSO-lite ensemble re-ranking — optimize each\n\
\x20 conformer at the OPT level, single-point at the SP\n\
\x20 level, Boltzmann weights from the composite energies\n\
\x20 (at most 6 conformers, the documented cap; the\n\
\x20 generation step is screened at RHF/STO-3G)\n\
\x20 --grid <0..4> DFT integration grid level [default: 3]\n\
\x20 (only valid with a functional method)\n\
\x20 --charge <int> net charge [default: 0]\n\
\x20 --spin <int> spin multiplicity 2S+1 [default: 1]\n\
\x20 --all-electron correlate all orbitals for mp2/ccsd/ccsd(t) (default:\n\
\x20 noble-gas frozen core)\n\
\x20 --opt optimize the geometry (analytic gradient for rhf/uhf,\n\
\x20 finite differences for rohf and DFT functionals)\n\
\x20 --ts locate a transition state (first-order saddle) by P-RFO\n\
\x20 eigenvector-following from the input geometry, which must\n\
\x20 be a guess near the saddle; verifies one imaginary mode\n\
\x20 and reports its frequency (rhf/uhf/rohf/DFT; mutually\n\
\x20 exclusive with --opt; no post-HF/RI/direct/COSX/X2C/\n\
\x20 smearing/implicit-solvent)\n\
\x20 --ts-irc after a converged TS, trace the IRC a short way downhill\n\
\x20 in both senses of the reaction mode and report the two\n\
\x20 endpoint energies (confirms the basins the saddle joins)\n\
\x20 --ts-algo <prfo|dimer> saddle-point algorithm [default: prfo]; both prfo\n\
\x20 (Hessian eigenvector-following) and dimer (Hessian-free,\n\
\x20 midpoint-gradient curvature estimate) are available\n\
\x20 --ts-max-iter <int> max saddle-search iterations [default: 300] (>= 1)\n\
\x20 --ts-trust <bohr> initial trust radius for the climbing step [default: 0.2]\n\
\x20 --ts-follow <int> P-RFO mode to follow uphill, 0 = softest [default: 0]\n\
\x20 --ts-recalc-hessian <int> recompute the FD Hessian every N accepted steps;\n\
\x20 0 = compute once, then Bofill update [default: 0]\n\
\x20 --ts-fd-step <bohr> finite-difference step for gradients/Hessian [default:\n\
\x20 5e-3]\n\
\x20 --ts-neg-tol <a.u.> eigenvalue cutoff for a negative (reaction) mode\n\
\x20 [default: 1e-4]\n\
\x20 (TS runs raise the SCF cap to 400 and add a 0.3 level\n\
\x20 shift for small-gap stability)\n\
\x20 (the --ts-* knobs take effect only with --ts)\n\
\x20 --direct integral-direct SCF: recompute ERIs each iteration\n\
\x20 instead of storing the nao^4 tensor (rhf/uhf/rohf\n\
\x20 single points; reaches larger systems, slower)\n\
\x20 --ri density-fitted (RI-JK) SCF: J/K from 3-center integrals\n\
\x20 over the def2-universal-jkfit auxiliary set (hf and all\n\
\x20 DFT functionals, single points; small fitting error\n\
\x20 ~1e-4 Eh; incompatible with --direct, --opt, post-HF,\n\
\x20 --properties, --freq)\n\
\x20 --ri-mp2 RI-MP2: density-fitted MP2 over the matching def2 /C\n\
\x20 (MP2-fit) auxiliary set (<basis>/c; bundled: def2-svp/c,\n\
\x20 def2-tzvp/c — errors if no /C partner exists, never\n\
\x20 falls back to jkfit). --method mp2 only (RHF or UHF\n\
\x20 reference by multiplicity); the SCF step\n\
\x20 keeps its backend (default in-core, or --ri). Same\n\
\x20 frozen-core convention and E_OS/E_SS split as\n\
\x20 conventional MP2 (fitting error ~1e-5 Eh at def2-svp)\n\
\x20 --cosx COSX semi-numerical exchange: K built on a coarse\n\
\x20 overlap-fitted grid (Neese et al. 2009 / Izsak & Neese\n\
\x20 2011); J keeps the configured path (in-core or --ri).\n\
\x20 HF, global-hybrid, and range-separated DFT single\n\
\x20 points (RS hybrids serve K_LR(omega) semi-numerically\n\
\x20 too, on either backend); rejected for post-HF,\n\
\x20 --direct, --opt, --freq, and --fod\n\
\x20 --x2c scalar-relativistic X2C-1e Hamiltonian: T+V replaced\n\
\x20 by the picture-changed spin-free exact-two-component\n\
\x20 one-electron Hamiltonian (2e integrals stay\n\
\x20 nonrelativistic). HF and DFT on any backend; rejected\n\
\x20 with ECP atoms (double-counts relativity), --opt,\n\
\x20 --freq, and post-HF; properties carry a picture-change\n\
\x20 caveat\n\
\x20 --solvent <name> C-PCM implicit solvation (electrostatics only) in a\n\
\x20 named solvent: water, acetonitrile, methanol, dmso,\n\
\x20 chloroform, toluene (SCF-level methods, incl. --ri and\n\
\x20 -d3; E_solv included in the total; --opt uses FD\n\
\x20 gradients; --freq in solvent is rejected)\n\
\x20 --eps <float> C-PCM with an explicit dielectric constant instead of\n\
\x20 a named solvent (mutually exclusive with --solvent)\n\
--smd <name> SMD universal solvation model (Marenich/Cramer/
\n Truhlar 2009): C-PCM electrostatics with SMD intrinsic
\n Coulomb radii + the CDS surface-tension term; reports
\n dG_EP, G_CDS, dG_solv (1 M -> 1 M standard state) vs a
\n gas-phase reference SCF. 20 bundled solvents (water,
\n methanol, ethanol, acetonitrile, dmso, acetone, thf,
\n dmf, toluene, benzene, chloroform, ...). Same scope as
\n C-PCM (SCF-level methods; --opt via FD; no --freq);
\n mutually exclusive with --solvent/--eps
\n --alpb <name> ALPB implicit solvation (Ehlert/Stahn/Spicher/
\n Grimme 2021), xtb GFN2 parameters. Post-SCF on
\n SCF Mulliken charges; reports G_born/G_hb/G_sasa/
\n G_shift and DeltaG_solv. Single points only.
\n GFN2-fit params on ab-initio charges (caveat).
\n --gbsa <name> GBSA implicit solvation (generalized Born + SASA),
\n xtb GFN2 parameters. Like --alpb but the Still
\n kernel and no Poisson-Boltzmann shape term.
\n --cosmo-file <path> run C-PCM at the ideal-conductor limit (eps=inf)
\n and write a .cosmo file (TURBOMOLE/COSMOtherm
\n format) for COSMO-RS. Mutually exclusive with the
\n other solvation options
\n \x20 --smear <K> Fermi-Dirac fractional-occupation smearing at the\n\
\x20 given electronic temperature (kelvin). Energy-only:\n\
\x20 rhf/uhf and DFT single points (no --opt, --freq, or\n\
\x20 post-HF). Reports occupations, T*S_el, and the free\n\
\x20 energy F = E - T*S_el\n\
\x20 --fod Grimme FOD multireference diagnostic (Angew. Chem. Int.\n\
\x20 Ed. 54, 12308 (2015)): Fermi-smeared SCF at\n\
\x20 T_el = 5000 K + 20000 K*a_x, reporting N_FOD (the\n\
\x20 integral of the fractional-occupation weighted density;\n\
\x20 N_FOD >= 1 signals strong static correlation). Implies\n\
\x20 --method tpss --basis def2-tzvp unless given explicitly\n\
\x20 (warned); --smear overrides the temperature. In-core\n\
\x20 single points only (no --direct/--ri/--opt/--freq or\n\
\x20 post-HF)\n\
\x20 --fod-cube <path> with --fod, export rho_FOD (alpha+beta) as a Gaussian\n\
\x20 cube file (bounding box + 4 bohr margin, 0.2 bohr\n\
\x20 spacing)\n\
\x20 --properties, --props compute one-electron properties: dipole (Debye),\n\
\x20 Mulliken/Löwdin charges, Mayer bond orders\n\
\x20 --freq compute harmonic frequencies (numerical Hessian of the\n\
\x20 gradient) and RRHO + quasi-RRHO (mRRHO, Grimme 2012)\n\
\x20 thermochemistry at 298.15 K. Works for rhf/uhf and DFT\n\
\x20 functionals incl. the 3c composites and ECP atoms;\n\
\x20 run at the optimized geometry for all-real modes\n\
\x20 --qrrho-w0 <cm-1> quasi-RRHO interpolation frequency w0 for the mRRHO\n\
\x20 entropy [default: 100]\n\
\x20 --sph single-point Hessian (Spicher & Grimme, JCTC 2021):\n\
\x20 with --freq, take the geometry as-is (skip the\n\
\x20 stationarity assumption) and project the gradient\n\
\x20 direction out of the Hessian. Reduces to ordinary\n\
\x20 frequencies at a true minimum; SPH frequencies are\n\
\x20 APPROXIMATE (gradient-projection variant, not xtb\n\
\x20 --bhess; see crate docs)\n\
\x20 --method gfn2-xtb run the external `xtb` binary (GFN2-xTB); --method\n\
\x20 --method gfn-ff gfn-ff runs GFN-FF. Detected on PATH or via\n\
\x20 HARTREE_XTB_PATH; honours --charge/--spin/--alpb/--opt.\n\
\x20 A clear 'xtb binary not found' error with install\n\
\x20 guidance is raised when absent\n\
\x20 --conformers [crest] conformer ensemble search. Bare --conformers uses the\n\
\x20 built-in RDKit-free torsion-driving fallback (screened\n\
\x20 with --method/--basis single points; default RHF/\n\
\x20 STO-3G); `--conformers crest` runs the external CREST\n\
\x20 (HARTREE_CREST_PATH or PATH). Prints relative energies\n\
\x20 and Boltzmann weights at 298.15 K\n\
\x20 --conformers-out <path> write the conformer ensemble as a multi-frame XYZ\n\
\x20 --symmetry-number <int> rotational symmetry number σ for RRHO entropy\n\
\x20 [default: 1; use 2 for H₂O, HF, CO; 12 for CH₄]\n\
\x20 --cp <nA> Boys-Bernardi counterpoise correction for a\n\
\x20 two-fragment complex: the first nA atoms of the XYZ\n\
\x20 (input order) are fragment A, the rest fragment B.\n\
\x20 Runs the 5 required single points (complex; each\n\
\x20 fragment with the partner's ghost basis; each\n\
\x20 fragment alone) and reports dE_int^CP, the\n\
\x20 uncorrected dE_int, and the BSSE estimate. Works\n\
\x20 with any single-point-capable method (no --opt/\n\
\x20 --freq/--props/--fod/solvent)\n\
\x20 --cp-charges <qA,qB> per-fragment net charges for --cp [default: 0,0];\n\
\x20 must sum to --charge\n\
\x20 --cp-mults <mA,mB> per-fragment spin multiplicities for --cp\n\
\x20 [default: 1,1]\n\
\x20 --gcp <set> standalone geometric counterpoise (gCP) correction\n\
\x20 added to the total energy (outside the 3c\n\
\x20 composites). Parameter sets (vendored from the\n\
\x20 reference implementation mctc-gcp): r2scan-3c\n\
\x20 (alias def2-mtzvpp), dft/sv(p) (alias b3lyp-3c),\n\
\x20 pbeh-3c (alias def2-msvp). Ghost atoms are excluded\n\
\x20 --recommend <task> print the recommended level of theory for a task and\n\
\x20 exit (no XYZ needed): general (r2scan-3c opt+freq),\n\
\x20 barriers / nci (wb97m-v/def2-tzvpp single point on a\n\
\x20 r2scan-3c geometry via the multi-level '//' workflow),\n\
\x20 thermochemistry (the multi-level composite free\n\
\x20 energy with --freq)\n\
\x20 --no-method-warnings suppress the method-quality assessment section\n\
\x20 (informational warnings: Pople/cc bases with DFT,\n\
\x20 missing dispersion where the functional metadata\n\
\x20 recommends one, minimal/unpolarized bases, pure-GGA\n\
\x20 barrier caveats, HF-without-correlation, coarse grids\n\
\x20 under grid-sensitive functionals) and the no-method\n\
\x20 recommendation pointer. The warnings stay attached to\n\
\x20 the library JobResult either way\n\
\n\
\x20 PERIODIC (solid-state GPW DFT; selected by --cell / --kpoints / --cutoff):\n\
\x20 --cell <spec> unit cell: 'cubic <a>', or 1/3/6/9 numbers (Å, deg):\n\
\x20 1=cubic, 3=a b c (orthorhombic), 6=a b c α β γ,\n\
\x20 9=three row vectors. Omit the spec (or '--cell file')\n\
\x20 to read a Lattice=\"…\" line from the extended-XYZ file\n\
\x20 --kpoints <n1 n2 n3> Monkhorst–Pack mesh (or 'gamma') [default: gamma]\n\
\x20 --cutoff <E> plane-wave density cutoff in hartree [default: 280]\n\
\x20 --basis <name> GTH basis [default: DZVP-GTH-PADE; also SZV-GTH]\n\
\x20 --xc <name> pade (GTH-PADE LDA) | lda (Slater+PW92) [default: pade]\n\
\x20 --pseudo <file> user GTH potential file (CP2K format) [default: bundled]\n\
\x20 --basis-file <file> user GTH basis file (.gbs) [default: bundled]\n\
\x20 --forces also compute analytic forces\n\
\x20 --stress also compute the analytic stress tensor\n\
\x20 --max-iter <n> --mixing <a> SCF controls\n\
\x20 (v1: GTH PPs, spin-restricted, insulators/semiconductors;\n\
\x20 bundled elements: H, C, O, Si)\n\
\x20 -h, --help show this help\n\n\
XYZ coordinates are in angstrom. Ghost atoms (counterpoise centers: basis\n\
functions only, no nuclear charge, no electrons) are written Gh(O) or with\n\
the @O shorthand in the element column; they are excluded from D3/D4/gCP/SRB\n\
and from geometry jobs (--opt/--freq), and rejected with --props/solvent."
);
}