shiozaki
commented
1 year ago Hi - it's in the commercial (QSimulate) version but it has not been ported to the open-source version (and at this point there is no immediate plan to do so). It is not a very complicated procedure. @bessvlai do you have a python script that you don't mind sharing??
PS: If you are skilled in C++, the relevant code is below. Just slot this in and recompile. if you are not familiar with C++ please don't, as I won't be able to do user support.
void Hess::compute_free_energy_() {
// evaluate free energies
const double R = boltzmann__ / au2joule__;
const double T = temperature_;
const double P = pressure_;
const double h = planck__;
// prepare rotational constants
const double mtot_amu = accumulate(geom_->atoms().begin(), geom_->atoms().end(), 0.0, [](double s, shared_ptr<const Atom> a) { return s+a->mass(); });
const double mtot = mtot_amu * amu2kilogram__;
const bool is_linear = geom_->is_linear();
VectorB rotconst(3);
{
array<double,3> center{{0.0, 0.0, 0.0}};
for (auto& i : geom_->atoms())
blas::ax_plus_y_n(i->mass() / mtot_amu, i->position().data(), 3, center.data());
XYZFile xyz(geom_->natom());
for (int i = 0; i != geom_->natom(); ++i) {
shared_ptr<const Atom> iatom = geom_->atoms(i);
xyz(0, i) += (iatom->position(0) - center[0]) * sqrt(iatom->mass());
xyz(1, i) += (iatom->position(1) - center[1]) * sqrt(iatom->mass());
xyz(2, i) += (iatom->position(2) - center[2]) * sqrt(iatom->mass());
}
Matrix inertia = xyz ^ xyz;
const double trace = inertia.trace();
inertia(0, 0) = trace - inertia(0, 0);
inertia(1, 1) = trace - inertia(1, 1);
inertia(2, 2) = trace - inertia(2, 2);
inertia.diagonalize(rotconst);
for (int i = 0; i != 3; ++i)
rotconst(i) = hbar__ / rotconst(i) / (4 * pi__ * amu2kilogram__ * pow(au2meter__, 2));
}
component_ = vector<string>{"electronic", "translational", "rotational (no symmetry)", "vibrational"};
const int ncomp = component_.size();
internal_.resize(ncomp);
entropy_.resize(ncomp);
heatv_.resize(ncomp);
// electronic contributions
internal_[0] = ref_->energy(0);
entropy_[0] = R * log(2*(ref_->noccA() - ref_->noccB())+1);
heatv_[0] = 0.0;
// translational contributions
const double kT = boltzmann__ * temperature_;
internal_[1] = 1.5 * R * T;
entropy_[1] = 2.5 * R + R * log(pow(2*pi__*kT*mtot / (h*h), 1.5) * kT / P);
heatv_[1] = 1.5 * R;
if (geom_->natom() != 1) {
assert(!is_linear || fabs(rotconst(1) - rotconst(2)) < 1.0e-8);
const double qrot = is_linear ? (boltzmann__ * T / (planck__ * rotconst(1)))
: pow(boltzmann__ * T / planck__, 1.5) * sqrt(pi__) / sqrt(rotconst(0)*rotconst(1)*rotconst(2));
// rotational contributions
internal_[2] = (is_linear ? 1.0 : 1.5) * R * T;
entropy_[2] = (is_linear ? 1.0 : 1.5) * R + log(qrot) * R;
heatv_[2] = (is_linear ? 1.0 : 1.5) * R;
// set vibrational contributions to NAN if there are negative modes larger than the cutoff
if (freq_[0] < -freq_cutoff_) {
internal_[3] = entropy_[3] = heatv_[3] = nan("");
} else {
// vibrational contributions
for (int i = (is_linear ? 5 : 6); i != freq_.size(); ++i) {
const double freq = fabs(freq_[i]) < freq_cutoff_ ? freq_cutoff_: freq_[i];
const double en = freq / au2wavenumber__ ;
const double tt = en / (R * T);
const double ee = exp(-tt);
internal_[3] += en * (0.5 + ee / (1.0 - ee));
entropy_[3] += R * (tt * ee / (1.0 - ee) - log(1.0 - ee));
heatv_[3] += R * ee * pow(tt / (1.0 - ee), 2);
}
}
}
bout << " * Free energy analysis" << endl;
auto print_one = [](const string component, const double G, const double U, const double S, const double Cv) {
bout << " == " << component << " ==" << endl;
bout << " Gibbs energy: " << setw(20) << setprecision(8) << G << endl;
bout << " Internal energy: " << setw(20) << setprecision(8) << U << endl;
bout << " Entropy: " << setw(20) << setprecision(8) << S << endl;
bout << " Heat capacity (V):" << setw(20) << setprecision(8) << Cv << endl << endl;
};
{
const double U = accumulate(internal_.begin(), internal_.end(), 0.0);
const double S = accumulate(entropy_.begin(), entropy_.end(), 0.0);
const double Cv = accumulate(heatv_.begin(), heatv_.end(), 0.0);
gibbs_tot_ = U + R * T - S * T;
entropy_tot_ = S;
print_one("total", gibbs_tot_, U, S, Cv);
}
for (int i = 0; i != ncomp; ++i)
print_one(component_[i], internal_[i] + (component_[i] == "translational" ? R*T : 0.0) - entropy_[i] * T, internal_[i], entropy_[i], heatv_[i]);
}
bessvlai
ABQTrap
Dear all,
Can BAGEL output thermochemistry e.g. enthalpy and gibbs free energy, after the hessian and frequency calculation? I couldn't find the thermochemistry information after the hessian calculation. Many thanks.