From ddf6512a3a06415bbfb25a0cc4ecdb6ce082630f Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Tue, 9 Mar 2021 14:33:55 +0000 Subject: [PATCH] Use constexpr for physical constants and move them into gmx namespace The physical constants that are part of the public API are converted from defines to constexprs and are moved into the gmx namespace. This means that many files are touched with tirvial changes. The benefit is that now the risk of symbol collision is reduced in cases where a library might depend on these constants. Additionally, in the same vein of library cleanup, the math/utility header is removed from units header and included only in the sources that actually need it. --- api/legacy/include/gromacs/math/units.h | 95 ++++++++++--------- src/gromacs/applied_forces/awh/awh.cpp | 4 +- src/gromacs/applied_forces/awh/dimparams.h | 4 +- src/gromacs/applied_forces/electricfield.cpp | 2 +- src/gromacs/ewald/ewald.cpp | 5 +- src/gromacs/ewald/long_range_correction.cpp | 8 +- src/gromacs/ewald/pme_gpu_internal.cpp | 2 +- src/gromacs/ewald/pme_gpu_types.h | 4 +- src/gromacs/ewald/pme_solve.cpp | 2 +- src/gromacs/fileio/pdbio.cpp | 14 +-- src/gromacs/fileio/tngio.cpp | 16 ++-- src/gromacs/gmxana/anadih.cpp | 10 +- src/gromacs/gmxana/angle_correction.h | 1 - src/gromacs/gmxana/gmx_angle.cpp | 14 +-- src/gromacs/gmxana/gmx_awh.cpp | 11 ++- src/gromacs/gmxana/gmx_bar.cpp | 8 +- src/gromacs/gmxana/gmx_bundle.cpp | 14 +-- src/gromacs/gmxana/gmx_chi.cpp | 14 +-- src/gromacs/gmxana/gmx_clustsize.cpp | 4 +- src/gromacs/gmxana/gmx_current.cpp | 13 +-- src/gromacs/gmxana/gmx_density.cpp | 4 +- src/gromacs/gmxana/gmx_densorder.cpp | 7 +- src/gromacs/gmxana/gmx_dipoles.cpp | 13 +-- src/gromacs/gmxana/gmx_dos.cpp | 30 +++--- src/gromacs/gmxana/gmx_enemat.cpp | 4 +- src/gromacs/gmxana/gmx_energy.cpp | 39 ++++---- src/gromacs/gmxana/gmx_hbond.cpp | 14 +-- src/gromacs/gmxana/gmx_nmeig.cpp | 40 ++++---- src/gromacs/gmxana/gmx_nmens.cpp | 4 +- src/gromacs/gmxana/gmx_nmtraj.cpp | 8 +- src/gromacs/gmxana/gmx_rama.cpp | 6 +- src/gromacs/gmxana/gmx_sham.cpp | 11 ++- src/gromacs/gmxana/gmx_sigeps.cpp | 6 +- src/gromacs/gmxana/gmx_spol.cpp | 4 +- src/gromacs/gmxana/gmx_tcaf.cpp | 9 +- src/gromacs/gmxana/gmx_traj.cpp | 4 +- src/gromacs/gmxana/gmx_velacc.cpp | 4 +- src/gromacs/gmxana/gmx_wham.cpp | 19 ++-- src/gromacs/gmxana/hxprops.cpp | 14 +-- src/gromacs/gmxana/pp2shift.cpp | 6 +- src/gromacs/gmxana/thermochemistry.cpp | 45 ++++----- src/gromacs/gmxpreprocess/calch.cpp | 8 +- src/gromacs/gmxpreprocess/convparm.cpp | 8 +- src/gromacs/gmxpreprocess/editconf.cpp | 16 ++-- .../gmxpreprocess/gen_maxwell_velocities.cpp | 12 +-- src/gromacs/gmxpreprocess/gen_vsite.cpp | 50 +++++----- src/gromacs/gmxpreprocess/genion.cpp | 4 +- src/gromacs/gmxpreprocess/grompp.cpp | 4 +- src/gromacs/gmxpreprocess/hizzie.cpp | 4 +- src/gromacs/gmxpreprocess/readir.cpp | 3 +- src/gromacs/gmxpreprocess/solvate.cpp | 4 +- src/gromacs/gmxpreprocess/topshake.cpp | 8 +- src/gromacs/gmxpreprocess/vsite_parm.cpp | 20 ++-- src/gromacs/gmxpreprocess/x2top.cpp | 4 +- src/gromacs/imd/imd.cpp | 8 +- src/gromacs/listed_forces/bonded.cpp | 48 +++++----- src/gromacs/listed_forces/restcbt.cpp | 6 +- src/gromacs/mdlib/calc_verletbuf.cpp | 6 +- src/gromacs/mdlib/calcmu.cpp | 10 +- src/gromacs/mdlib/coupling.cpp | 52 +++++----- src/gromacs/mdlib/dispersioncorrection.cpp | 2 +- src/gromacs/mdlib/energyoutput.cpp | 10 +- src/gromacs/mdlib/expanded.cpp | 17 ++-- src/gromacs/mdlib/forcerec.cpp | 2 +- src/gromacs/mdlib/lincs.cpp | 6 +- src/gromacs/mdlib/rf_util.cpp | 2 +- src/gromacs/mdlib/sim_util.cpp | 2 +- src/gromacs/mdlib/update.cpp | 8 +- src/gromacs/mdlib/updategroups.cpp | 7 +- src/gromacs/mdrun/replicaexchange.cpp | 6 +- src/gromacs/mdrun/shellfc.cpp | 5 +- src/gromacs/mdrun/tpi.cpp | 2 +- src/gromacs/mimic/communicator.cpp | 30 +++--- .../parrinellorahmanbarostat.cpp | 4 +- .../velocityscalingtemperaturecoupling.cpp | 2 +- src/gromacs/pbcutil/pbc.cpp | 4 +- src/gromacs/pulling/pull.cpp | 5 +- src/gromacs/selection/params.cpp | 3 +- src/gromacs/selection/sm_insolidangle.cpp | 6 +- src/gromacs/tables/forcetable.cpp | 6 +- src/gromacs/tools/check.cpp | 4 +- src/gromacs/tools/pme_error.cpp | 6 +- .../trajectoryanalysis/modules/angle.cpp | 2 +- .../trajectoryanalysis/modules/freevolume.cpp | 2 +- .../trajectoryanalysis/modules/sasa.cpp | 2 +- 85 files changed, 492 insertions(+), 464 deletions(-) diff --git a/api/legacy/include/gromacs/math/units.h b/api/legacy/include/gromacs/math/units.h index a83a006c84..cf869851ef 100644 --- a/api/legacy/include/gromacs/math/units.h +++ b/api/legacy/include/gromacs/math/units.h @@ -38,69 +38,70 @@ #ifndef GMX_MATH_UNITS_H #define GMX_MATH_UNITS_H +#include + /* * Physical constants to be used in Gromacs. * No constants (apart from 0, 1 or 2) should * be anywhere else in the code. */ -#include "gromacs/math/utilities.h" - -#define ANGSTROM (1e-10) /* Old... */ -#define KILO (1e3) /* Thousand */ -#define NANO (1e-9) /* A Number */ -#define PICO (1e-12) /* A Number */ -#define A2NM (ANGSTROM / NANO) /* NANO */ -#define NM2A (NANO / ANGSTROM) /* 10.0 */ -#define RAD2DEG (180.0 / M_PI) /* Conversion */ -#define DEG2RAD (M_PI / 180.0) /* id */ -#define CAL2JOULE (4.184) /* Exact definition of the calorie */ -#define E_CHARGE (1.602176634e-19) /* Exact definition, Coulomb NIST 2018 CODATA */ - -#define AMU (1.66053906660e-27) /* kg, NIST 2018 CODATA */ -#define BOLTZMANN (1.380649e-23) /* (J/K, Exact definition, NIST 2018 CODATA */ -#define AVOGADRO (6.02214076e23) /* 1/mol, Exact definition, NIST 2018 CODATA */ -#define RGAS (BOLTZMANN * AVOGADRO) /* (J/(mol K)) */ -#define BOLTZ (RGAS / KILO) /* (kJ/(mol K)) */ -#define FARADAY (E_CHARGE * AVOGADRO) /* (C/mol) */ -#define ELECTRONVOLT (E_CHARGE * AVOGADRO / KILO) /* (kJ/mol) */ -#define PLANCK1 (6.62607015e-34) /* J/Hz, Exact definition, NIST 2018 CODATA */ -#define PLANCK (PLANCK1 * AVOGADRO / (PICO * KILO)) /* (kJ/mol) ps */ - -#define EPSILON0_SI (8.8541878128e-12) /* F/m, NIST 2018 CODATA */ +namespace gmx +{ + +constexpr double c_angstrom = 1e-10; +constexpr double c_kilo = 1e3; +constexpr double c_nano = 1e-9; +constexpr double c_pico = 1e-12; +constexpr double c_nm2A = c_nano / c_angstrom; +constexpr double c_rad2Deg = 180.0 / M_PI; +constexpr double c_deg2Rad = M_PI / 180.0; +constexpr double c_cal2Joule = 4.184; /* Exact definition of the calorie */ +constexpr double c_electronCharge = 1.602176634e-19; /* Exact definition, Coulomb NIST 2018 CODATA */ + +constexpr double c_amu = 1.66053906660e-27; /* kg, NIST 2018 CODATA */ +constexpr double c_boltzmann = 1.380649e-23; /* (J/K, Exact definition, NIST 2018 CODATA */ +constexpr double c_avogadro = 6.02214076e23; /* 1/mol, Exact definition, NIST 2018 CODATA */ +constexpr double c_universalGasConstant = c_boltzmann * c_avogadro; /* (J/(mol K)) */ +constexpr double c_boltz = c_universalGasConstant / c_kilo; /* (kJ/(mol K)) */ +constexpr double c_faraday = c_electronCharge * c_avogadro; /* (C/mol) */ +constexpr double c_planck1 = 6.62607015e-34; /* J/Hz, Exact definition, NIST 2018 CODATA */ +constexpr double c_planck = (c_planck1 * c_avogadro / (c_pico * c_kilo)); /* (kJ/mol) ps */ + +constexpr double c_epsilon0Si = 8.8541878128e-12; /* F/m, NIST 2018 CODATA */ /* Epsilon in our MD units: (e^2 / Na (kJ nm)) == (e^2 mol/(kJ nm)) */ -#define EPSILON0 ((EPSILON0_SI * NANO * KILO) / (E_CHARGE * E_CHARGE * AVOGADRO)) +constexpr double c_epsilon0 = + ((c_epsilon0Si * c_nano * c_kilo) / (c_electronCharge * c_electronCharge * c_avogadro)); -#define SPEED_OF_LIGHT (2.99792458e05) /* units of nm/ps, Exact definition, NIST 2018 CODATA */ -#define ATOMICMASS_keV (931494.10242) /* Atomic mass in keV, NIST 2018 CODATA */ -#define ELECTRONMASS_keV (510.99895000) /* Electron mas in keV, NIST 2018 CODATA */ +constexpr double c_speedOfLight = + 2.99792458e05; /* units of nm/ps, Exact definition, NIST 2018 CODATA */ -#define RYDBERG (1.0973731568160e-02) /* nm^-1, NIST 2018 CODATA */ +constexpr double c_rydberg = 1.0973731568160e-02; /* nm^-1, NIST 2018 CODATA */ -#define ONE_4PI_EPS0 (1.0 / (4.0 * M_PI * EPSILON0)) -#define FACEL (10.0 * ONE_4PI_EPS0) +constexpr double c_one4PiEps0 = (1.0 / (4.0 * M_PI * c_epsilon0)); /* Pressure in MD units is: - * 1 bar = 1e5 Pa = 1e5 kg m^-1 s^-2 = 1e-28 kg nm^-1 ps^-2 = 1e-28 / AMU amu nm^1 ps ^2 + * 1 bar = 1e5 Pa = 1e5 kg m^-1 s^-2 = 1e-28 kg nm^-1 ps^-2 = 1e-28 / amu amu nm^1 ps ^2 */ -#define BAR_MDUNITS (1e5 * NANO * PICO * PICO / AMU) -#define PRESFAC (1.0 / BAR_MDUNITS) +constexpr double c_barMdunits = (1e5 * c_nano * c_pico * c_pico / c_amu); +constexpr double c_presfac = 1.0 / c_barMdunits; -/* DEBYE2ENM should be (1e-21*PICO)/(SPEED_OF_LIGHT*E_CHARGE*NANO*NANO), +/* c_debye2Enm should be (1e-21*c_pico)/(c_speedOfLight*c_electronCharge*c_nano*c_nano), * but we need to factor out some of the exponents to avoid single-precision overflows. */ -#define DEBYE2ENM (1e-15 / (SPEED_OF_LIGHT * E_CHARGE)) -#define ENM2DEBYE (1.0 / DEBYE2ENM) +constexpr double c_debye2Enm = (1e-15 / (c_speedOfLight * c_electronCharge)); +constexpr double c_enm2Debye = 1.0 / c_debye2Enm; /* to convert from a acceleration in (e V)/(amu nm) */ -/* FIELDFAC is also Faraday's constant and E_CHARGE/(1e6 AMU) */ -#define FIELDFAC (FARADAY / KILO) +/* c_fieldfac is also Faraday's constant and c_electronCharge/(1e6 amu) */ +constexpr double c_fieldfac = c_faraday / c_kilo; /* to convert AU to MD units: */ -#define HARTREE2KJ ((2.0 * RYDBERG * PLANCK * SPEED_OF_LIGHT) / AVOGADRO) -#define BOHR2NM (0.0529177210903) /* nm^-1, NIST 2018 CODATA */ -#define HARTREE_BOHR2MD (HARTREE2KJ * AVOGADRO / BOHR2NM) +constexpr double c_hartree2Kj = ((2.0 * c_rydberg * c_planck * c_speedOfLight) / c_avogadro); +constexpr double c_bohr2Nm = 0.0529177210903; /* nm^-1, NIST 2018 CODATA */ +constexpr double c_hartreeBohr2Md = (c_hartree2Kj * c_avogadro / c_bohr2Nm); +} // namespace gmx /* The four basic units */ #define unit_length "nm" @@ -108,16 +109,16 @@ #define unit_mass "u" #define unit_energy "kJ/mol" -/* Temperature unit, T in this unit times BOLTZ give energy in unit_energy */ +/* Temperature unit, T in this unit times c_boltz give energy in unit_energy */ #define unit_temp_K "K" -/* Charge unit, electron charge, involves ONE_4PI_EPS0 */ +/* Charge unit, electron charge, involves c_one4PiEps0 */ #define unit_charge_e "e" -/* Pressure unit, pressure in basic units times PRESFAC gives this unit */ +/* Pressure unit, pressure in basic units times c_presfac gives this unit */ #define unit_pres_bar "bar" -/* Dipole unit, debye, conversion from the unit_charge_e involves ENM2DEBYE */ +/* Dipole unit, debye, conversion from the unit_charge_e involves c_enm2Debye */ #define unit_dipole_D "D" /* Derived units from basic units only */ @@ -128,7 +129,7 @@ /* Other derived units */ #define unit_surft_bar unit_pres_bar " " unit_length -/* SI units, conversion from basic units involves NANO, PICO and AMU */ +/* SI units, conversion from basic units involves c_nano, c_pico and amu */ #define unit_length_SI "m" #define unit_time_SI "s" #define unit_mass_SI "kg" diff --git a/src/gromacs/applied_forces/awh/awh.cpp b/src/gromacs/applied_forces/awh/awh.cpp index 8c80087391..d55853d933 100644 --- a/src/gromacs/applied_forces/awh/awh.cpp +++ b/src/gromacs/applied_forces/awh/awh.cpp @@ -215,7 +215,7 @@ Awh::Awh(FILE* fplog, } /* Initialize all the biases */ - const double beta = 1 / (BOLTZ * inputRecord.opts.ref_t[0]); + const double beta = 1 / (gmx::c_boltz * inputRecord.opts.ref_t[0]); for (int k = 0; k < awhParams.numBias; k++) { const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams[k]; @@ -240,7 +240,7 @@ Awh::Awh(FILE* fplog, GMX_THROW(InvalidInputError( "Pull geometry 'direction-periodic' is not supported by AWH")); } - double conversionFactor = pull_coordinate_is_angletype(&pullCoord) ? DEG2RAD : 1; + double conversionFactor = pull_coordinate_is_angletype(&pullCoord) ? gmx::c_deg2Rad : 1; pullCoordIndex.push_back(awhDimParams.coordIndex); dimParams.push_back(DimParams::pullDimParams( conversionFactor, awhDimParams.forceConstant, beta)); diff --git a/src/gromacs/applied_forces/awh/dimparams.h b/src/gromacs/applied_forces/awh/dimparams.h index 2b3ac01720..11fc93f83c 100644 --- a/src/gromacs/applied_forces/awh/dimparams.h +++ b/src/gromacs/applied_forces/awh/dimparams.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -103,7 +103,7 @@ public: * Builder function for pull dimension parameters. * * \param[in] conversionFactor Conversion factor from user coordinate units to bias internal - * units (=DEG2RAD for angles). + * units (=c_deg2Rad for angles). * \param[in] forceConstant The harmonic force constant. * \param[in] beta 1/(k_B T). */ diff --git a/src/gromacs/applied_forces/electricfield.cpp b/src/gromacs/applied_forces/electricfield.cpp index 5be37f1de1..6c4381747a 100644 --- a/src/gromacs/applied_forces/electricfield.cpp +++ b/src/gromacs/applied_forces/electricfield.cpp @@ -312,7 +312,7 @@ void ElectricField::calculateForces(const ForceProviderInput& forceProviderInput for (int m = 0; (m < DIM); m++) { - const real fieldStrength = FIELDFAC * field(m, t); + const real fieldStrength = gmx::c_fieldfac * field(m, t); if (fieldStrength != 0) { diff --git a/src/gromacs/ewald/ewald.cpp b/src/gromacs/ewald/ewald.cpp index 5f67078b68..26185de74a 100644 --- a/src/gromacs/ewald/ewald.cpp +++ b/src/gromacs/ewald/ewald.cpp @@ -173,7 +173,8 @@ real do_ewald(const t_inputrec& ir, } /* 1/(Vol*e0) */ - real scaleRecip = 4.0 * M_PI / (boxDiag[XX] * boxDiag[YY] * boxDiag[ZZ]) * ONE_4PI_EPS0 / ir.epsilon_r; + real scaleRecip = + 4.0 * M_PI / (boxDiag[XX] * boxDiag[YY] * boxDiag[ZZ]) * gmx::c_one4PiEps0 / ir.epsilon_r; snew(eir, et->kmax); for (n = 0; n < et->kmax; n++) @@ -328,7 +329,7 @@ real ewald_charge_correction(const t_commrec* cr, /* Apply charge correction */ vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ]; - fac = M_PI * ONE_4PI_EPS0 + fac = M_PI * gmx::c_one4PiEps0 / (fr->ic->epsilon_r * 2.0 * vol * vol * gmx::square(fr->ic->ewaldcoeff_q)); qs2A = fr->qsum[0] * fr->qsum[0]; diff --git a/src/gromacs/ewald/long_range_correction.cpp b/src/gromacs/ewald/long_range_correction.cpp index a746ae8bc8..ce2625717c 100644 --- a/src/gromacs/ewald/long_range_correction.cpp +++ b/src/gromacs/ewald/long_range_correction.cpp @@ -98,7 +98,7 @@ void ewald_LRcorrection(const int numAtomsLocal, EwaldBoxZScaler boxScaler(ir); boxScaler.scaleBox(box, scaledBox); - one_4pi_eps = ONE_4PI_EPS0 / fr.ic->epsilon_r; + one_4pi_eps = gmx::c_one4PiEps0 / fr.ic->epsilon_r; Vexcl_q = 0; dvdl_excl_q = 0; Vdipole[0] = 0; @@ -109,8 +109,8 @@ void ewald_LRcorrection(const int numAtomsLocal, */ for (i = 0; (i < DIM); i++) { - mutot[0][i] = mu_tot[0][i] * DEBYE2ENM; - mutot[1][i] = mu_tot[1][i] * DEBYE2ENM; + mutot[0][i] = mu_tot[0][i] * gmx::c_debye2Enm; + mutot[1][i] = mu_tot[1][i] * gmx::c_debye2Enm; dipcorrA[i] = 0; dipcorrB[i] = 0; } @@ -122,7 +122,7 @@ void ewald_LRcorrection(const int numAtomsLocal, case EwaldGeometry::ThreeD: if (ir.epsilon_surface != 0) { - dipole_coeff = 2 * M_PI * ONE_4PI_EPS0 + dipole_coeff = 2 * M_PI * gmx::c_one4PiEps0 / ((2 * ir.epsilon_surface + fr.ic->epsilon_r) * boxVolume); for (i = 0; (i < DIM); i++) { diff --git a/src/gromacs/ewald/pme_gpu_internal.cpp b/src/gromacs/ewald/pme_gpu_internal.cpp index 0d7ee39e23..028a66a35a 100644 --- a/src/gromacs/ewald/pme_gpu_internal.cpp +++ b/src/gromacs/ewald/pme_gpu_internal.cpp @@ -898,7 +898,7 @@ static void pme_gpu_init(gmx_pme_t* pme, GMX_ASSERT(pmeGpu->common->epsilon_r != 0.0F, "PME GPU: bad electrostatic coefficient"); auto* kernelParamsPtr = pme_gpu_get_kernel_params_base_ptr(pmeGpu); - kernelParamsPtr->constants.elFactor = ONE_4PI_EPS0 / pmeGpu->common->epsilon_r; + kernelParamsPtr->constants.elFactor = gmx::c_one4PiEps0 / pmeGpu->common->epsilon_r; } void pme_gpu_get_real_grid_sizes(const PmeGpu* pmeGpu, gmx::IVec* gridSize, gmx::IVec* paddedGridSize) diff --git a/src/gromacs/ewald/pme_gpu_types.h b/src/gromacs/ewald/pme_gpu_types.h index 4274f095d9..abf7a17ed0 100644 --- a/src/gromacs/ewald/pme_gpu_types.h +++ b/src/gromacs/ewald/pme_gpu_types.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2016,2017,2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -100,7 +100,7 @@ static_assert(sizeof(DeviceBuffer) == 8, */ struct PmeGpuConstParams { - /*! \brief Electrostatics coefficient = ONE_4PI_EPS0 / pme->epsilon_r */ + /*! \brief Electrostatics coefficient = c_one4PiEps0 / pme->epsilon_r */ float elFactor; /*! \brief Virial and energy GPU array. Size is c_virialAndEnergyCount (7) floats. * The element order is virxx, viryy, virzz, virxy, virxz, viryz, energy. */ diff --git a/src/gromacs/ewald/pme_solve.cpp b/src/gromacs/ewald/pme_solve.cpp index 03ca35bfee..416f67b3e1 100644 --- a/src/gromacs/ewald/pme_solve.cpp +++ b/src/gromacs/ewald/pme_solve.cpp @@ -352,7 +352,7 @@ int solve_pme_yzx(const gmx_pme_t* pme, t_complex* grid, real vol, bool computeE ivec local_ndata, local_offset, local_size; real elfac; - elfac = ONE_4PI_EPS0 / pme->epsilon_r; + elfac = gmx::c_one4PiEps0 / pme->epsilon_r; nx = pme->nkx; ny = pme->nky; diff --git a/src/gromacs/fileio/pdbio.cpp b/src/gromacs/fileio/pdbio.cpp index 530459e2ab..2025fd17ae 100644 --- a/src/gromacs/fileio/pdbio.cpp +++ b/src/gromacs/fileio/pdbio.cpp @@ -104,7 +104,7 @@ void gmx_write_pdb_box(FILE* out, PbcType pbcType, const matrix box) if (norm2(box[YY]) * norm2(box[ZZ]) != 0) { - alpha = RAD2DEG * gmx_angle(box[YY], box[ZZ]); + alpha = gmx::c_rad2Deg * gmx_angle(box[YY], box[ZZ]); } else { @@ -112,7 +112,7 @@ void gmx_write_pdb_box(FILE* out, PbcType pbcType, const matrix box) } if (norm2(box[XX]) * norm2(box[ZZ]) != 0) { - beta = RAD2DEG * gmx_angle(box[XX], box[ZZ]); + beta = gmx::c_rad2Deg * gmx_angle(box[XX], box[ZZ]); } else { @@ -120,7 +120,7 @@ void gmx_write_pdb_box(FILE* out, PbcType pbcType, const matrix box) } if (norm2(box[XX]) * norm2(box[YY]) != 0) { - gamma = RAD2DEG * gmx_angle(box[XX], box[YY]); + gamma = gmx::c_rad2Deg * gmx_angle(box[XX], box[YY]); } else { @@ -206,7 +206,7 @@ static void read_cryst1(char* line, PbcType* pbcType, matrix box) { if (alpha != 90.0) { - cosa = std::cos(alpha * DEG2RAD); + cosa = std::cos(alpha * gmx::c_deg2Rad); } else { @@ -214,7 +214,7 @@ static void read_cryst1(char* line, PbcType* pbcType, matrix box) } if (beta != 90.0) { - cosb = std::cos(beta * DEG2RAD); + cosb = std::cos(beta * gmx::c_deg2Rad); } else { @@ -222,8 +222,8 @@ static void read_cryst1(char* line, PbcType* pbcType, matrix box) } if (gamma != 90.0) { - cosg = std::cos(gamma * DEG2RAD); - sing = std::sin(gamma * DEG2RAD); + cosg = std::cos(gamma * gmx::c_deg2Rad); + sing = std::sin(gamma * gmx::c_deg2Rad); } else { diff --git a/src/gromacs/fileio/tngio.cpp b/src/gromacs/fileio/tngio.cpp index 90615238cc..485d3e35b2 100644 --- a/src/gromacs/fileio/tngio.cpp +++ b/src/gromacs/fileio/tngio.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2013,2014,2015,2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -549,7 +549,7 @@ void gmx_tng_prepare_md_writing(gmx_tng_trajectory_t gmx_tng, const gmx_mtop_t* #if GMX_USE_TNG gmx_tng_add_mtop(gmx_tng, mtop); set_writing_intervals(gmx_tng, FALSE, ir); - tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * PICO); + tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * gmx::c_pico); gmx_tng->timePerFrameIsSet = true; #else GMX_UNUSED_VALUE(gmx_tng); @@ -765,7 +765,7 @@ void gmx_tng_prepare_low_prec_writing(gmx_tng_trajectory_t gmx_tng, const gmx_mt gmx_tng_add_mtop(gmx_tng, mtop); add_selection_groups(gmx_tng, mtop); set_writing_intervals(gmx_tng, TRUE, ir); - tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * PICO); + tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * gmx::c_pico); gmx_tng->timePerFrameIsSet = true; gmx_tng_set_compression_precision(gmx_tng, ir->x_compression_precision); #else @@ -801,7 +801,7 @@ void gmx_fwrite_tng(gmx_tng_trajectory_t gmx_tng, # else static write_data_func_pointer write_data = tng_util_generic_with_time_write; # endif - double elapsedSeconds = elapsedPicoSeconds * PICO; + double elapsedSeconds = elapsedPicoSeconds * gmx::c_pico; int64_t nParticles; char compression; @@ -1004,7 +1004,7 @@ float gmx_tng_get_time_of_final_frame(gmx_tng_trajectory_t gmx_tng) tng_num_frames_get(tng, &nFrames); tng_util_time_of_frame_get(tng, nFrames - 1, &time); - fTime = time / PICO; + fTime = time / gmx::c_pico; return fTime; #else GMX_UNUSED_VALUE(gmx_tng); @@ -1271,8 +1271,8 @@ real getDistanceScaleFactor(gmx_tng_trajectory_t in) // GROMACS expects distances in nm switch (exp) { - case 9: distanceScaleFactor = NANO / NANO; break; - case 10: distanceScaleFactor = NANO / ANGSTROM; break; + case 9: distanceScaleFactor = gmx::c_nano / gmx::c_nano; break; + case 10: distanceScaleFactor = gmx::c_nano / gmx::c_angstrom; break; default: distanceScaleFactor = pow(10.0, exp + 9.0); } @@ -1545,7 +1545,7 @@ gmx_bool gmx_read_next_tng_frame(gmx_tng_trajectory_t gmx_tng_input, fr->bStep = TRUE; // Convert the time to ps - fr->time = frameTime / PICO; + fr->time = frameTime / gmx::c_pico; fr->bTime = (frameTime > 0); // TODO This does not leak, but is not exception safe. diff --git a/src/gromacs/gmxana/anadih.cpp b/src/gromacs/gmxana/anadih.cpp index e08ccbc075..4e1f89690b 100644 --- a/src/gromacs/gmxana/anadih.cpp +++ b/src/gromacs/gmxana/anadih.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -107,7 +107,7 @@ static int calc_RBbin(real phi, int gmx_unused multiplicity, real gmx_unused cor static int calc_Nbin(real phi, int multiplicity, real core_frac) { - static const real r360 = 360 * DEG2RAD; + static const real r360 = 360 * gmx::c_deg2Rad; real rot_width, core_width, core_offset, low, hi; int bin; /* with multiplicity 3 and core_frac 0.5 @@ -127,8 +127,8 @@ static int calc_Nbin(real phi, int multiplicity, real core_frac) { low = ((bin - 1) * rot_width) + core_offset; hi = ((bin - 1) * rot_width) + core_offset + core_width; - low *= DEG2RAD; - hi *= DEG2RAD; + low *= gmx::c_deg2Rad; + hi *= gmx::c_deg2Rad; if ((phi > low) && (phi < hi)) { return bin; @@ -716,7 +716,7 @@ static real calc_fraction(const real angles[], int nangles) for (i = 0; i < nangles; i++) { - angle = angles[i] * RAD2DEG; + angle = angles[i] * gmx::c_rad2Deg; if (angle > 135 && angle < 225) { diff --git a/src/gromacs/gmxana/angle_correction.h b/src/gromacs/gmxana/angle_correction.h index a1019ad1da..feed719fcf 100644 --- a/src/gromacs/gmxana/angle_correction.h +++ b/src/gromacs/gmxana/angle_correction.h @@ -35,7 +35,6 @@ #ifndef GMX_GMXANA_ANGLE_CORRECTION_H #define GMX_GMXANA_ANGLE_CORRECTION_H -#include "gromacs/math/units.h" #include "gromacs/utility/real.h" /*! \brief diff --git a/src/gromacs/gmxana/gmx_angle.cpp b/src/gromacs/gmxana/gmx_angle.cpp index 4fe720ba9d..684a7b15c3 100644 --- a/src/gromacs/gmxana/gmx_angle.cpp +++ b/src/gromacs/gmxana/gmx_angle.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -312,7 +312,7 @@ int gmx_g_angle(int argc, char* argv[]) out = xvgropen(opt2fn("-ov", NFILE, fnm), title, "Time (ps)", "Angle (degrees)", oenv); for (i = 0; (i < nframes); i++) { - fprintf(out, "%10.5f %8.3f", time[i], aver_angle[i] * RAD2DEG); + fprintf(out, "%10.5f %8.3f", time[i], aver_angle[i] * gmx::c_rad2Deg); if (bALL) { for (j = 0; (j < nangles); j++) @@ -320,11 +320,11 @@ int gmx_g_angle(int argc, char* argv[]) if (bPBC) { real dd = dih[j][i]; - fprintf(out, " %8.3f", std::atan2(std::sin(dd), std::cos(dd)) * RAD2DEG); + fprintf(out, " %8.3f", std::atan2(std::sin(dd), std::cos(dd)) * gmx::c_rad2Deg); } else { - fprintf(out, " %8.3f", dih[j][i] * RAD2DEG); + fprintf(out, " %8.3f", dih[j][i] * gmx::c_rad2Deg); } } } @@ -372,7 +372,7 @@ int gmx_g_angle(int argc, char* argv[]) if (bChandler) { - real dval, sixty = DEG2RAD * 60; + real dval, sixty = gmx::c_deg2Rad * 60; gmx_bool bTest; for (i = 0; (i < nangles); i++) @@ -455,8 +455,8 @@ int gmx_g_angle(int argc, char* argv[]) } aver /= nframes; double aversig = correctRadianAngleRange(aver); - aversig *= RAD2DEG; - aver *= RAD2DEG; + aversig *= gmx::c_rad2Deg; + aver *= gmx::c_rad2Deg; printf(" < angle > = %g\n", aversig); if (mult == 3) diff --git a/src/gromacs/gmxana/gmx_awh.cpp b/src/gromacs/gmxana/gmx_awh.cpp index 00b0b9a542..690d932eea 100644 --- a/src/gromacs/gmxana/gmx_awh.cpp +++ b/src/gromacs/gmxana/gmx_awh.cpp @@ -2,7 +2,7 @@ * This file is part of the GROMACS molecular simulation package. * * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -630,8 +630,13 @@ int gmx_awh(int argc, char* argv[]) AwhGraphSelection awhGraphSelection = (moreGraphs ? AwhGraphSelection::All : AwhGraphSelection::Pmf); EnergyUnit energyUnit = (kTUnit ? EnergyUnit::KT : EnergyUnit::KJPerMol); - awhReader = std::make_unique( - ir.awhParams, nfile, fnm, awhGraphSelection, energyUnit, BOLTZ * ir.opts.ref_t[0], block); + awhReader = std::make_unique(ir.awhParams, + nfile, + fnm, + awhGraphSelection, + energyUnit, + gmx::c_boltz * ir.opts.ref_t[0], + block); } awhReader->processAwhFrame(*block, frame->t, oenv); diff --git a/src/gromacs/gmxana/gmx_bar.cpp b/src/gromacs/gmxana/gmx_bar.cpp index b70ffedf36..fa6c81011d 100644 --- a/src/gromacs/gmxana/gmx_bar.cpp +++ b/src/gromacs/gmxana/gmx_bar.cpp @@ -1640,7 +1640,7 @@ static double calc_bar_lowlevel(sample_coll_t* ca, sample_coll_t* cb, double tem double DG0, DG1, DG2, dDG1; double n1, n2; /* numbers of samples as doubles */ - kT = BOLTZ * temp; + kT = gmx::c_boltz * temp; beta = 1 / kT; /* count the numbers of samples */ @@ -1771,7 +1771,7 @@ static void calc_rel_entropy(sample_coll_t* ca, sample_coll_t* cb, double temp, double Wfac1, Wfac2; double n1, n2; - kT = BOLTZ * temp; + kT = gmx::c_boltz * temp; beta = 1 / kT; /* count the numbers of samples */ @@ -1885,7 +1885,7 @@ static void calc_dg_stddev(sample_coll_t* ca, sample_coll_t* cb, double temp, do double Wfac1, Wfac2; double n1, n2; - kT = BOLTZ * temp; + kT = gmx::c_boltz * temp; beta = 1 / kT; /* count the numbers of samples */ @@ -3622,7 +3622,7 @@ int gmx_bar(int argc, char* argv[]) } /* print results in kT */ - kT = BOLTZ * temp; + kT = gmx::c_boltz * temp; printf("\nTemperature: %g K\n", temp); diff --git a/src/gromacs/gmxana/gmx_bundle.cpp b/src/gromacs/gmxana/gmx_bundle.cpp index 07564ee2d5..c4fa316f92 100644 --- a/src/gromacs/gmxana/gmx_bundle.cpp +++ b/src/gromacs/gmxana/gmx_bundle.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -391,27 +391,27 @@ int gmx_bundle(int argc, char* argv[]) fprintf(flen, " %6g", bun.len[i]); fprintf(fdist, " %6g", norm(bun.mid[i])); fprintf(fz, " %6g", bun.mid[i][ZZ]); - fprintf(ftilt, " %6g", RAD2DEG * acos(bun.dir[i][ZZ])); + fprintf(ftilt, " %6g", gmx::c_rad2Deg * acos(bun.dir[i][ZZ])); comp = bun.mid[i][XX] * bun.dir[i][XX] + bun.mid[i][YY] * bun.dir[i][YY]; - fprintf(ftiltr, " %6g", RAD2DEG * std::asin(comp / std::hypot(comp, bun.dir[i][ZZ]))); + fprintf(ftiltr, " %6g", gmx::c_rad2Deg * std::asin(comp / std::hypot(comp, bun.dir[i][ZZ]))); comp = bun.mid[i][YY] * bun.dir[i][XX] - bun.mid[i][XX] * bun.dir[i][YY]; - fprintf(ftiltl, " %6g", RAD2DEG * std::asin(comp / std::hypot(comp, bun.dir[i][ZZ]))); + fprintf(ftiltl, " %6g", gmx::c_rad2Deg * std::asin(comp / std::hypot(comp, bun.dir[i][ZZ]))); if (bKink) { rvec_sub(bun.end[0][i], bun.end[2][i], va); rvec_sub(bun.end[2][i], bun.end[1][i], vb); unitv(va, va); unitv(vb, vb); - fprintf(fkink, " %6g", RAD2DEG * acos(iprod(va, vb))); + fprintf(fkink, " %6g", gmx::c_rad2Deg * acos(iprod(va, vb))); cprod(va, vb, vc); copy_rvec(bun.mid[i], vr); vr[ZZ] = 0; unitv(vr, vr); - fprintf(fkinkr, " %6g", RAD2DEG * std::asin(iprod(vc, vr))); + fprintf(fkinkr, " %6g", gmx::c_rad2Deg * std::asin(iprod(vc, vr))); vl[XX] = vr[YY]; vl[YY] = -vr[XX]; vl[ZZ] = 0; - fprintf(fkinkl, " %6g", RAD2DEG * std::asin(iprod(vc, vl))); + fprintf(fkinkl, " %6g", gmx::c_rad2Deg * std::asin(iprod(vc, vl))); } } fprintf(flen, "\n"); diff --git a/src/gromacs/gmxana/gmx_chi.cpp b/src/gromacs/gmxana/gmx_chi.cpp index 9cafe37eac..a23ff444a6 100644 --- a/src/gromacs/gmxana/gmx_chi.cpp +++ b/src/gromacs/gmxana/gmx_chi.cpp @@ -130,7 +130,7 @@ static gmx_bool bAllowed(real phi, real psi) #define NPP asize(map) int x, y; -#define INDEX(ppp) (((static_cast(360 + (ppp)*RAD2DEG)) % 360) / 6) +#define INDEX(ppp) (((static_cast(360 + (ppp)*gmx::c_rad2Deg)) % 360) / 6) x = INDEX(phi); y = INDEX(psi); #undef INDEX @@ -932,16 +932,18 @@ static void do_rama(int nf, Psi = dlist[i].j0[edPsi]; for (j = 0; (j < nf); j++) { - phi = RAD2DEG * dih[Phi][j]; - psi = RAD2DEG * dih[Psi][j]; + phi = gmx::c_rad2Deg * dih[Phi][j]; + psi = gmx::c_rad2Deg * dih[Psi][j]; fprintf(fp, "%10g %10g\n", phi, psi); if (bViol) { - fprintf(gp, "%d\n", static_cast(!bAllowed(dih[Phi][j], RAD2DEG * dih[Psi][j]))); + fprintf(gp, + "%d\n", + static_cast(!bAllowed(dih[Phi][j], gmx::c_rad2Deg * dih[Psi][j]))); } if (bOm) { - omega = RAD2DEG * dih[Om][j]; + omega = gmx::c_rad2Deg * dih[Om][j]; mat[static_cast(((phi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))] [static_cast(((psi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))] += omega; } @@ -1023,7 +1025,7 @@ static void do_rama(int nf, Xi2 = dlist[i].j0[edChi2]; for (j = 0; (j < nf); j++) { - fprintf(fp, "%10g %10g\n", RAD2DEG * dih[Xi1][j], RAD2DEG * dih[Xi2][j]); + fprintf(fp, "%10g %10g\n", gmx::c_rad2Deg * dih[Xi1][j], gmx::c_rad2Deg * dih[Xi2][j]); } xvgrclose(fp); } diff --git a/src/gromacs/gmxana/gmx_clustsize.cpp b/src/gromacs/gmxana/gmx_clustsize.cpp index b98b5b1db9..44ec18ba8c 100644 --- a/src/gromacs/gmxana/gmx_clustsize.cpp +++ b/src/gromacs/gmxana/gmx_clustsize.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2007, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -354,7 +354,7 @@ static void clust_size(const char* ndx, ekin += 0.5 * m * iprod(v[ai], v[ai]); } } - temp = (ekin * 2.0) / (3.0 * tfac * max_clust_size * BOLTZ); + temp = (ekin * 2.0) / (3.0 * tfac * max_clust_size * gmx::c_boltz); fprintf(tp, "%10.3f %10.3f\n", frameTime, temp); } } diff --git a/src/gromacs/gmxana/gmx_current.cpp b/src/gromacs/gmxana/gmx_current.cpp index 228ddaba0c..78c58596a9 100644 --- a/src/gromacs/gmxana/gmx_current.cpp +++ b/src/gromacs/gmxana/gmx_current.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 2008,2009,2010,2011,2012 by the GROMACS development team. * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -58,7 +58,8 @@ #include "gromacs/utility/gmxassert.h" #include "gromacs/utility/smalloc.h" -#define EPSI0 (EPSILON0 * E_CHARGE * E_CHARGE * AVOGADRO / (KILO * NANO)) /* EPSILON0 in SI units */ +constexpr double EPSI0 = (gmx::c_epsilon0 * gmx::c_electronCharge * gmx::c_electronCharge * gmx::c_avogadro + / (gmx::c_kilo * gmx::c_nano)); /* c_epsilon0 in SI units */ static void index_atom2mol(int* n, int* index, t_block* mols) { @@ -640,11 +641,11 @@ static void dielectric(FILE* fmj, volume_av /= refr; prefactor = 1.0; - prefactor /= 3.0 * EPSILON0 * volume_av * BOLTZ * temp; + prefactor /= 3.0 * gmx::c_epsilon0 * volume_av * gmx::c_boltz * temp; - prefactorav = E_CHARGE * E_CHARGE; - prefactorav /= volume_av * BOLTZMANN * temp * NANO * 6.0; + prefactorav = gmx::c_electronCharge * gmx::c_electronCharge; + prefactorav /= volume_av * gmx::c_boltzmann * temp * gmx::c_nano * 6.0; fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav); @@ -694,7 +695,7 @@ static void dielectric(FILE* fmj, { printf("\nCalculating current autocorrelation ... \n"); - sgk = calc_cacf(outf, prefactorav / PICO, cacf, time, nvfr, vfr, ie, nshift); + sgk = calc_cacf(outf, prefactorav / gmx::c_pico, cacf, time, nvfr, vfr, ie, nshift); if (ie > ii) { diff --git a/src/gromacs/gmxana/gmx_density.cpp b/src/gromacs/gmxana/gmx_density.cpp index 950e69b231..cade69fe5d 100644 --- a/src/gromacs/gmxana/gmx_density.cpp +++ b/src/gromacs/gmxana/gmx_density.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -589,7 +589,7 @@ static void plot_density(double* slDensity[], } if (dens_opt[0][0] == 'm') { - fprintf(den, " %12g", ddd * AMU / (NANO * NANO * NANO)); + fprintf(den, " %12g", ddd * gmx::c_amu / (gmx::c_nano * gmx::c_nano * gmx::c_nano)); } else { diff --git a/src/gromacs/gmxana/gmx_densorder.cpp b/src/gromacs/gmxana/gmx_densorder.cpp index f512024c4d..b9bebf3c8b 100644 --- a/src/gromacs/gmxana/gmx_densorder.cpp +++ b/src/gromacs/gmxana/gmx_densorder.cpp @@ -2,7 +2,7 @@ * This file is part of the GROMACS molecular simulation package. * * Copyright (c) 2010-2018, The GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -228,8 +228,9 @@ static void density_in_time(const char* fn, (*Densdevel)[*tblock] = Densslice; } - dscale = (*xslices) * (*yslices) * (*zslices) * AMU - / (box[ax1][ax1] * box[ax2][ax2] * box[axis][axis] * nsttblock * (NANO * NANO * NANO)); + dscale = (*xslices) * (*yslices) * (*zslices) * gmx::c_amu + / (box[ax1][ax1] * box[ax2][ax2] * box[axis][axis] * nsttblock + * (gmx::c_nano * gmx::c_nano * gmx::c_nano)); if (bCenter) { diff --git a/src/gromacs/gmxana/gmx_dipoles.cpp b/src/gromacs/gmxana/gmx_dipoles.cpp index 7e335775a4..5f1d64283b 100644 --- a/src/gromacs/gmxana/gmx_dipoles.cpp +++ b/src/gromacs/gmxana/gmx_dipoles.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -74,9 +74,9 @@ #include "gromacs/utility/gmxassert.h" #include "gromacs/utility/smalloc.h" -#define e2d(x) ENM2DEBYE*(x) -#define EANG2CM (E_CHARGE * 1.0e-10) /* e Angstrom to Coulomb meter */ -#define CM2D (SPEED_OF_LIGHT * 1.0e+24) /* Coulomb meter to Debye */ +#define e2d(x) gmx::c_enm2Debye*(x) +constexpr double EANG2CM = gmx::c_electronCharge * 1.0e-10; /* e Angstrom to Coulomb meter */ +constexpr double CM2D = gmx::c_speedOfLight * 1.0e+24; /* Coulomb meter to Debye */ typedef struct { @@ -440,7 +440,7 @@ static void print_gkrbin(const char* fn, t_gkrbin* gb, int ngrp, int nframes, re { cosav = 0; } - ener = -0.5 * cosav * ONE_4PI_EPS0 / (x1 * x1 * x1); + ener = -0.5 * cosav * gmx::c_one4PiEps0 / (x1 * x1 * x1); fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n", x1, Gkr, cosav, hOO, gOO, ener); @@ -673,7 +673,8 @@ static real calc_eps(double M_diff, double volume, double epsRF, double temp) double eps_0 = 8.854187817e-12; /* epsilon_0 in C^2 J^-1 m^-1 */ double fac = 1.112650021e-59; /* converts Debye^2 to C^2 m^2 */ - A = M_diff * fac / (3 * eps_0 * volume * NANO * NANO * NANO * BOLTZMANN * temp); + A = M_diff * fac + / (3 * eps_0 * volume * gmx::c_nano * gmx::c_nano * gmx::c_nano * gmx::c_boltzmann * temp); if (epsRF == 0.0) { diff --git a/src/gromacs/gmxana/gmx_dos.cpp b/src/gromacs/gmxana/gmx_dos.cpp index 1c4483c069..25d2acc2f2 100644 --- a/src/gromacs/gmxana/gmx_dos.cpp +++ b/src/gromacs/gmxana/gmx_dos.cpp @@ -2,7 +2,7 @@ * This file is part of the GROMACS molecular simulation package. * * Copyright (c) 2011-2018, The GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -190,12 +190,12 @@ static double calc_Shs(double f, double y) { double fy = f * y; - return BOLTZ * (std::log(calc_compress(fy)) + fy * (3 * fy - 4) / gmx::square(1 - fy)); + return gmx::c_boltz * (std::log(calc_compress(fy)) + fy * (3 * fy - 4) / gmx::square(1 - fy)); } static real wCsolid(real nu, real beta) { - real bhn = beta * PLANCK * nu; + real bhn = beta * gmx::c_planck * nu; real ebn, koko; if (bhn == 0) @@ -212,7 +212,7 @@ static real wCsolid(real nu, real beta) static real wSsolid(real nu, real beta) { - real bhn = beta * PLANCK * nu; + real bhn = beta * gmx::c_planck * nu; if (bhn == 0) { @@ -226,7 +226,7 @@ static real wSsolid(real nu, real beta) static real wAsolid(real nu, real beta) { - real bhn = beta * PLANCK * nu; + real bhn = beta * gmx::c_planck * nu; if (bhn == 0) { @@ -240,7 +240,7 @@ static real wAsolid(real nu, real beta) static real wEsolid(real nu, real beta) { - real bhn = beta * PLANCK * nu; + real bhn = beta * gmx::c_planck * nu; if (bhn == 0) { @@ -341,7 +341,7 @@ int gmx_dos(int argc, char* argv[]) return 0; } - beta = 1 / (Temp * BOLTZ); + beta = 1 / (Temp * gmx::c_boltz); fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w"); fprintf(fplog, "Doing density of states analysis based on trajectory.\n"); @@ -527,15 +527,16 @@ int gmx_dos(int argc, char* argv[]) DoS0 = dos[DOS][0]; /* Note this eqn. is incorrect in Pascal2011a! */ - Delta = ((2 * DoS0 / (9 * Natom)) * std::sqrt(M_PI * BOLTZ * Temp * Natom / tmass) + Delta = ((2 * DoS0 / (9 * Natom)) * std::sqrt(M_PI * gmx::c_boltz * Temp * Natom / tmass) * std::pow((Natom / V), 1.0 / 3.0) * std::pow(6.0 / M_PI, 2.0 / 3.0)); f = calc_fluidicity(Delta, toler); y = calc_y(f, Delta, toler); z = calc_compress(y); - Sig = BOLTZ - * (5.0 / 2.0 + std::log(2 * M_PI * BOLTZ * Temp / (gmx::square(PLANCK)) * V / (f * Natom))); + Sig = gmx::c_boltz + * (5.0 / 2.0 + + std::log(2 * M_PI * gmx::c_boltz * Temp / (gmx::square(gmx::c_planck)) * V / (f * Natom))); Shs = Sig + calc_Shs(f, y); - rho = (tmass * AMU) / (V * NANO * NANO * NANO); + rho = (tmass * gmx::c_amu) / (V * gmx::c_nano * gmx::c_nano * gmx::c_nano); sigHS = std::cbrt(6 * y * V / (M_PI * Natom)); fprintf(fplog, "System = \"%s\"\n", *top.name); @@ -567,7 +568,7 @@ int gmx_dos(int argc, char* argv[]) "\\f{4}S(\\f{12}n\\f{4})", oenv); xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv); - recip_fac = bRecip ? (1e7 / SPEED_OF_LIGHT) : 1.0; + recip_fac = bRecip ? (1e7 / gmx::c_speedOfLight) : 1.0; for (j = 0; (j < nframes / 4); j++) { dos[DOS_DIFF][j] = DoS0 / (1 + gmx::square(DoS0 * M_PI * nu[j] / (6 * f * Natom))); @@ -583,7 +584,7 @@ int gmx_dos(int argc, char* argv[]) /* Finally analyze the results! */ wCdiff = 0.5; - wSdiff = Shs / (3 * BOLTZ); /* Is this correct? */ + wSdiff = Shs / (3 * gmx::c_boltz); /* Is this correct? */ wEdiff = 0.5; wAdiff = wEdiff - wSdiff; for (j = 0; (j < nframes / 4); j++) @@ -598,7 +599,8 @@ int gmx_dos(int argc, char* argv[]) fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n", DiffCoeff); fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n", 1000 * DoS0 / (12 * tmass * beta)); - cP = BOLTZ * evaluate_integral(nframes / 4, nu, dos[DOS_CP], nullptr, int{ nframes / 4 }, &stddev); + cP = gmx::c_boltz + * evaluate_integral(nframes / 4, nu, dos[DOS_CP], nullptr, int{ nframes / 4 }, &stddev); fprintf(fplog, "Heat capacity %g J/mol K\n", 1000 * cP / Nmol); fprintf(fplog, "\nArrivederci!\n"); gmx_fio_fclose(fplog); diff --git a/src/gromacs/gmxana/gmx_enemat.cpp b/src/gromacs/gmxana/gmx_enemat.cpp index c03904be6d..fab46dce52 100644 --- a/src/gromacs/gmxana/gmx_enemat.cpp +++ b/src/gromacs/gmxana/gmx_enemat.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -419,7 +419,7 @@ int gmx_enemat(int argc, char* argv[]) } eaver[i] /= nenergy; } - beta = 1.0 / (BOLTZ * reftemp); + beta = 1.0 / (gmx::c_boltz * reftemp); snew(efree, ngroups); snew(edif, ngroups); for (i = 0; (i < ngroups); i++) diff --git a/src/gromacs/gmxana/gmx_energy.cpp b/src/gromacs/gmxana/gmx_energy.cpp index 47dcde2f60..b13a137695 100644 --- a/src/gromacs/gmxana/gmx_energy.cpp +++ b/src/gromacs/gmxana/gmx_energy.cpp @@ -359,7 +359,8 @@ static void einstein_visco(const char* fn, } } /* Convert to SI for the viscosity */ - fac = (V * NANO * NANO * NANO * PICO * 1e10) / (2 * BOLTZMANN * T) / (nint - i); + fac = (V * gmx::c_nano * gmx::c_nano * gmx::c_nano * gmx::c_pico * 1e10) + / (2 * gmx::c_boltzmann * T) / (nint - i); fprintf(fp0, "%10g", i * dt); for (m = 0; (m <= nsets); m++) { @@ -737,7 +738,7 @@ static void calc_fluctuation_props(FILE* fp, const char* my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" }; int ii[eNR]; - NANO3 = NANO * NANO * NANO; + NANO3 = gmx::c_nano * gmx::c_nano * gmx::c_nano; if (!bDriftCorr) { fprintf(fp, @@ -770,15 +771,15 @@ if ((ii[eVol] < nset) && (ii[eTemp] < nset)) { vv = edat->s[ii[eVol]].av * NANO3; varv = gmx::square(edat->s[ii[eVol]].rmsd * NANO3); - kappa = (varv / vv) / (BOLTZMANN * tt); + kappa = (varv / vv) / (gmx::c_boltzmann * tt); } /* Enthalpy */ hh = varh = NOTSET; if ((ii[eEnth] < nset) && (ii[eTemp] < nset)) { - hh = KILO * edat->s[ii[eEnth]].av / AVOGADRO; - varh = gmx::square(KILO * edat->s[ii[eEnth]].rmsd / AVOGADRO); - cp = AVOGADRO * ((varh / nmol) / (BOLTZMANN * tt * tt)); + hh = gmx::c_kilo * edat->s[ii[eEnth]].av / gmx::c_avogadro; + varh = gmx::square(gmx::c_kilo * edat->s[ii[eEnth]].rmsd / gmx::c_avogadro); + cp = gmx::c_avogadro * ((varh / nmol) / (gmx::c_boltzmann * tt * tt)); } /* Total energy */ if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET)) @@ -787,7 +788,7 @@ if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET)) by checking whether the enthalpy was computed. */ varet = gmx::square(edat->s[ii[eEtot]].rmsd); - cv = KILO * ((varet / nmol) / (BOLTZ * tt * tt)); + cv = gmx::c_kilo * ((varet / nmol) / (gmx::c_boltz * tt * tt)); } /* Alpha, dcp */ if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset)) @@ -797,7 +798,7 @@ if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset)) for (j = 0; (j < edat->nframes); j++) { v = edat->s[ii[eVol]].ener[j] * NANO3; - h = KILO * edat->s[ii[eEnth]].ener[j] / AVOGADRO; + h = gmx::c_kilo * edat->s[ii[eEnth]].ener[j] / gmx::c_avogadro; v_sum += v; h_sum += h; vh_sum += (v * h); @@ -805,8 +806,8 @@ if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset)) vh_aver = vh_sum / edat->nframes; v_aver = v_sum / edat->nframes; h_aver = h_sum / edat->nframes; - alpha = (vh_aver - v_aver * h_aver) / (v_aver * BOLTZMANN * tt * tt); - dcp = (v_aver * AVOGADRO / nmol) * tt * gmx::square(alpha) / (kappa); + alpha = (vh_aver - v_aver * h_aver) / (v_aver * gmx::c_boltzmann * tt * tt); + dcp = (v_aver * gmx::c_avogadro / nmol) * tt * gmx::square(alpha) / (kappa); } if (tt != NOTSET) @@ -827,16 +828,18 @@ if (tt != NOTSET) { if (varv != NOTSET) { - fprintf(fp, "varv = %10g (m^6)\n", varv * AVOGADRO / nmol); + fprintf(fp, "varv = %10g (m^6)\n", varv * gmx::c_avogadro / nmol); } } if (vv != NOTSET) { - fprintf(fp, "Volume = %10g m^3/mol\n", vv * AVOGADRO / nmol); + fprintf(fp, "Volume = %10g m^3/mol\n", vv * gmx::c_avogadro / nmol); } if (varh != NOTSET) { - fprintf(fp, "Enthalpy = %10g kJ/mol\n", hh * AVOGADRO / (KILO * nmol)); + fprintf(fp, + "Enthalpy = %10g kJ/mol\n", + hh * gmx::c_avogadro / (gmx::c_kilo * nmol)); } if (alpha != NOTSET) { @@ -1001,7 +1004,7 @@ static void analyse_ener(gmx_bool bCorr, expEtot = 0; if (bFee) { - beta = 1.0 / (BOLTZ * reftemp); + beta = 1.0 / (gmx::c_boltz * reftemp); snew(fee, nset); } for (i = 0; (i < nset); i++) @@ -1206,7 +1209,7 @@ static void analyse_ener(gmx_bool bCorr, 0.0, 0); - factor = (Vaver * 1e-26 / (BOLTZMANN * Temp)) * Dt; + factor = (Vaver * 1e-26 / (gmx::c_boltzmann * Temp)) * Dt; fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv); xvgr_legend(fp, asize(leg), leg, oenv); @@ -1370,7 +1373,7 @@ static void fec(const char* ene2fn, } fprintf(stdout, "\n%-24s %10s\n", "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A"); sum = 0; - beta = 1.0 / (BOLTZ * reftemp); + beta = 1.0 / (gmx::c_boltz * reftemp); for (i = 0; i < nset; i++) { if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0) @@ -1383,10 +1386,10 @@ static void fec(const char* ene2fn, sum += std::exp(-dE * beta); if (fp) { - fprintf(fp, "%10g %10g %10g\n", time[j], dE, -BOLTZ * reftemp * std::log(sum / (j + 1))); + fprintf(fp, "%10g %10g %10g\n", time[j], dE, -gmx::c_boltz * reftemp * std::log(sum / (j + 1))); } } - aver = -BOLTZ * reftemp * std::log(sum / nenergy); + aver = -gmx::c_boltz * reftemp * std::log(sum / nenergy); fprintf(stdout, "%-24s %10g\n", leg[i], aver); } if (fp) diff --git a/src/gromacs/gmxana/gmx_hbond.cpp b/src/gromacs/gmxana/gmx_hbond.cpp index b4795b68f5..51f34b5172 100644 --- a/src/gromacs/gmxana/gmx_hbond.cpp +++ b/src/gromacs/gmxana/gmx_hbond.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2008, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -1741,14 +1741,14 @@ static real calc_dg(real tau, real temp) { real kbt; - kbt = BOLTZ * temp; + kbt = gmx::c_boltz * temp; if (tau <= 0) { return -666; } else { - return kbt * std::log(kbt * tau / PLANCK); + return kbt * std::log(kbt * tau / gmx::c_planck); } } @@ -1865,7 +1865,7 @@ void analyse_corr(int n, compute_weighted_rates( n, t, ct, nt, kt, sigma_ct, sigma_nt, sigma_kt, &k, &kp, &sigma_k, &sigma_kp, fit_start); Q = 0; /* quality_of_fit(chi2, 2);*/ - ddg = BOLTZ * temp * sigma_k / k; + ddg = gmx::c_boltz * temp * sigma_k / k; printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n", chi2, Q); printf("The Rate and Delta G are followed by an error estimate\n"); printf("----------------------------------------------------------\n" @@ -1876,7 +1876,7 @@ void analyse_corr(int n, 1 / k, calc_dg(1 / k, temp), ddg); - ddg = BOLTZ * temp * sigma_kp / kp; + ddg = gmx::c_boltz * temp * sigma_kp / kp; printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n", kp, sigma_kp, @@ -2695,7 +2695,7 @@ int gmx_hbond(int argc, char* argv[]) /* process input */ bSelected = FALSE; - ccut = std::cos(acut * DEG2RAD); + ccut = std::cos(acut * gmx::c_deg2Rad); if (bContact) { @@ -3180,7 +3180,7 @@ int gmx_hbond(int argc, char* argv[]) "for an hbond: %f", dist); } - ang *= RAD2DEG; + ang *= gmx::c_rad2Deg; __ADIST[static_cast(ang / abin)]++; __RDIST[static_cast(dist / rbin)]++; if (!bTwo) diff --git a/src/gromacs/gmxana/gmx_nmeig.cpp b/src/gromacs/gmxana/gmx_nmeig.cpp index 560425e38e..7f308a16f0 100644 --- a/src/gromacs/gmxana/gmx_nmeig.cpp +++ b/src/gromacs/gmxana/gmx_nmeig.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -71,31 +71,31 @@ static double cv_corr(double nu, double T) { - double x = PLANCK * nu / (BOLTZ * T); + double x = gmx::c_planck * nu / (gmx::c_boltz * T); double ex = std::exp(x); if (nu <= 0) { - return BOLTZ * KILO; + return gmx::c_boltz * gmx::c_kilo; } else { - return BOLTZ * KILO * (ex * gmx::square(x) / gmx::square(ex - 1) - 1); + return gmx::c_boltz * gmx::c_kilo * (ex * gmx::square(x) / gmx::square(ex - 1) - 1); } } static double u_corr(double nu, double T) { - double x = PLANCK * nu / (BOLTZ * T); + double x = gmx::c_planck * nu / (gmx::c_boltz * T); double ex = std::exp(x); if (nu <= 0) { - return BOLTZ * T; + return gmx::c_boltz * T; } else { - return BOLTZ * T * (0.5 * x - 1 + x / (ex - 1)); + return gmx::c_boltz * T * (0.5 * x - 1 + x / (ex - 1)); } } @@ -287,14 +287,14 @@ static real* allocateEigenvectors(int nrow, int first, int last, bool ignoreBegi */ static double calcTranslationalHeatCapacity() { - return RGAS * 1.5; + return gmx::c_universalGasConstant * 1.5; } /*! \brief Compute internal energy due to translational motion */ static double calcTranslationalInternalEnergy(double T) { - return BOLTZ * T * 1.5; + return gmx::c_boltz * T * 1.5; } /*! \brief Compute heat capacity due to rotational motion @@ -307,11 +307,11 @@ static double calcRotationalInternalEnergy(gmx_bool linear, double T) { if (linear) { - return BOLTZ * T; + return gmx::c_boltz * T; } else { - return BOLTZ * T * 1.5; + return gmx::c_boltz * T * 1.5; } } @@ -324,11 +324,11 @@ static double calcRotationalHeatCapacity(gmx_bool linear) { if (linear) { - return RGAS; + return gmx::c_universalGasConstant; } else { - return RGAS * 1.5; + return gmx::c_universalGasConstant * 1.5; } } @@ -361,9 +361,9 @@ static void analyzeThermochemistry(FILE* fp, principal_comp(index.size(), index.data(), top.atoms.atom, as_rvec_array(x_com.data()), trans, inertia); bool linear = (inertia[XX] / inertia[YY] < linear_toler && inertia[XX] / inertia[ZZ] < linear_toler); // (kJ/mol ps)^2/(Dalton nm^2 kJ/mol K) = - // KILO kg m^2 ps^2/(s^2 mol g/mol nm^2 K) = - // KILO^2 10^18 / 10^24 K = 1/K - double rot_const = gmx::square(PLANCK) / (8 * gmx::square(M_PI) * BOLTZ); + // c_kilo kg m^2 ps^2/(s^2 mol g/mol nm^2 K) = + // c_kilo^2 10^18 / 10^24 K = 1/K + double rot_const = gmx::square(gmx::c_planck) / (8 * gmx::square(M_PI) * gmx::c_boltz); // Rotational temperature (1/K) rvec theta = { 0, 0, 0 }; if (linear) @@ -729,8 +729,8 @@ int gmx_nmeig(int argc, char* argv[]) * light. Do this by first converting to omega^2 (units 1/s), take the square * root, and finally divide by the speed of light (nm/ps in gromacs). */ - factor_gmx_to_omega2 = 1.0E21 / (AVOGADRO * AMU); - factor_omega_to_wavenumber = 1.0E-5 / (2.0 * M_PI * SPEED_OF_LIGHT); + factor_gmx_to_omega2 = 1.0E21 / (gmx::c_avogadro * gmx::c_amu); + factor_omega_to_wavenumber = 1.0E-5 / (2.0 * M_PI * gmx::c_speedOfLight); value = 0; for (i = begin; (i <= end); i++) @@ -758,8 +758,8 @@ int gmx_nmeig(int argc, char* argv[]) qu = u_corr(nu, T); if (i > end - nharm) { - qcv += BOLTZ * KILO; - qu += BOLTZ * T; + qcv += gmx::c_boltz * gmx::c_kilo; + qu += gmx::c_boltz * T; } fprintf(qc, "%6d %15g %15g\n", i, qcv, qu); qcvtot += qcv; diff --git a/src/gromacs/gmxana/gmx_nmens.cpp b/src/gromacs/gmxana/gmx_nmens.cpp index 2ddcc6a234..63ae3e5226 100644 --- a/src/gromacs/gmxana/gmx_nmens.cpp +++ b/src/gromacs/gmxana/gmx_nmens.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -232,7 +232,7 @@ int gmx_nmens(int argc, char* argv[]) /* (r-0.5) n times: var_n = n * var_1 = n/12 n=4: var_n = 1/3, so multiply with 3 */ - rfac = std::sqrt(3.0 * BOLTZ * temp / eigval[iout[j]]); + rfac = std::sqrt(3.0 * gmx::c_boltz * temp / eigval[iout[j]]); rhalf = 2.0 * rfac; rfac = rfac / im; diff --git a/src/gromacs/gmxana/gmx_nmtraj.cpp b/src/gromacs/gmxana/gmx_nmtraj.cpp index 006f66b2f9..1ee9de463e 100644 --- a/src/gromacs/gmxana/gmx_nmtraj.cpp +++ b/src/gromacs/gmxana/gmx_nmtraj.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2017,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2017,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -231,7 +231,7 @@ int gmx_nmtraj(int argc, char* argv[]) /* Derive amplitude from temperature and eigenvalue if we can */ /* Convert eigenvalue to angular frequency, in units s^(-1) */ - omega = std::sqrt(eigval[kmode] * 1.0E21 / (AVOGADRO * AMU)); + omega = std::sqrt(eigval[kmode] * 1.0E21 / (gmx::c_avogadro * gmx::c_amu)); /* Harmonic motion will be x=x0 + A*sin(omega*t)*eigenvec. * The velocity is thus: * @@ -261,10 +261,10 @@ int gmx_nmtraj(int argc, char* argv[]) /* Convert Ekin from amu*(nm/s)^2 to J, i.e., kg*(m/s)^2 * This will also be proportional to A^2 */ - Ekin *= AMU * 1E-18; + Ekin *= gmx::c_amu * 1E-18; /* Set the amplitude so the energy is kT/2 */ - amplitude[i] = std::sqrt(0.5 * BOLTZMANN * temp / Ekin); + amplitude[i] = std::sqrt(0.5 * gmx::c_boltzmann * temp / Ekin); } else { diff --git a/src/gromacs/gmxana/gmx_rama.cpp b/src/gromacs/gmxana/gmx_rama.cpp index c10b877468..24c554cb8c 100644 --- a/src/gromacs/gmxana/gmx_rama.cpp +++ b/src/gromacs/gmxana/gmx_rama.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2017,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2017,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -58,8 +58,8 @@ static void plot_rama(FILE* out, t_xrama* xr) for (i = 0; (i < xr->npp); i++) { - phi = xr->dih[xr->pp[i].iphi].ang * RAD2DEG; - psi = xr->dih[xr->pp[i].ipsi].ang * RAD2DEG; + phi = xr->dih[xr->pp[i].iphi].ang * gmx::c_rad2Deg; + psi = xr->dih[xr->pp[i].ipsi].ang * gmx::c_rad2Deg; fprintf(out, "%g %g %s\n", phi, psi, xr->pp[i].label); } } diff --git a/src/gromacs/gmxana/gmx_sham.cpp b/src/gromacs/gmxana/gmx_sham.cpp index 2fe286abe1..44dffff46d 100644 --- a/src/gromacs/gmxana/gmx_sham.cpp +++ b/src/gromacs/gmxana/gmx_sham.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -52,6 +52,7 @@ #include "gromacs/gmxana/gstat.h" #include "gromacs/math/functions.h" #include "gromacs/math/units.h" +#include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/topology/topology.h" #include "gromacs/utility/arraysize.h" @@ -528,7 +529,7 @@ static void do_sham(const char* fn, bfac[i] = ibox[i] / (max_eig[i] - min_eig[i]); } /* Do the binning */ - bref = 1 / (BOLTZ * Tref); + bref = 1 / (gmx::c_boltz * Tref); snew(bE, n); if (bGE || nenerT == 2) { @@ -542,7 +543,7 @@ static void do_sham(const char* fn, } else { - bE[j] = (bref - 1 / (BOLTZ * enerT[1][j])) * enerT[0][j]; + bE[j] = (bref - 1 / (gmx::c_boltz * enerT[1][j])) * enerT[0][j]; } Emin = std::min(Emin, static_cast(bE[j])); } @@ -605,7 +606,7 @@ static void do_sham(const char* fn, } else if (idim[i] == -1) { - efac /= std::sin(DEG2RAD * eig[i][j]); + efac /= std::sin(gmx::c_deg2Rad * eig[i][j]); } } /* Update the probability */ @@ -635,7 +636,7 @@ static void do_sham(const char* fn, if (P[i] != 0) { Pmax = std::max(P[i], Pmax); - W[i] = -BOLTZ * Tref * std::log(P[i]); + W[i] = -gmx::c_boltz * Tref * std::log(P[i]); if (W[i] < Wmin) { Wmin = W[i]; diff --git a/src/gromacs/gmxana/gmx_sigeps.cpp b/src/gromacs/gmxana/gmx_sigeps.cpp index e50ca4ee42..6b102699f4 100644 --- a/src/gromacs/gmxana/gmx_sigeps.cpp +++ b/src/gromacs/gmxana/gmx_sigeps.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -53,7 +53,7 @@ static real pot(real x, real qq, real c6, real cn, int npow) { - return cn * pow(x, -npow) - c6 / gmx::power6(x) + qq * ONE_4PI_EPS0 / x; + return cn * pow(x, -npow) - c6 / gmx::power6(x) + qq * gmx::c_one4PiEps0 / x; } static real bhpot(real x, real A, real B, real C) @@ -64,7 +64,7 @@ static real bhpot(real x, real A, real B, real C) static real dpot(real x, real qq, real c6, real cn, int npow) { return -(npow * cn * std::pow(x, -npow - 1) - 6 * c6 / (x * gmx::power6(x)) - + qq * ONE_4PI_EPS0 / gmx::square(x)); + + qq * gmx::c_one4PiEps0 / gmx::square(x)); } int gmx_sigeps(int argc, char* argv[]) diff --git a/src/gromacs/gmxana/gmx_spol.cpp b/src/gromacs/gmxana/gmx_spol.cpp index a787a78aeb..d4321ed96f 100644 --- a/src/gromacs/gmxana/gmx_spol.cpp +++ b/src/gromacs/gmxana/gmx_spol.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -335,7 +335,7 @@ int gmx_spol(int argc, char* argv[]) } unitv(dir, dir); - svmul(ENM2DEBYE, dip, dip); + svmul(gmx::c_enm2Debye, dip, dip); dip2 = norm2(dip); sdip += std::sqrt(dip2); sdip2 += dip2; diff --git a/src/gromacs/gmxana/gmx_tcaf.cpp b/src/gromacs/gmxana/gmx_tcaf.cpp index a0d0ce2ae6..92f543e294 100644 --- a/src/gromacs/gmxana/gmx_tcaf.cpp +++ b/src/gromacs/gmxana/gmx_tcaf.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -241,7 +241,8 @@ static void process_tcaf(int nframes, fitparms[0] = 1; fitparms[1] = 1; do_lmfit(ncorr, tcaf[k], sig, dt, nullptr, 0, ncorr * dt, oenv, bDebugMode(), effnVAC, fitparms, 0, nullptr); - eta = 1000 * fitparms[1] * rho / (4 * fitparms[0] * PICO * norm2(kfac[k]) / (NANO * NANO)); + eta = 1000 * fitparms[1] * rho + / (4 * fitparms[0] * gmx::c_pico * norm2(kfac[k]) / (gmx::c_nano * gmx::c_nano)); fprintf(stdout, "k %6.3f tau %6.3f eta %8.5f 10^-3 kg/(m s)\n", norm(kfac[k]), fitparms[0], eta); fprintf(fp_vk, "%6.3f %g\n", norm(kfac[k]), eta); for (i = 0; i < ncorr; i++) @@ -264,7 +265,7 @@ static void process_tcaf(int nframes, fitparms[1] = 1; do_lmfit(ncorr, tcafc[k], sig, dt, nullptr, 0, ncorr * dt, oenv, bDebugMode(), effnVAC, fitparms, 0, nullptr); eta = 1000 * fitparms[1] * rho - / (4 * fitparms[0] * PICO * norm2(kfac[kset_c[k]]) / (NANO * NANO)); + / (4 * fitparms[0] * gmx::c_pico * norm2(kfac[kset_c[k]]) / (gmx::c_nano * gmx::c_nano)); fprintf(stdout, "k %6.3f tau %6.3f Omega %6.3f eta %8.5f 10^-3 kg/(m s)\n", norm(kfac[kset_c[k]]), @@ -508,7 +509,7 @@ int gmx_tcaf(int argc, char* argv[]) dt = (t1 - t0) / (nframes - 1); - rho *= sysmass / nframes * AMU / (NANO * NANO * NANO); + rho *= sysmass / nframes * gmx::c_amu / (gmx::c_nano * gmx::c_nano * gmx::c_nano); fprintf(stdout, "Density = %g (kg/m^3)\n", rho); process_tcaf(nframes, dt, diff --git a/src/gromacs/gmxana/gmx_traj.cpp b/src/gromacs/gmxana/gmx_traj.cpp index 53c71154a0..6526ccd7a2 100644 --- a/src/gromacs/gmxana/gmx_traj.cpp +++ b/src/gromacs/gmxana/gmx_traj.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -389,7 +389,7 @@ static real temp(rvec v[], const real mass[], int isize, const int index[]) ekin2 += mass[j] * norm2(v[j]); } - return ekin2 / (3 * isize * BOLTZ); + return ekin2 / (3 * isize * gmx::c_boltz); } static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[]) diff --git a/src/gromacs/gmxana/gmx_velacc.cpp b/src/gromacs/gmxana/gmx_velacc.cpp index 20c21c5b7a..b12cfbb021 100644 --- a/src/gromacs/gmxana/gmx_velacc.cpp +++ b/src/gromacs/gmxana/gmx_velacc.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -161,7 +161,7 @@ static void calc_spectrum(int n, const real c[], real dt, const char* fn, gmx_ou * The timestep between saving the trajectory is * 1e7 is to convert nanometer to cm */ - recip_fac = bRecip ? (1e7 / SPEED_OF_LIGHT) : 1.0; + recip_fac = bRecip ? (1e7 / gmx::c_speedOfLight) : 1.0; for (i = 0; (i < n); i += 2) { nu = i / (2 * dt); diff --git a/src/gromacs/gmxana/gmx_wham.cpp b/src/gromacs/gmxana/gmx_wham.cpp index 04d3106b36..9f263dc3a0 100644 --- a/src/gromacs/gmxana/gmx_wham.cpp +++ b/src/gromacs/gmxana/gmx_wham.cpp @@ -580,8 +580,8 @@ static void setup_acc_wham(const double* profile, t_UmbrellaWindow* window, int } } /* Note: there are two contributions to bin k in the wham equations: - i) N[j]*exp(- U/(BOLTZ*opt->Temperature) + window[i].z[j]) - ii) exp(- U/(BOLTZ*opt->Temperature)) + i) N[j]*exp(- U/(c_boltz*opt->Temperature) + window[i].z[j]) + ii) exp(- U/(c_boltz*opt->Temperature)) where U is the umbrella potential If any of these number is larger wham_contrib_lim, I set contrib=TRUE */ @@ -594,8 +594,9 @@ static void setup_acc_wham(const double* profile, t_UmbrellaWindow* window, int { U = tabulated_pot(distance, opt); /* Use tabulated potential */ } - contrib1 = profile[k] * std::exp(-U / (BOLTZ * opt->Temperature)); - contrib2 = window[i].N[j] * std::exp(-U / (BOLTZ * opt->Temperature) + window[i].z[j]); + contrib1 = profile[k] * std::exp(-U / (gmx::c_boltz * opt->Temperature)); + contrib2 = window[i].N[j] + * std::exp(-U / (gmx::c_boltz * opt->Temperature) + window[i].z[j]); window[i].bContrib[j][k] = (contrib1 > wham_contrib_lim || contrib2 > wham_contrib_lim); bAnyContrib = bAnyContrib || window[i].bContrib[j][k]; if (window[i].bContrib[j][k]) @@ -689,7 +690,7 @@ static void calc_profile(double* profile, t_UmbrellaWindow* window, int nWindows U = tabulated_pot(distance, opt); /* Use tabulated potential */ } denom += invg * window[j].N[k] - * std::exp(-U / (BOLTZ * opt->Temperature) + window[j].z[k]); + * std::exp(-U / (gmx::c_boltz * opt->Temperature) + window[j].z[k]); } } profile[i] = num / denom; @@ -755,7 +756,7 @@ static double calc_z(const double* profile, t_UmbrellaWindow* window, int nWindo { U = tabulated_pot(distance, opt); /* Use tabulated potential */ } - total += profile[k] * std::exp(-U / (BOLTZ * opt->Temperature)); + total += profile[k] * std::exp(-U / (gmx::c_boltz * opt->Temperature)); } /* Avoid floating point exception if window is far outside min and max */ if (total != 0.0) @@ -852,11 +853,11 @@ static void prof_normalization_and_unit(double* profile, t_UmbrellaOptions* opt) } else if (opt->unit == en_kJ) { - unit_factor = BOLTZ * opt->Temperature; + unit_factor = gmx::c_boltz * opt->Temperature; } else if (opt->unit == en_kCal) { - unit_factor = (BOLTZ / CAL2JOULE) * opt->Temperature; + unit_factor = (gmx::c_boltz / gmx::c_cal2Joule) * opt->Temperature; } else { @@ -2711,7 +2712,7 @@ static void guessPotByIntegration(t_UmbrellaWindow* window, int nWindows, t_Umbr */ for (j = 0; j < opt->bins; ++j) { - pot[j] = std::exp(-pot[j] / (BOLTZ * opt->Temperature)); + pot[j] = std::exp(-pot[j] / (gmx::c_boltz * opt->Temperature)); } calc_z(pot, window, nWindows, opt, TRUE); diff --git a/src/gromacs/gmxana/hxprops.cpp b/src/gromacs/gmxana/hxprops.cpp index 66e4ddb135..a51457a369 100644 --- a/src/gromacs/gmxana/hxprops.cpp +++ b/src/gromacs/gmxana/hxprops.cpp @@ -143,7 +143,7 @@ static real rot(rvec x1, const rvec x2) xx = cp * x2[XX] + sp * x2[YY]; yy = -sp * x2[XX] + cp * x2[YY]; - dphi = RAD2DEG * std::atan2(yy, xx); + dphi = gmx::c_rad2Deg * std::atan2(yy, xx); return dphi; } @@ -190,7 +190,7 @@ real ca_phi(int gnx, const int index[], rvec x[]) aj = index[i + 1]; ak = index[i + 2]; al = index[i + 3]; - phi = RAD2DEG + phi = gmx::c_rad2Deg * dih_angle(x[ai], x[aj], x[ak], x[al], nullptr, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3); phitot += phi; } @@ -497,7 +497,7 @@ void calc_hxprops(int nres, t_bb bb[], const rvec x[]) bb[i].d5 = norm(dx); } - bb[i].phi = RAD2DEG + bb[i].phi = gmx::c_rad2Deg * dih_angle(x[bb[i].Cprev], x[bb[i].N], x[bb[i].CA], @@ -511,7 +511,7 @@ void calc_hxprops(int nres, t_bb bb[], const rvec x[]) &t1, &t2, &t3); - bb[i].psi = RAD2DEG + bb[i].psi = gmx::c_rad2Deg * dih_angle(x[bb[i].N], x[bb[i].CA], x[bb[i].C], @@ -527,9 +527,9 @@ void calc_hxprops(int nres, t_bb bb[], const rvec x[]) &t3); bb[i].pprms2 = gmx::square(bb[i].phi - PHI_AHX) + gmx::square(bb[i].psi - PSI_AHX); - bb[i].jcaha += 1.4 * std::sin((bb[i].psi + 138.0) * DEG2RAD) - - 4.1 * std::cos(2.0 * DEG2RAD * (bb[i].psi + 138.0)) - + 2.0 * std::cos(2.0 * DEG2RAD * (bb[i].phi + 30.0)); + bb[i].jcaha += 1.4 * std::sin((bb[i].psi + 138.0) * gmx::c_deg2Rad) + - 4.1 * std::cos(2.0 * gmx::c_deg2Rad * (bb[i].psi + 138.0)) + + 2.0 * std::cos(2.0 * gmx::c_deg2Rad * (bb[i].phi + 30.0)); } } diff --git a/src/gromacs/gmxana/pp2shift.cpp b/src/gromacs/gmxana/pp2shift.cpp index 707673c662..da8ae11dba 100644 --- a/src/gromacs/gmxana/pp2shift.cpp +++ b/src/gromacs/gmxana/pp2shift.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -119,13 +119,13 @@ static void dump_sd(const char* fn, t_shiftdata* sd) { snew(newdata[i], nny); phi = i * 2 * M_PI / (nnx - 1); - x_phi[i] = phi * RAD2DEG - 180; + x_phi[i] = phi * gmx::c_rad2Deg - 180; for (j = 0; (j < nny); j++) { psi = j * 2 * M_PI / (nny - 1); if (i == 0) { - y_psi[j] = psi * RAD2DEG - 180; + y_psi[j] = psi * gmx::c_rad2Deg - 180; } /*if (((i % nfac) == 0) && ((j % nfac) == 0)) newdata[i][j] = sd->data[i/nfac][j/nfac]; diff --git a/src/gromacs/gmxana/thermochemistry.cpp b/src/gromacs/gmxana/thermochemistry.cpp index 20bae48db8..2b5d68988f 100644 --- a/src/gromacs/gmxana/thermochemistry.cpp +++ b/src/gromacs/gmxana/thermochemistry.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -46,14 +46,14 @@ static double eigval_to_frequency(double eigval) { - double factor_gmx_to_omega2 = 1.0E21 / (AVOGADRO * AMU); + double factor_gmx_to_omega2 = 1.0E21 / (gmx::c_avogadro * gmx::c_amu); return std::sqrt(std::max(0.0, eigval) * factor_gmx_to_omega2); } double calcZeroPointEnergy(gmx::ArrayRef eigval, real scale_factor) { // Convert frequency (ps^-1) to energy (kJ/mol) - double factor = PLANCK * PICO / (2.0 * M_PI); + double factor = gmx::c_planck * gmx::c_pico / (2.0 * M_PI); double zpe = 0; for (auto& r : eigval) { @@ -67,13 +67,13 @@ double calcVibrationalInternalEnergy(gmx::ArrayRef eigval, real temp { size_t nskip = linear ? 5 : 6; double Evib = 0; - double hbar = PLANCK1 / (2 * M_PI); + double hbar = gmx::c_planck1 / (2 * M_PI); for (gmx::index i = nskip; i < eigval.ssize(); i++) { if (eigval[i] > 0) { double omega = scale_factor * eigval_to_frequency(eigval[i]); - double hwkT = (hbar * omega) / (BOLTZMANN * temperature); + double hwkT = (hbar * omega) / (gmx::c_boltzmann * temperature); // Prevent overflow by checking for unreasonably large numbers. if (hwkT < 100) { @@ -92,20 +92,20 @@ double calcVibrationalInternalEnergy(gmx::ArrayRef eigval, real temp } } } - return temperature * BOLTZ * Evib; + return temperature * gmx::c_boltz * Evib; } double calcVibrationalHeatCapacity(gmx::ArrayRef eigval, real temperature, gmx_bool linear, real scale_factor) { size_t nskip = linear ? 5 : 6; double cv = 0; - double hbar = PLANCK1 / (2 * M_PI); + double hbar = gmx::c_planck1 / (2 * M_PI); for (gmx::index i = nskip; i < eigval.ssize(); i++) { if (eigval[i] > 0) { double omega = scale_factor * eigval_to_frequency(eigval[i]); - double hwkT = (hbar * omega) / (BOLTZMANN * temperature); + double hwkT = (hbar * omega) / (gmx::c_boltzmann * temperature); // Prevent overflow by checking for unreasonably large numbers. if (hwkT < 100) { @@ -124,20 +124,21 @@ double calcVibrationalHeatCapacity(gmx::ArrayRef eigval, real temper } } } - return RGAS * cv; + return gmx::c_universalGasConstant * cv; } double calcTranslationalEntropy(real mass, real temperature, real pressure) { - double kT = BOLTZ * temperature; + double kT = gmx::c_boltz * temperature; GMX_RELEASE_ASSERT(mass > 0, "Molecular mass should be larger than zero"); GMX_RELEASE_ASSERT(pressure > 0, "Pressure should be larger than zero"); GMX_RELEASE_ASSERT(temperature > 0, "Temperature should be larger than zero"); // Convert bar to Pascal - double P = pressure * 1e5; - double qT = (std::pow(2 * M_PI * mass * kT / gmx::square(PLANCK), 1.5) * (kT / P) * (1e30 / AVOGADRO)); - return RGAS * (std::log(qT) + 2.5); + double P = pressure * 1e5; + double qT = (std::pow(2 * M_PI * mass * kT / gmx::square(gmx::c_planck), 1.5) * (kT / P) + * (1e30 / gmx::c_avogadro)); + return gmx::c_universalGasConstant * (std::log(qT) + 2.5); } double calcRotationalEntropy(real temperature, int natom, gmx_bool linear, const rvec theta, real sigma_r) @@ -152,14 +153,14 @@ double calcRotationalEntropy(real temperature, int natom, gmx_bool linear, const { GMX_RELEASE_ASSERT(theta[0] > 0, "Theta should be larger than zero"); double qR = temperature / (sigma_r * theta[0]); - sR = RGAS * (std::log(qR) + 1); + sR = gmx::c_universalGasConstant * (std::log(qR) + 1); } else { double Q = theta[XX] * theta[YY] * theta[ZZ]; GMX_RELEASE_ASSERT(Q > 0, "Q should be larger than zero"); double qR = std::sqrt(M_PI * std::pow(temperature, 3) / Q) / sigma_r; - sR = RGAS * (std::log(qR) + 1.5); + sR = gmx::c_universalGasConstant * (std::log(qR) + 1.5); } } return sR; @@ -169,13 +170,13 @@ double calcQuasiHarmonicEntropy(gmx::ArrayRef eigval, real temperatu { size_t nskip = bLinear ? 5 : 6; double S = 0; - double hbar = PLANCK1 / (2 * M_PI); + double hbar = gmx::c_planck1 / (2 * M_PI); for (gmx::index i = nskip; (i < eigval.ssize()); i++) { if (eigval[i] > 0) { double omega = scale_factor * eigval_to_frequency(eigval[i]); - double hwkT = (hbar * omega) / (BOLTZMANN * temperature); + double hwkT = (hbar * omega) / (gmx::c_boltzmann * temperature); double dS = (hwkT / std::expm1(hwkT) - std::log1p(-std::exp(-hwkT))); S += dS; if (debug) @@ -194,16 +195,16 @@ double calcQuasiHarmonicEntropy(gmx::ArrayRef eigval, real temperatu fprintf(debug, "eigval[%d] = %g\n", static_cast(i + 1), static_cast(eigval[i])); } } - return S * RGAS; + return S * gmx::c_universalGasConstant; } double calcSchlitterEntropy(gmx::ArrayRef eigval, real temperature, gmx_bool bLinear) { size_t nskip = bLinear ? 5 : 6; - double hbar = PLANCK1 / (2 * M_PI); // J s - double kt = BOLTZMANN * temperature; // J + double hbar = gmx::c_planck1 / (2 * M_PI); // J s + double kt = gmx::c_boltzmann * temperature; // J double kteh = kt * std::exp(2.0) / (hbar * hbar); // 1/(J s^2) = 1/(kg m^2) - double evcorr = NANO * NANO * AMU; + double evcorr = gmx::c_nano * gmx::c_nano * gmx::c_amu; if (debug) { fprintf(debug, "n = %td, kteh = %g evcorr = %g\n", ssize(eigval), kteh, evcorr); @@ -214,5 +215,5 @@ double calcSchlitterEntropy(gmx::ArrayRef eigval, real temperature, double dd = 1 + kteh * eigval[i] * evcorr; deter += std::log(dd); } - return 0.5 * RGAS * deter; + return 0.5 * gmx::c_universalGasConstant * deter; } diff --git a/src/gromacs/gmxpreprocess/calch.cpp b/src/gromacs/gmxpreprocess/calch.cpp index 24572cd44b..18a2043ba2 100644 --- a/src/gromacs/gmxpreprocess/calch.cpp +++ b/src/gromacs/gmxpreprocess/calch.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2010,2014,2015,2016,2019, by the GROMACS development team, led by + * Copyright (c) 2010,2014,2015,2016,2019,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -97,9 +97,9 @@ void calc_h_pos(int nht, rvec xa[], rvec xh[], int* l) #define alfaHpl (2 * M_PI / 3) /* 120 degrees */ #define distH 0.1 -#define alfaCOM (DEG2RAD * 117) -#define alfaCO (DEG2RAD * 121) -#define alfaCOA (DEG2RAD * 115) +#define alfaCOM (gmx::c_deg2Rad * 117) +#define alfaCO (gmx::c_deg2Rad * 121) +#define alfaCOA (gmx::c_deg2Rad * 115) #define distO 0.123 #define distOA 0.125 diff --git a/src/gromacs/gmxpreprocess/convparm.cpp b/src/gromacs/gmxpreprocess/convparm.cpp index f06620e5b4..1712ec3c16 100644 --- a/src/gromacs/gmxpreprocess/convparm.cpp +++ b/src/gromacs/gmxpreprocess/convparm.cpp @@ -145,9 +145,9 @@ static int assign_param(t_functype ftype, { case F_G96ANGLES: /* Post processing of input data: store cosine iso angle itself */ - newparam->harmonic.rA = cos(old[0] * DEG2RAD); + newparam->harmonic.rA = cos(old[0] * gmx::c_deg2Rad); newparam->harmonic.krA = old[1]; - newparam->harmonic.rB = cos(old[2] * DEG2RAD); + newparam->harmonic.rB = cos(old[2] * gmx::c_deg2Rad); newparam->harmonic.krB = old[3]; break; case F_G96BONDS: @@ -430,8 +430,8 @@ static int assign_param(t_functype ftype, newparam->vsite.f = old[5]; break; case F_VSITE3FAD: - newparam->vsite.a = old[1] * cos(DEG2RAD * old[0]); - newparam->vsite.b = old[1] * sin(DEG2RAD * old[0]); + newparam->vsite.a = old[1] * cos(gmx::c_deg2Rad * old[0]); + newparam->vsite.b = old[1] * sin(gmx::c_deg2Rad * old[0]); newparam->vsite.c = old[2]; newparam->vsite.d = old[3]; newparam->vsite.e = old[4]; diff --git a/src/gromacs/gmxpreprocess/editconf.cpp b/src/gromacs/gmxpreprocess/editconf.cpp index 085907a90b..7a64f7abdc 100644 --- a/src/gromacs/gmxpreprocess/editconf.cpp +++ b/src/gromacs/gmxpreprocess/editconf.cpp @@ -970,9 +970,9 @@ int gmx_editconf(int argc, char* argv[]) printf(" center :%7.3f%7.3f%7.3f (nm)\n", gc[XX], gc[YY], gc[ZZ]); printf(" box vectors :%7.3f%7.3f%7.3f (nm)\n", norm(box[XX]), norm(box[YY]), norm(box[ZZ])); printf(" box angles :%7.2f%7.2f%7.2f (degrees)\n", - norm2(box[ZZ]) == 0 ? 0 : RAD2DEG * gmx_angle(box[YY], box[ZZ]), - norm2(box[ZZ]) == 0 ? 0 : RAD2DEG * gmx_angle(box[XX], box[ZZ]), - norm2(box[YY]) == 0 ? 0 : RAD2DEG * gmx_angle(box[XX], box[YY])); + norm2(box[ZZ]) == 0 ? 0 : gmx::c_rad2Deg * gmx_angle(box[YY], box[ZZ]), + norm2(box[ZZ]) == 0 ? 0 : gmx::c_rad2Deg * gmx_angle(box[XX], box[ZZ]), + norm2(box[YY]) == 0 ? 0 : gmx::c_rad2Deg * gmx_angle(box[XX], box[YY])); printf(" box volume :%7.2f (nm^3)\n", det(box)); } @@ -1005,7 +1005,7 @@ int gmx_editconf(int argc, char* argv[]) real vol, dens; vol = det(box); - dens = (mass * AMU) / (vol * NANO * NANO * NANO); + dens = (mass * gmx::c_amu) / (vol * gmx::c_nano * gmx::c_nano * gmx::c_nano); fprintf(stderr, "Volume of input %g (nm^3)\n", vol); fprintf(stderr, "Mass of input %g (a.m.u.)\n", mass); fprintf(stderr, "Density of input %g (g/l)\n", dens); @@ -1125,7 +1125,7 @@ int gmx_editconf(int argc, char* argv[]) rotangles[ZZ]); for (i = 0; i < DIM; i++) { - rotangles[i] *= DEG2RAD; + rotangles[i] *= gmx::c_deg2Rad; } rotate_conf(natom, x, v, rotangles[XX], rotangles[YY], rotangles[ZZ]); } @@ -1235,9 +1235,9 @@ int gmx_editconf(int argc, char* argv[]) { printf("new box vectors :%7.3f%7.3f%7.3f (nm)\n", norm(box[XX]), norm(box[YY]), norm(box[ZZ])); printf("new box angles :%7.2f%7.2f%7.2f (degrees)\n", - norm2(box[ZZ]) == 0 ? 0 : RAD2DEG * gmx_angle(box[YY], box[ZZ]), - norm2(box[ZZ]) == 0 ? 0 : RAD2DEG * gmx_angle(box[XX], box[ZZ]), - norm2(box[YY]) == 0 ? 0 : RAD2DEG * gmx_angle(box[XX], box[YY])); + norm2(box[ZZ]) == 0 ? 0 : gmx::c_rad2Deg * gmx_angle(box[YY], box[ZZ]), + norm2(box[ZZ]) == 0 ? 0 : gmx::c_rad2Deg * gmx_angle(box[XX], box[ZZ]), + norm2(box[YY]) == 0 ? 0 : gmx::c_rad2Deg * gmx_angle(box[XX], box[YY])); printf("new box volume :%7.2f (nm^3)\n", det(box)); } diff --git a/src/gromacs/gmxpreprocess/gen_maxwell_velocities.cpp b/src/gromacs/gmxpreprocess/gen_maxwell_velocities.cpp index 2893816f38..b02735d60d 100644 --- a/src/gromacs/gmxpreprocess/gen_maxwell_velocities.cpp +++ b/src/gromacs/gmxpreprocess/gen_maxwell_velocities.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -55,13 +55,11 @@ static void low_mspeed(real tempi, gmx_mtop_t* mtop, rvec v[], gmx::ThreeFry2x64<>* rng, const gmx::MDLogger& logger) { int nrdf; - real boltz; real ekin, temp; gmx::TabulatedNormalDistribution normalDist; - boltz = BOLTZ * tempi; - ekin = 0.0; - nrdf = 0; + ekin = 0.0; + nrdf = 0; for (const AtomProxy atomP : AtomRange(*mtop)) { const t_atom& local = atomP.atom(); @@ -70,7 +68,7 @@ static void low_mspeed(real tempi, gmx_mtop_t* mtop, rvec v[], gmx::ThreeFry2x64 if (mass > 0) { rng->restart(i, 0); - real sd = std::sqrt(boltz / mass); + real sd = std::sqrt(gmx::c_boltz * tempi / mass); for (int m = 0; (m < DIM); m++) { v[i][m] = sd * normalDist(*rng); @@ -79,7 +77,7 @@ static void low_mspeed(real tempi, gmx_mtop_t* mtop, rvec v[], gmx::ThreeFry2x64 nrdf += DIM; } } - temp = (2.0 * ekin) / (nrdf * BOLTZ); + temp = (2.0 * ekin) / (nrdf * gmx::c_boltz); if (temp > 0) { real scal = std::sqrt(tempi / temp); diff --git a/src/gromacs/gmxpreprocess/gen_vsite.cpp b/src/gromacs/gmxpreprocess/gen_vsite.cpp index 8c6ac3352e..201a412db0 100644 --- a/src/gromacs/gmxpreprocess/gen_vsite.cpp +++ b/src/gromacs/gmxpreprocess/gen_vsite.cpp @@ -754,7 +754,7 @@ static void add_vsites(gmx::ArrayRef plist, } /* for i */ } -#define ANGLE_6RING (DEG2RAD * 120) +#define ANGLE_6RING (gmx::c_deg2Rad * 120) /* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */ /* get a^2 when a, b and alpha are given: */ @@ -1019,21 +1019,21 @@ static int gen_vsites_trp(PreprocessingAtomTypes* atype, b_CE3_HE3 = get_ddb_bond(vsitetop, "TRP", "CE3", "HE3"); b_CZ3_HZ3 = get_ddb_bond(vsitetop, "TRP", "CZ3", "HZ3"); - a_NE1_CE2_CD2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "NE1", "CE2", "CD2"); - a_CE2_CD2_CG = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CG"); - a_CB_CG_CD2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CB", "CG", "CD2"); - a_CD2_CG_CD1 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CG", "CD1"); - - a_CE2_CD2_CE3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CE3"); - a_CD2_CE2_CZ2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CE2", "CZ2"); - a_CD2_CE3_CZ3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "CZ3"); - a_CE3_CZ3_HZ3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE3", "CZ3", "HZ3"); - a_CZ2_CH2_HH2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CZ2", "CH2", "HH2"); - a_CE2_CZ2_HZ2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "HZ2"); - a_CE2_CZ2_CH2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "CH2"); - a_CG_CD1_HD1 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CG", "CD1", "HD1"); - a_HE1_NE1_CE2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "HE1", "NE1", "CE2"); - a_CD2_CE3_HE3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "HE3"); + a_NE1_CE2_CD2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "NE1", "CE2", "CD2"); + a_CE2_CD2_CG = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CG"); + a_CB_CG_CD2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CB", "CG", "CD2"); + a_CD2_CG_CD1 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CD2", "CG", "CD1"); + + a_CE2_CD2_CE3 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CE3"); + a_CD2_CE2_CZ2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CD2", "CE2", "CZ2"); + a_CD2_CE3_CZ3 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "CZ3"); + a_CE3_CZ3_HZ3 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CE3", "CZ3", "HZ3"); + a_CZ2_CH2_HH2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CZ2", "CH2", "HH2"); + a_CE2_CZ2_HZ2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "HZ2"); + a_CE2_CZ2_CH2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "CH2"); + a_CG_CD1_HD1 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CG", "CD1", "HD1"); + a_HE1_NE1_CE2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "HE1", "NE1", "CE2"); + a_CD2_CE3_HE3 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "HE3"); /* Calculate local coordinates. * y-axis (x=0) is the bond CD2-CE2. @@ -1291,7 +1291,7 @@ static int gen_vsites_tyr(PreprocessingAtomTypes* atype, bond_ch = get_ddb_bond(vsitetop, "TYR", "CD1", "HD1"); bond_co = get_ddb_bond(vsitetop, "TYR", "CZ", "OH"); bond_oh = get_ddb_bond(vsitetop, "TYR", "OH", "HH"); - angle_coh = DEG2RAD * get_ddb_angle(vsitetop, "TYR", "CZ", "OH", "HH"); + angle_coh = gmx::c_deg2Rad * get_ddb_angle(vsitetop, "TYR", "CZ", "OH", "HH"); xi[atCG] = -bond_cc + bond_cc * std::cos(ANGLE_6RING); xi[atCD1] = -bond_cc; @@ -1471,30 +1471,30 @@ static int gen_vsites_his(t_atoms* at, b_CE1_NE2 = get_ddb_bond(vsitetop, resname, "CE1", "NE2"); b_CG_CD2 = get_ddb_bond(vsitetop, resname, "CG", "CD2"); b_CD2_NE2 = get_ddb_bond(vsitetop, resname, "CD2", "NE2"); - a_CG_ND1_CE1 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CG", "ND1", "CE1"); - a_CG_CD2_NE2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CG", "CD2", "NE2"); - a_ND1_CE1_NE2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "ND1", "CE1", "NE2"); - a_CE1_NE2_CD2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CE1", "NE2", "CD2"); + a_CG_ND1_CE1 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, resname, "CG", "ND1", "CE1"); + a_CG_CD2_NE2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, resname, "CG", "CD2", "NE2"); + a_ND1_CE1_NE2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, resname, "ND1", "CE1", "NE2"); + a_CE1_NE2_CD2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, resname, "CE1", "NE2", "CD2"); if (ats[atHD1] != NOTSET) { b_ND1_HD1 = get_ddb_bond(vsitetop, resname, "ND1", "HD1"); - a_CE1_ND1_HD1 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CE1", "ND1", "HD1"); + a_CE1_ND1_HD1 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, resname, "CE1", "ND1", "HD1"); } if (ats[atHE2] != NOTSET) { b_NE2_HE2 = get_ddb_bond(vsitetop, resname, "NE2", "HE2"); - a_CE1_NE2_HE2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CE1", "NE2", "HE2"); + a_CE1_NE2_HE2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, resname, "CE1", "NE2", "HE2"); } if (ats[atHD2] != NOTSET) { b_CD2_HD2 = get_ddb_bond(vsitetop, resname, "CD2", "HD2"); - a_NE2_CD2_HD2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "NE2", "CD2", "HD2"); + a_NE2_CD2_HD2 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, resname, "NE2", "CD2", "HD2"); } if (ats[atHE1] != NOTSET) { b_CE1_HE1 = get_ddb_bond(vsitetop, resname, "CE1", "HE1"); - a_NE2_CE1_HE1 = DEG2RAD * get_ddb_angle(vsitetop, resname, "NE2", "CE1", "HE1"); + a_NE2_CE1_HE1 = gmx::c_deg2Rad * get_ddb_angle(vsitetop, resname, "NE2", "CE1", "HE1"); } /* constraints between CG, CE1 and NE1 */ diff --git a/src/gromacs/gmxpreprocess/genion.cpp b/src/gromacs/gmxpreprocess/genion.cpp index 4d98124057..049921fb46 100644 --- a/src/gromacs/gmxpreprocess/genion.cpp +++ b/src/gromacs/gmxpreprocess/genion.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -540,7 +540,7 @@ int gmx_genion(int argc, char* argv[]) { /* Compute number of ions to be added */ vol = det(box); - nsalt = gmx::roundToInt(conc * vol * AVOGADRO / 1e24); + nsalt = gmx::roundToInt(conc * vol * gmx::c_avogadro / 1e24); p_num = abs(nsalt * n_q); n_num = abs(nsalt * p_q); } diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index fe22cd5b0f..67abdfc7ec 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -1384,7 +1384,7 @@ static real calc_temp(const gmx_mtop_t* mtop, const t_inputrec* ir, rvec* v) nrdf += ir->opts.nrdf[g]; } - return sum_mv2 / (nrdf * BOLTZ); + return sum_mv2 / (nrdf * gmx::c_boltz); } static real get_max_reference_temp(const t_inputrec* ir, warninp* wi) @@ -2264,7 +2264,7 @@ int gmx_grompp(int argc, char* argv[]) * to be on the safe side with constraints. */ const real totalEnergyDriftPerAtomPerPicosecond = - 2 * BOLTZ * buffer_temp / (ir->nsteps * ir->delta_t); + 2 * gmx::c_boltz * buffer_temp / (ir->nsteps * ir->delta_t); if (ir->verletbuf_tol > 1.1 * driftTolerance * totalEnergyDriftPerAtomPerPicosecond) { diff --git a/src/gromacs/gmxpreprocess/hizzie.cpp b/src/gromacs/gmxpreprocess/hizzie.cpp index 83ba164b66..9278f09b5f 100644 --- a/src/gromacs/gmxpreprocess/hizzie.cpp +++ b/src/gromacs/gmxpreprocess/hizzie.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -145,7 +145,7 @@ static bool chk_hbonds(int i, t_atoms* pdba, rvec x[], const bool ad[], bool hbo d2 = distance2(x[i], x[j]); rvec_sub(x[i], xh, nh); rvec_sub(x[aj], xh, oh); - a = RAD2DEG * acos(cos_angle(nh, oh)); + a = gmx::c_rad2Deg * acos(cos_angle(nh, oh)); if ((d2 < dist2) && (a > angle)) { hbond[i] = TRUE; diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index c115cbeb9c..eaf605d85b 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -57,6 +57,7 @@ #include "gromacs/gmxpreprocess/toputil.h" #include "gromacs/math/functions.h" #include "gromacs/math/units.h" +#include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/calc_verletbuf.h" #include "gromacs/mdrun/mdmodules.h" @@ -4371,7 +4372,7 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_ * of errors. The factor 0.5 is because energy distributes * equally over Ekin and Epot. */ - max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * BOLTZ * T); + max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * gmx::c_boltz * T); if (max_T_error > T_error_warn) { sprintf(warn_buf, diff --git a/src/gromacs/gmxpreprocess/solvate.cpp b/src/gromacs/gmxpreprocess/solvate.cpp index 4d747607a4..62ab86c576 100644 --- a/src/gromacs/gmxpreprocess/solvate.cpp +++ b/src/gromacs/gmxpreprocess/solvate.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -762,7 +762,7 @@ static void update_top(t_atoms* atoms, vol = det(box); fprintf(stderr, "Volume : %10g (nm^3)\n", vol); - fprintf(stderr, "Density : %10g (g/l)\n", (mtot * 1e24) / (AVOGADRO * vol)); + fprintf(stderr, "Density : %10g (g/l)\n", (mtot * 1e24) / (gmx::c_avogadro * vol)); fprintf(stderr, "Number of solvent molecules: %5d \n\n", nsol); /* open topology file and append sol molecules */ diff --git a/src/gromacs/gmxpreprocess/topshake.cpp b/src/gromacs/gmxpreprocess/topshake.cpp index acea8377c5..f26e044453 100644 --- a/src/gromacs/gmxpreprocess/topshake.cpp +++ b/src/gromacs/gmxpreprocess/topshake.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -162,9 +162,9 @@ void make_shake(gmx::ArrayRef plist, t_atoms* atoms, int nsh } if (bFound) { - real param = std::sqrt(b_ij * b_ij + b_jk * b_jk - - 2.0 * b_ij * b_jk - * cos(DEG2RAD * ang->c0())); + real param = std::sqrt( + b_ij * b_ij + b_jk * b_jk + - 2.0 * b_ij * b_jk * cos(gmx::c_deg2Rad * ang->c0())); std::vector forceParm = { param, param }; if (ftype == F_CONNBONDS || ftype_a == F_CONNBONDS) { diff --git a/src/gromacs/gmxpreprocess/vsite_parm.cpp b/src/gromacs/gmxpreprocess/vsite_parm.cpp index 2c4133119b..36a77cd04c 100644 --- a/src/gromacs/gmxpreprocess/vsite_parm.cpp +++ b/src/gromacs/gmxpreprocess/vsite_parm.cpp @@ -376,7 +376,7 @@ static real get_angle(gmx::ArrayRef angles, t_iato if (((ai == ang.ai()) && (aj == ang.aj()) && (ak == ang.ak())) || ((ai == ang.ak()) && (aj == ang.aj()) && (ak == ang.ai()))) { - angle = DEG2RAD * ang.parameterValue(); + angle = gmx::c_deg2Rad * ang.parameterValue(); break; } } @@ -535,7 +535,7 @@ static bool calc_vsite3fad_param(InteractionOfType* vsi bError = (bij == NOTSET) || (aijk == NOTSET); vsite->setForceParameter(1, bij); - vsite->setForceParameter(0, RAD2DEG * aijk); + vsite->setForceParameter(0, gmx::c_rad2Deg * aijk); if (bSwapParity) { @@ -601,8 +601,8 @@ static bool calc_vsite3out_param(PreprocessingAtomTypes* aty dH = bCN - bNH * std::cos(aCNH); rH = bNH * std::sin(aCNH); /* we assume the H's are symmetrically distributed */ - rHx = rH * std::cos(DEG2RAD * 30); - rHy = rH * std::sin(DEG2RAD * 30); + rHx = rH * std::cos(gmx::c_deg2Rad * 30); + rHy = rH * std::sin(gmx::c_deg2Rad * 30); rM = 0.5 * bMM; dM = std::sqrt(gmx::square(bCM) - gmx::square(rM)); a = 0.5 * ((dH / dM) - (rHy / rM)); @@ -682,9 +682,9 @@ static bool calc_vsite4fd_param(InteractionOfType* vsit .appendTextFormatted( "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f", vsite->ai() + 1, - RAD2DEG * aijk, - RAD2DEG * aijl, - RAD2DEG * aijm); + gmx::c_rad2Deg * aijk, + gmx::c_rad2Deg * aijl, + gmx::c_rad2Deg * aijm); gmx_fatal(FARGS, "invalid construction in calc_vsite4fd for atom %d: " "cosakl=%f, cosakm=%f\n", @@ -752,9 +752,9 @@ static bool calc_vsite4fdn_param(InteractionOfType* vsi .appendTextFormatted( "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f", vsite->ai() + 1, - RAD2DEG * aijk, - RAD2DEG * aijl, - RAD2DEG * aijm); + gmx::c_rad2Deg * aijk, + gmx::c_rad2Deg * aijl, + gmx::c_rad2Deg * aijm); gmx_fatal(FARGS, "invalid construction in calc_vsite4fdn for atom %d: " "pl=%f, pm=%f\n", diff --git a/src/gromacs/gmxpreprocess/x2top.cpp b/src/gromacs/gmxpreprocess/x2top.cpp index b439d0ca45..876afdef94 100644 --- a/src/gromacs/gmxpreprocess/x2top.cpp +++ b/src/gromacs/gmxpreprocess/x2top.cpp @@ -270,7 +270,7 @@ static void calc_angles_dihs(InteractionsOfType* ang, InteractionsOfType* dih, c int ai = angle.ai(); int aj = angle.aj(); int ak = angle.ak(); - real th = RAD2DEG + real th = gmx::c_rad2Deg * bond_angle(x[ai], x[aj], x[ak], bPBC ? &pbc : nullptr, r_ij, r_kj, &costh, &t1, &t2); angle.setForceParameter(0, th); } @@ -281,7 +281,7 @@ static void calc_angles_dihs(InteractionsOfType* ang, InteractionsOfType* dih, c int ak = dihedral.ak(); int al = dihedral.al(); real ph = - RAD2DEG + gmx::c_rad2Deg * dih_angle(x[ai], x[aj], x[ak], x[al], bPBC ? &pbc : nullptr, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3); dihedral.setForceParameter(0, ph); } diff --git a/src/gromacs/imd/imd.cpp b/src/gromacs/imd/imd.cpp index cf1a9d81da..aaec29e3ad 100644 --- a/src/gromacs/imd/imd.cpp +++ b/src/gromacs/imd/imd.cpp @@ -570,9 +570,9 @@ static int imd_send_rvecs(IMDSocket* socket, int nat, rvec* x, char* buffer) fill_header(reinterpret_cast(buffer), IMD_FCOORDS, static_cast(nat)); for (i = 0; i < nat; i++) { - sendx[0] = static_cast(x[i][0]) * NM2A; - sendx[1] = static_cast(x[i][1]) * NM2A; - sendx[2] = static_cast(x[i][2]) * NM2A; + sendx[0] = static_cast(x[i][0]) * gmx::c_nm2A; + sendx[1] = static_cast(x[i][1]) * gmx::c_nm2A; + sendx[2] = static_cast(x[i][2]) * gmx::c_nm2A; memcpy(buffer + c_headerSize + i * tuplesize, sendx, tuplesize); } @@ -734,7 +734,7 @@ void ImdSession::Impl::prepareMDForces() void ImdSession::Impl::copyToMDForces() { int i; - real conversion = CAL2JOULE * NM2A; + real conversion = gmx::c_cal2Joule * gmx::c_nm2A; for (i = 0; i < nforces; i++) diff --git a/src/gromacs/listed_forces/bonded.cpp b/src/gromacs/listed_forces/bonded.cpp index 80ae6741de..5460db9b6a 100644 --- a/src/gromacs/listed_forces/bonded.cpp +++ b/src/gromacs/listed_forces/bonded.cpp @@ -734,7 +734,7 @@ real polarize(int nbonds, type = forceatoms[i++]; ai = forceatoms[i++]; aj = forceatoms[i++]; - ksh = gmx::square(md->chargeA[aj]) * ONE_4PI_EPS0 / forceparams[type].polarize.alpha; + ksh = gmx::square(md->chargeA[aj]) * gmx::c_one4PiEps0 / forceparams[type].polarize.alpha; ki = pbc_rvec_sub(pbc, x[ai], x[aj], dx); /* 3 */ dr2 = iprod(dx, dx); /* 5 */ @@ -781,7 +781,7 @@ real anharm_polarize(int nbonds, type = forceatoms[i++]; ai = forceatoms[i++]; aj = forceatoms[i++]; - ksh = gmx::square(md->chargeA[aj]) * ONE_4PI_EPS0 / forceparams[type].anharm_polarize.alpha; /* 7*/ + ksh = gmx::square(md->chargeA[aj]) * gmx::c_one4PiEps0 / forceparams[type].anharm_polarize.alpha; /* 7*/ khyp = forceparams[type].anharm_polarize.khyp; drcut = forceparams[type].anharm_polarize.drcut; @@ -841,9 +841,9 @@ real water_pol(int nbonds, type0 = forceatoms[0]; aS = forceatoms[5]; qS = md->chargeA[aS]; - kk[XX] = gmx::square(qS) * ONE_4PI_EPS0 / forceparams[type0].wpol.al_x; - kk[YY] = gmx::square(qS) * ONE_4PI_EPS0 / forceparams[type0].wpol.al_y; - kk[ZZ] = gmx::square(qS) * ONE_4PI_EPS0 / forceparams[type0].wpol.al_z; + kk[XX] = gmx::square(qS) * gmx::c_one4PiEps0 / forceparams[type0].wpol.al_x; + kk[YY] = gmx::square(qS) * gmx::c_one4PiEps0 / forceparams[type0].wpol.al_y; + kk[ZZ] = gmx::square(qS) * gmx::c_one4PiEps0 / forceparams[type0].wpol.al_z; r_HH = 1.0 / forceparams[type0].wpol.rHH; for (i = 0; (i < nbonds); i += 6) { @@ -935,7 +935,7 @@ do_1_thole(const rvec xi, const rvec xj, rvec fi, rvec fj, const t_pbc* pbc, rea r12sq = iprod(r12, r12); /* 5 */ r12_1 = gmx::invsqrt(r12sq); /* 5 */ r12bar = afac / r12_1; /* 5 */ - v0 = qq * ONE_4PI_EPS0 * r12_1; /* 2 */ + v0 = qq * gmx::c_one4PiEps0 * r12_1; /* 2 */ ebar = std::exp(-r12bar); /* 5 */ v1 = (1 - (1 + 0.5 * r12bar) * ebar); /* 4 */ fscal = ((v0 * r12_1) * v1 - v0 * 0.5 * afac * ebar * (r12bar + 1)) * r12_1; /* 9 */ @@ -1041,8 +1041,8 @@ angles(int nbonds, *dvdlambda += harmonic(forceparams[type].harmonic.krA, forceparams[type].harmonic.krB, - forceparams[type].harmonic.rA * DEG2RAD, - forceparams[type].harmonic.rB * DEG2RAD, + forceparams[type].harmonic.rA * gmx::c_deg2Rad, + forceparams[type].harmonic.rB * gmx::c_deg2Rad, theta, lambda, &va, @@ -1126,7 +1126,7 @@ angles(int nbonds, alignas(GMX_SIMD_ALIGNMENT) std::int32_t aj[GMX_SIMD_REAL_WIDTH]; alignas(GMX_SIMD_ALIGNMENT) std::int32_t ak[GMX_SIMD_REAL_WIDTH]; alignas(GMX_SIMD_ALIGNMENT) real coeff[2 * GMX_SIMD_REAL_WIDTH]; - SimdReal deg2rad_S(DEG2RAD); + SimdReal deg2rad_S(gmx::c_deg2Rad); SimdReal xi_S, yi_S, zi_S; SimdReal xj_S, yj_S, zj_S; SimdReal xk_S, yk_S, zk_S; @@ -1364,11 +1364,11 @@ urey_bradley(int nbonds, ai = forceatoms[i++]; aj = forceatoms[i++]; ak = forceatoms[i++]; - th0A = forceparams[type].u_b.thetaA * DEG2RAD; + th0A = forceparams[type].u_b.thetaA * gmx::c_deg2Rad; kthA = forceparams[type].u_b.kthetaA; r13A = forceparams[type].u_b.r13A; kUBA = forceparams[type].u_b.kUBA; - th0B = forceparams[type].u_b.thetaB * DEG2RAD; + th0B = forceparams[type].u_b.thetaB * gmx::c_deg2Rad; kthB = forceparams[type].u_b.kthetaB; r13B = forceparams[type].u_b.r13B; kUBB = forceparams[type].u_b.kUBB; @@ -1527,7 +1527,7 @@ urey_bradley(int nbonds, SimdReal rikz_S = zi_S - zk_S; const SimdReal ktheta_S = load(coeff); - const SimdReal theta0_S = load(coeff + GMX_SIMD_REAL_WIDTH) * DEG2RAD; + const SimdReal theta0_S = load(coeff + GMX_SIMD_REAL_WIDTH) * gmx::c_deg2Rad; const SimdReal kUB_S = load(coeff + 2 * GMX_SIMD_REAL_WIDTH); const SimdReal r13_S = load(coeff + 3 * GMX_SIMD_REAL_WIDTH); @@ -1621,7 +1621,7 @@ real quartic_angles(int nbonds, theta = bond_angle(x[ai], x[aj], x[ak], pbc, r_ij, r_kj, &cos_theta, &t1, &t2); /* 41 */ - dt = theta - forceparams[type].qangle.theta * DEG2RAD; /* 2 */ + dt = theta - forceparams[type].qangle.theta * gmx::c_deg2Rad; /* 2 */ dVdt = 0; va = forceparams[type].qangle.c[0]; @@ -1913,8 +1913,8 @@ template real dopdihs(real cpA, real cpB, real phiA, real phiB, int mult, real phi, real lambda, real* V, real* dvdlambda) { const real L1 = 1.0 - lambda; - const real ph0 = (L1 * phiA + lambda * phiB) * DEG2RAD; - const real dph0 = (phiB - phiA) * DEG2RAD; + const real ph0 = (L1 * phiA + lambda * phiB) * gmx::c_deg2Rad; + const real dph0 = (phiB - phiA) * gmx::c_deg2Rad; const real cp = L1 * cpA + lambda * cpB; const real mdphi = mult * phi - ph0; @@ -1937,8 +1937,8 @@ real dopdihs_min(real cpA, real cpB, real phiA, real phiB, int mult, real phi, r { real v, dvdlambda, mdphi, v1, sdphi, ddphi; real L1 = 1.0 - lambda; - real ph0 = (L1 * phiA + lambda * phiB) * DEG2RAD; - real dph0 = (phiB - phiA) * DEG2RAD; + real ph0 = (L1 * phiA + lambda * phiB) * gmx::c_deg2Rad; + real dph0 = (phiB - phiA) * gmx::c_deg2Rad; real cp = L1 * cpA + lambda * cpB; mdphi = mult * (phi - ph0); @@ -2046,7 +2046,7 @@ pdihs(int nbonds, alignas(GMX_SIMD_ALIGNMENT) std::int32_t al[GMX_SIMD_REAL_WIDTH]; alignas(GMX_SIMD_ALIGNMENT) real buf[3 * GMX_SIMD_REAL_WIDTH]; real * cp, *phi0, *mult; - SimdReal deg2rad_S(DEG2RAD); + SimdReal deg2rad_S(gmx::c_deg2Rad); SimdReal p_S, q_S; SimdReal phi0_S, phi_S; SimdReal mx_S, my_S, mz_S; @@ -2313,8 +2313,8 @@ real idihs(int nbonds, pB = forceparams[type].harmonic.rB; kk = L1 * kA + lambda * kB; - phi0 = (L1 * pA + lambda * pB) * DEG2RAD; - dphi0 = (pB - pA) * DEG2RAD; + phi0 = (L1 * pA + lambda * pB) * gmx::c_deg2Rad; + dphi0 = (pB - pA) * gmx::c_deg2Rad; dp = phi - phi0; @@ -2495,7 +2495,7 @@ real dihres(int nbonds, L1 = 1.0 - lambda; - d2r = DEG2RAD; + d2r = gmx::c_deg2Rad; for (i = 0; (i < nbonds);) { @@ -3323,8 +3323,8 @@ real cmap_dihs(int nbonds, /* Switch to degrees */ dx = 360.0 / cmap_grid->grid_spacing; - xphi1 = xphi1 * RAD2DEG; - xphi2 = xphi2 * RAD2DEG; + xphi1 = xphi1 * gmx::c_rad2Deg; + xphi2 = xphi2 * gmx::c_rad2Deg; for (i = 0; i < 4; i++) /* 16 */ { @@ -3361,7 +3361,7 @@ real cmap_dihs(int nbonds, df2 = tt * df2 + (3.0 * tc[i * 4 + 3] * tu + 2.0 * tc[i * 4 + 2]) * tu + tc[i * 4 + 1]; } - fac = RAD2DEG / dx; + fac = gmx::c_rad2Deg / dx; df1 = df1 * fac; df2 = df2 * fac; diff --git a/src/gromacs/listed_forces/restcbt.cpp b/src/gromacs/listed_forces/restcbt.cpp index 2da9dcf8fe..68326ebca3 100644 --- a/src/gromacs/listed_forces/restcbt.cpp +++ b/src/gromacs/listed_forces/restcbt.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2014,2015,2017,2019, by the GROMACS development team, led by + * Copyright (c) 2014,2015,2017,2019,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -78,7 +78,7 @@ void compute_factors_restangles(int type, double term_theta_theta_equil; k_bending = forceparams[type].harmonic.krA; - theta_equil = forceparams[type].harmonic.rA * DEG2RAD; + theta_equil = forceparams[type].harmonic.rA * gmx::c_deg2Rad; theta_equil = M_PI - theta_equil; cosine_theta_equil = cos(theta_equil); @@ -135,7 +135,7 @@ void compute_factors_restrdihs(int type, real norm_phi; /* Read parameters phi0 and k_torsion */ - phi0 = forceparams[type].pdihs.phiA * DEG2RAD; + phi0 = forceparams[type].pdihs.phiA * gmx::c_deg2Rad; cosine_phi0 = cos(phi0); k_torsion = forceparams[type].pdihs.cpA; diff --git a/src/gromacs/mdlib/calc_verletbuf.cpp b/src/gromacs/mdlib/calc_verletbuf.cpp index e7cd5759d3..f29ff7303b 100644 --- a/src/gromacs/mdlib/calc_verletbuf.cpp +++ b/src/gromacs/mdlib/calc_verletbuf.cpp @@ -792,7 +792,7 @@ static real displacementVariance(const t_inputrec& ir, real temperature, real ti * should be negligible (unless nstlist is extremely large, which * you wouldn't do anyhow). */ - kT_fac = 2 * BOLTZ * temperature * timePeriod; + kT_fac = 2 * gmx::c_boltz * temperature * timePeriod; if (ir.bd_fric > 0) { /* This is directly sigma^2 of the displacement */ @@ -813,7 +813,7 @@ static real displacementVariance(const t_inputrec& ir, real temperature, real ti } else { - kT_fac = BOLTZ * temperature * gmx::square(timePeriod); + kT_fac = gmx::c_boltz * temperature * gmx::square(timePeriod); } return kT_fac; @@ -981,7 +981,7 @@ real calcVerletBufferSize(const gmx_mtop_t& mtop, "interactions"); } - elfac = ONE_4PI_EPS0 / ir.epsilon_r; + elfac = gmx::c_one4PiEps0 / ir.epsilon_r; // Determine the 1st and 2nd derivative for the electostatics pot_derivatives_t elec = { 0, 0, 0 }; diff --git a/src/gromacs/mdlib/calcmu.cpp b/src/gromacs/mdlib/calcmu.cpp index 83a61bed78..982faf33e0 100644 --- a/src/gromacs/mdlib/calcmu.cpp +++ b/src/gromacs/mdlib/calcmu.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2012,2014,2015,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -79,7 +79,7 @@ void calc_mu(int start, for (m = 0; (m < DIM); m++) { - mu[m] *= ENM2DEBYE; + mu[m] *= gmx::c_enm2Debye; } if (nChargePerturbed) @@ -94,9 +94,9 @@ void calc_mu(int start, mu_y += qB[i] * x[i][YY]; mu_z += qB[i] * x[i][ZZ]; } - mu_B[XX] = mu_x * ENM2DEBYE; - mu_B[YY] = mu_y * ENM2DEBYE; - mu_B[ZZ] = mu_z * ENM2DEBYE; + mu_B[XX] = mu_x * gmx::c_enm2Debye; + mu_B[YY] = mu_y * gmx::c_enm2Debye; + mu_B[ZZ] = mu_z * gmx::c_enm2Debye; } else { diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp index 9c3f8ef0cc..005376d122 100644 --- a/src/gromacs/mdlib/coupling.cpp +++ b/src/gromacs/mdlib/coupling.cpp @@ -446,7 +446,7 @@ static void NHC_trotter(const t_grpopts* opts, Ekin = 2 * trace(tcstat->ekinh) * tcstat->ekinscaleh_nhc; } } - kT = BOLTZ * reft; + kT = gmx::c_boltz * reft; for (mi = 0; mi < mstepsi; mi++) { @@ -572,7 +572,8 @@ static void boxv_trotter(const t_inputrec* ir, pscal = calc_pres(ir->pbcType, nwall, box, ekinmod, vir, localpres) + pcorr; vol = det(box); - GW = (vol * (MassQ->Winv / PRESFAC)) * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */ + GW = (vol * (MassQ->Winv / gmx::c_presfac)) + * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */ *veta += 0.5 * dt * GW; } @@ -601,7 +602,7 @@ real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, * het systeem... */ - fac = PRESFAC * 2.0 / det(box); + fac = gmx::c_presfac * 2.0 / det(box); for (n = 0; (n < DIM); n++) { for (m = 0; (m < DIM); m++) @@ -625,7 +626,7 @@ real calc_temp(real ekin, real nrdf) { if (nrdf > 0) { - return (2.0 * ekin) / (nrdf * BOLTZ); + return (2.0 * ekin) / (nrdf * gmx::c_boltz); } else { @@ -633,7 +634,7 @@ real calc_temp(real ekin, real nrdf) } } -/*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */ +/*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: c_presfac is not included, so not in GROMACS units! */ static void calcParrinelloRahmanInvMass(const t_inputrec* ir, const matrix box, tensor wInv) { real maxBoxLength; @@ -696,7 +697,7 @@ void parrinellorahman_pcoupl(FILE* fplog, if (!bFirstStep) { - /* Note that PRESFAC does not occur here. + /* Note that c_presfac does not occur here. * The pressure and compressibility always occur as a product, * therefore the pressure unit drops out. */ @@ -972,7 +973,7 @@ void calculateScalingMatrixImplDetail(const t_inputr } real gauss = 0; real gauss2 = 0; - real kt = ir->opts.ref_t[0] * BOLTZ; + real kt = ir->opts.ref_t[0] * gmx::c_boltz; if (kt < 0.0) { kt = 0.0; @@ -986,7 +987,7 @@ void calculateScalingMatrixImplDetail(const t_inputr { const real factor = compressibilityFactor(d, d, ir, dt); mu[d][d] = std::exp(-factor * (ir->ref_p[d][d] - scalar_pressure) / DIM - + std::sqrt(2.0 * kt * factor * PRESFAC / vol) * gauss / DIM); + + std::sqrt(2.0 * kt * factor * gmx::c_presfac / vol) * gauss / DIM); } break; case PressureCouplingType::SemiIsotropic: @@ -996,13 +997,13 @@ void calculateScalingMatrixImplDetail(const t_inputr { const real factor = compressibilityFactor(d, d, ir, dt); mu[d][d] = std::exp(-factor * (ir->ref_p[d][d] - xy_pressure) / DIM - + std::sqrt((DIM - 1) * 2.0 * kt * factor * PRESFAC / vol / DIM) + + std::sqrt((DIM - 1) * 2.0 * kt * factor * gmx::c_presfac / vol / DIM) / (DIM - 1) * gauss); } { const real factor = compressibilityFactor(ZZ, ZZ, ir, dt); mu[ZZ][ZZ] = std::exp(-factor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM - + std::sqrt(2.0 * kt * factor * PRESFAC / vol / DIM) * gauss2); + + std::sqrt(2.0 * kt * factor * gmx::c_presfac / vol / DIM) * gauss2); } break; case PressureCouplingType::SurfaceTension: @@ -1014,12 +1015,12 @@ void calculateScalingMatrixImplDetail(const t_inputr /* Notice: we here use ref_p[ZZ][ZZ] as isotropic pressure and ir->ref_p[d][d] as surface tension */ mu[d][d] = std::exp( -factor * (ir->ref_p[ZZ][ZZ] - ir->ref_p[d][d] / box[ZZ][ZZ] - xy_pressure) / DIM - + std::sqrt(4.0 / 3.0 * kt * factor * PRESFAC / vol) / (DIM - 1) * gauss); + + std::sqrt(4.0 / 3.0 * kt * factor * gmx::c_presfac / vol) / (DIM - 1) * gauss); } { const real factor = compressibilityFactor(ZZ, ZZ, ir, dt); mu[ZZ][ZZ] = std::exp(-factor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM - + std::sqrt(2.0 / 3.0 * kt * factor * PRESFAC / vol) * gauss2); + + std::sqrt(2.0 / 3.0 * kt * factor * gmx::c_presfac / vol) * gauss2); } break; default: @@ -1491,16 +1492,16 @@ extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* Mas /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */ /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */ - MassQ->Winv = (PRESFAC * trace(ir->compress) * BOLTZ * opts->ref_t[0]) + MassQ->Winv = (gmx::c_presfac * trace(ir->compress) * gmx::c_boltz * opts->ref_t[0]) / (DIM * state->vol0 * gmx::square(ir->tau_p / M_2PI)); /* An alternate mass definition, from Tuckerman et al. */ - /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */ + /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*c_boltz*opts->ref_t[0]); */ for (d = 0; d < DIM; d++) { for (n = 0; n < DIM; n++) { - MassQ->Winvm[d][n] = - PRESFAC * ir->compress[d][n] / (state->vol0 * gmx::square(ir->tau_p / M_2PI)); + MassQ->Winvm[d][n] = gmx::c_presfac * ir->compress[d][n] + / (state->vol0 * gmx::square(ir->tau_p / M_2PI)); /* not clear this is correct yet for the anisotropic case. Will need to reevaluate before using MTTK for anisotropic states.*/ } @@ -1518,7 +1519,7 @@ extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* Mas { reft = std::max(0, opts->ref_t[i]); nd = opts->nrdf[i]; - kT = BOLTZ * reft; + kT = gmx::c_boltz * reft; for (j = 0; j < nh; j++) { if (j == 0) @@ -1695,7 +1696,7 @@ init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool b if ((ir->tau_p > 0) && (opts->ref_t[0] > 0)) { reft = std::max(0, opts->ref_t[0]); - kT = BOLTZ * reft; + kT = gmx::c_boltz * reft; for (i = 0; i < nnhpres; i++) { for (j = 0; j < nh; j++) @@ -1740,7 +1741,7 @@ static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t real nd = ir->opts.nrdf[i]; real reft = std::max(ir->opts.ref_t[i], 0); - real kT = BOLTZ * reft; + real kT = gmx::c_boltz * reft; if (nd > 0.0) { @@ -1768,7 +1769,7 @@ static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t } else /* Other non Trotter temperature NH control -- no chains yet. */ { - energy += 0.5 * BOLTZ * nd * gmx::square(ivxi[0]) / iQinv[0]; + energy += 0.5 * gmx::c_boltz * nd * gmx::square(ivxi[0]) / iQinv[0]; energy += nd * ixi[0] * kT; } } @@ -1788,7 +1789,7 @@ static real energyPressureMTTK(const t_inputrec* ir, const t_state* state, const { /* note -- assumes only one degree of freedom that is thermostatted in barostat */ real reft = std::max(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */ - real kT = BOLTZ * reft; + real kT = gmx::c_boltz * reft; for (int j = 0; j < nh; j++) { @@ -1849,7 +1850,8 @@ real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* Mas { if (invMass[d][n] > 0) { - energyNPT += 0.5 * gmx::square(state->boxv[d][n]) / (invMass[d][n] * PRESFAC); + energyNPT += 0.5 * gmx::square(state->boxv[d][n]) + / (invMass[d][n] * gmx::c_presfac); } } } @@ -1861,7 +1863,7 @@ real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* Mas * track of unwrapped box diagonal elements. This case is * excluded in integratorHasConservedEnergyQuantity(). */ - energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC); + energyNPT += vol * trace(ir->ref_p) / (DIM * gmx::c_presfac); break; } case PressureCoupling::Mttk: @@ -1869,7 +1871,7 @@ real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* Mas energyNPT += 0.5 * gmx::square(state->veta) / MassQ->Winv; /* contribution from the PV term */ - energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC); + energyNPT += vol * trace(ir->ref_p) / (DIM * gmx::c_presfac); if (ir->epc == PressureCoupling::Mttk) { @@ -2010,7 +2012,7 @@ void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0) { - Ek_ref1 = 0.5 * opts->ref_t[i] * BOLTZ; + Ek_ref1 = 0.5 * opts->ref_t[i] * gmx::c_boltz; Ek_ref = Ek_ref1 * opts->nrdf[i]; Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i], opts->tau_t[i] / dt, step, ir->ld_seed); diff --git a/src/gromacs/mdlib/dispersioncorrection.cpp b/src/gromacs/mdlib/dispersioncorrection.cpp index 481c8597e7..2c8976f89a 100644 --- a/src/gromacs/mdlib/dispersioncorrection.cpp +++ b/src/gromacs/mdlib/dispersioncorrection.cpp @@ -634,7 +634,7 @@ DispersionCorrection::Correction DispersionCorrection::calculate(const matrix bo corr.virial += numCorr * density * avctwelve * iParams_.virdifftwelve_ / 3.0; } /* The factor 2 is because of the Gromacs virial definition */ - corr.pressure = -2.0 * invvol * corr.virial * PRESFAC; + corr.pressure = -2.0 * invvol * corr.virial * gmx::c_presfac; } if (eFep_ != FreeEnergyPerturbationType::No) diff --git a/src/gromacs/mdlib/energyoutput.cpp b/src/gromacs/mdlib/energyoutput.cpp index 571edc54c8..1f75ff4781 100644 --- a/src/gromacs/mdlib/energyoutput.cpp +++ b/src/gromacs/mdlib/energyoutput.cpp @@ -908,7 +908,7 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, nboxs = boxs_nm.size(); } vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ]; - dens = (tmass * AMU) / (vol * NANO * NANO * NANO); + dens = (tmass * gmx::c_amu) / (vol * gmx::c_nano * gmx::c_nano * gmx::c_nano); add_ebin(ebin_, ib_, nboxs, bs, bSum); add_ebin(ebin_, ivol_, 1, &vol, bSum); add_ebin(ebin_, idens_, 1, &dens, bSum); @@ -917,7 +917,7 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, { /* This is pV (in kJ/mol). The pressure is the reference pressure, not the instantaneous pressure */ - pv = vol * ref_p_ / PRESFAC; + pv = vol * ref_p_ / gmx::c_presfac; add_ebin(ebin_, ipv_, 1, &pv, bSum); enthalpy = pv + enerd->term[F_ETOT]; @@ -953,12 +953,12 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, if (ekind && ekind->cosacc.cos_accel != 0) { vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ]; - dens = (tmass * AMU) / (vol * NANO * NANO * NANO); + dens = (tmass * gmx::c_amu) / (vol * gmx::c_nano * gmx::c_nano * gmx::c_nano); add_ebin(ebin_, ivcos_, 1, &(ekind->cosacc.vcos), bSum); /* 1/viscosity, unit 1/(kg m^-1 s^-1) */ tmp = 1 - / (ekind->cosacc.cos_accel / (ekind->cosacc.vcos * PICO) * dens - * gmx::square(box[ZZ][ZZ] * NANO / (2 * M_PI))); + / (ekind->cosacc.cos_accel / (ekind->cosacc.vcos * gmx::c_pico) * dens + * gmx::square(box[ZZ][ZZ] * gmx::c_nano / (2 * M_PI))); add_ebin(ebin_, ivisc_, 1, &tmp, bSum); } if (nE_ > 1) diff --git a/src/gromacs/mdlib/expanded.cpp b/src/gromacs/mdlib/expanded.cpp index 965d92b441..2a6837f6bd 100644 --- a/src/gromacs/mdlib/expanded.cpp +++ b/src/gromacs/mdlib/expanded.cpp @@ -52,6 +52,7 @@ #include "gromacs/listed_forces/orires.h" #include "gromacs/math/functions.h" #include "gromacs/math/units.h" +#include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/calcmu.h" #include "gromacs/mdlib/constr.h" @@ -1402,15 +1403,16 @@ int ExpandedEnsembleDynamics(FILE* log, if (ir->bSimTemp) { /* Note -- this assumes no mass changes, since kinetic energy is not added . . . */ - scaled_lamee[i] = enerd->foreignLambdaTerms.deltaH(i) / (simtemp->temperatures[i] * BOLTZ) - + enerd->term[F_EPOT] - * (1.0 / (simtemp->temperatures[i]) - - 1.0 / (simtemp->temperatures[fep_state])) - / BOLTZ; + scaled_lamee[i] = + enerd->foreignLambdaTerms.deltaH(i) / (simtemp->temperatures[i] * gmx::c_boltz) + + enerd->term[F_EPOT] + * (1.0 / (simtemp->temperatures[i]) + - 1.0 / (simtemp->temperatures[fep_state])) + / gmx::c_boltz; } else { - scaled_lamee[i] = enerd->foreignLambdaTerms.deltaH(i) / (expand->mc_temp * BOLTZ); + scaled_lamee[i] = enerd->foreignLambdaTerms.deltaH(i) / (expand->mc_temp * gmx::c_boltz); /* mc_temp is currently set to the system reft unless otherwise defined */ } @@ -1427,7 +1429,8 @@ int ExpandedEnsembleDynamics(FILE* log, { scaled_lamee[i] = enerd->term[F_EPOT] - * (1.0 / simtemp->temperatures[i] - 1.0 / simtemp->temperatures[fep_state]) / BOLTZ; + * (1.0 / simtemp->temperatures[i] - 1.0 / simtemp->temperatures[fep_state]) + / gmx::c_boltz; } } } diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index e231cb4b12..cda7d5743e 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -893,7 +893,7 @@ static interaction_const_t init_interaction_const(FILE* fp, /* Set the Coulomb energy conversion factor */ if (interactionConst.epsilon_r != 0) { - interactionConst.epsfac = ONE_4PI_EPS0 / interactionConst.epsilon_r; + interactionConst.epsfac = gmx::c_one4PiEps0 / interactionConst.epsilon_r; } else { diff --git a/src/gromacs/mdlib/lincs.cpp b/src/gromacs/mdlib/lincs.cpp index 5b24e8be3f..b5fe54be92 100644 --- a/src/gromacs/mdlib/lincs.cpp +++ b/src/gromacs/mdlib/lincs.cpp @@ -1096,7 +1096,7 @@ static void do_lincs(ArrayRefWithPadding xPadded, real wfac; - wfac = std::cos(DEG2RAD * wangle); + wfac = std::cos(gmx::c_deg2Rad * wangle); wfac = wfac * wfac; for (int iter = 0; iter < lincsd->nIter; iter++) @@ -2151,7 +2151,7 @@ static void lincs_warning(gmx_domdec_t* dd, int maxwarn, int* warncount) { - real wfac = std::cos(DEG2RAD * wangle); + real wfac = std::cos(gmx::c_deg2Rad * wangle); fprintf(stderr, "bonds that rotated more than %g degrees:\n" @@ -2183,7 +2183,7 @@ static void lincs_warning(gmx_domdec_t* dd, " %6d %6d %5.1f %8.4f %8.4f %8.4f\n", ddglatnr(dd, i), ddglatnr(dd, j), - RAD2DEG * std::acos(cosine), + gmx::c_rad2Deg * std::acos(cosine), d0, d1, bllen[b]); diff --git a/src/gromacs/mdlib/rf_util.cpp b/src/gromacs/mdlib/rf_util.cpp index 4584876236..3eea60f678 100644 --- a/src/gromacs/mdlib/rf_util.cpp +++ b/src/gromacs/mdlib/rf_util.cpp @@ -74,7 +74,7 @@ void calc_rffac(FILE* fplog, real eps_r, real eps_rf, real Rc, real* krf, real* Rc, *krf, *crf, - ONE_4PI_EPS0 / eps_r); + gmx::c_one4PiEps0 / eps_r); // Make sure we don't lose resolution in pow() by casting real arg to double real rmin = gmx::invcbrt(static_cast(*krf * 2.0)); fprintf(fplog, "The electrostatics potential has its minimum at r = %g\n", rmin); diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index b9e9bbc4a2..0bcf0c5451 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -485,7 +485,7 @@ static real averageKineticEnergyEstimate(const t_grpopts& groupOptions) if (groupOptions.tau_t[g] >= 0) { nrdfCoupled += groupOptions.nrdf[g]; - kineticEnergy += groupOptions.nrdf[g] * 0.5 * groupOptions.ref_t[g] * BOLTZ; + kineticEnergy += groupOptions.nrdf[g] * 0.5 * groupOptions.ref_t[g] * gmx::c_boltz; } else { diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp index 6b695e5c52..16c722e915 100644 --- a/src/gromacs/mdlib/update.cpp +++ b/src/gromacs/mdlib/update.cpp @@ -891,7 +891,7 @@ gmx_stochd_t::gmx_stochd_t(const t_inputrec& inputRecord) && (reft > 0)) /* tau_t or ref_t = 0 means that no randomization is done */ { randomize_group[gt] = true; - boltzfac[gt] = BOLTZ * opts->ref_t[gt]; + boltzfac[gt] = gmx::c_boltz * opts->ref_t[gt]; } else { @@ -909,7 +909,7 @@ void Update::Impl::update_temperature_constants(const t_inputrec& inputRecord) { for (int gt = 0; gt < inputRecord.opts.ngtc; gt++) { - sd_.bd_rf[gt] = std::sqrt(2.0 * BOLTZ * inputRecord.opts.ref_t[gt] + sd_.bd_rf[gt] = std::sqrt(2.0 * gmx::c_boltz * inputRecord.opts.ref_t[gt] / (inputRecord.bd_fric * inputRecord.delta_t)); } } @@ -917,7 +917,7 @@ void Update::Impl::update_temperature_constants(const t_inputrec& inputRecord) { for (int gt = 0; gt < inputRecord.opts.ngtc; gt++) { - sd_.bd_rf[gt] = std::sqrt(2.0 * BOLTZ * inputRecord.opts.ref_t[gt]); + sd_.bd_rf[gt] = std::sqrt(2.0 * gmx::c_boltz * inputRecord.opts.ref_t[gt]); } } } @@ -925,7 +925,7 @@ void Update::Impl::update_temperature_constants(const t_inputrec& inputRecord) { for (int gt = 0; gt < inputRecord.opts.ngtc; gt++) { - real kT = BOLTZ * inputRecord.opts.ref_t[gt]; + real kT = gmx::c_boltz * inputRecord.opts.ref_t[gt]; /* The mass is accounted for later, since this differs per atom */ sd_.sdsig[gt].V = std::sqrt(kT * (1 - sd_.sdc[gt].em * sd_.sdc[gt].em)); } diff --git a/src/gromacs/mdlib/updategroups.cpp b/src/gromacs/mdlib/updategroups.cpp index b8ed8dcadb..3896550d99 100644 --- a/src/gromacs/mdlib/updategroups.cpp +++ b/src/gromacs/mdlib/updategroups.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -522,9 +522,10 @@ static real constraintGroupRadius(const gmx_moltype_t& molty /* Set number of stddevs such that change of exceeding < 10^-9 */ constexpr real c_numSigma = 6.0; /* Compute the maximally stretched angle */ - const real eqAngle = angleParams.harmonic.rA * DEG2RAD; + const real eqAngle = angleParams.harmonic.rA * gmx::c_deg2Rad; const real fc = angleParams.harmonic.krA; - const real maxAngle = eqAngle + c_numSigma * BOLTZ * temperature / ((numPartnerAtoms - 1) * fc); + const real maxAngle = + eqAngle + c_numSigma * gmx::c_boltz * temperature / ((numPartnerAtoms - 1) * fc); if (maxAngle >= M_PI) { return -1; diff --git a/src/gromacs/mdrun/replicaexchange.cpp b/src/gromacs/mdrun/replicaexchange.cpp index d813cfdc97..67e19f212e 100644 --- a/src/gromacs/mdrun/replicaexchange.cpp +++ b/src/gromacs/mdrun/replicaexchange.cpp @@ -845,7 +845,7 @@ static real calc_delta(FILE* fplog, gmx_bool bPrint, struct gmx_repl_ex* re, int if (re->bNPT) { /* revist the calculation for 5.0. Might be some improvements. */ - dpV = (beta[ap] * re->pres[ap] - beta[bp] * re->pres[bp]) * (Vol[b] - Vol[a]) / PRESFAC; + dpV = (beta[ap] * re->pres[ap] - beta[bp] * re->pres[bp]) * (Vol[b] - Vol[a]) / gmx::c_presfac; if (bPrint) { fprintf(fplog, " dpV = %10.3e d = %10.3e\n", dpV, delta + dpV); @@ -899,14 +899,14 @@ static void test_for_replica_exchange(FILE* fplog, /* temperatures of different states*/ for (i = 0; i < re->nrepl; i++) { - re->beta[i] = 1.0 / (re->q[ereTEMP][i] * BOLTZ); + re->beta[i] = 1.0 / (re->q[ereTEMP][i] * gmx::c_boltz); } } else { for (i = 0; i < re->nrepl; i++) { - re->beta[i] = 1.0 / (re->temp * BOLTZ); /* we have a single temperature */ + re->beta[i] = 1.0 / (re->temp * gmx::c_boltz); /* we have a single temperature */ } } if (re->type == ereLAMBDA || re->type == ereTL) diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index cd008e3b9d..dfd8df643c 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -54,6 +54,7 @@ #include "gromacs/gmxlib/network.h" #include "gromacs/math/functions.h" #include "gromacs/math/units.h" +#include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/math/vecdump.h" #include "gromacs/mdlib/constr.h" @@ -454,7 +455,7 @@ gmx_shellfc_t* init_shell_flexcon(FILE* fplog, aS + 1, mb + 1); } - shell[nsi].k += gmx::square(qS) * ONE_4PI_EPS0 + shell[nsi].k += gmx::square(qS) * gmx::c_one4PiEps0 / ffparams->iparams[type].polarize.alpha; break; case F_WATER_POL: @@ -472,7 +473,7 @@ gmx_shellfc_t* init_shell_flexcon(FILE* fplog, + ffparams->iparams[type].wpol.al_y + ffparams->iparams[type].wpol.al_z) / 3.0; - shell[nsi].k += gmx::square(qS) * ONE_4PI_EPS0 / alpha; + shell[nsi].k += gmx::square(qS) * gmx::c_one4PiEps0 / alpha; break; default: gmx_fatal(FARGS, "Death Horror: %s, %d", __FILE__, __LINE__); } diff --git a/src/gromacs/mdrun/tpi.cpp b/src/gromacs/mdrun/tpi.cpp index 2eda488785..debaef5c31 100644 --- a/src/gromacs/mdrun/tpi.cpp +++ b/src/gromacs/mdrun/tpi.cpp @@ -277,7 +277,7 @@ void LegacySimulator::do_tpi() } fprintf(fplog, "\n The temperature for test particle insertion is %.3f K\n\n", temp); } - beta = 1.0 / (BOLTZ * temp); + beta = 1.0 / (gmx::c_boltz * temp); /* Number of insertions per frame */ nsteps = inputrec->nsteps; diff --git a/src/gromacs/mimic/communicator.cpp b/src/gromacs/mimic/communicator.cpp index 12dcf874d3..177dc8813b 100644 --- a/src/gromacs/mimic/communicator.cpp +++ b/src/gromacs/mimic/communicator.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -141,7 +141,7 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector(mtop->ffparams.iparams[contype].constr.dA) - / BOHR2NM); + / gmx::c_bohr2Nm); } for (int ncon = 0; ncon < nsettle; ++ncon) @@ -165,11 +165,11 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector(mtop->ffparams.iparams[contype].constr.dA) - / BOHR2NM); + / gmx::c_bohr2Nm); bondLengths.push_back(static_cast(mtop->ffparams.iparams[contype].constr.dA) - / BOHR2NM); + / gmx::c_bohr2Nm); bondLengths.push_back(static_cast(mtop->ffparams.iparams[contype].constr.dB) - / BOHR2NM); + / gmx::c_bohr2Nm); } nAtomsMol.push_back(type->atoms.nr); @@ -234,9 +234,9 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector convertedCoords; for (auto& coord : coords) { - convertedCoords.push_back(static_cast(coord[0]) / BOHR2NM); - convertedCoords.push_back(static_cast(coord[1]) / BOHR2NM); - convertedCoords.push_back(static_cast(coord[2]) / BOHR2NM); + convertedCoords.push_back(static_cast(coord[0]) / gmx::c_bohr2Nm); + convertedCoords.push_back(static_cast(coord[1]) / gmx::c_bohr2Nm); + convertedCoords.push_back(static_cast(coord[2]) / gmx::c_bohr2Nm); } // sending array of coordinates @@ -256,15 +256,15 @@ void gmx::MimicCommunicator::getCoords(PaddedHostVector* x, const int nato MCL_receive(&*coords.begin(), 3 * natoms, TYPE_DOUBLE, 0); for (int j = 0; j < natoms; ++j) { - (*x)[j][0] = static_cast(coords[j * 3] * BOHR2NM); - (*x)[j][1] = static_cast(coords[j * 3 + 1] * BOHR2NM); - (*x)[j][2] = static_cast(coords[j * 3 + 2] * BOHR2NM); + (*x)[j][0] = static_cast(coords[j * 3] * gmx::c_bohr2Nm); + (*x)[j][1] = static_cast(coords[j * 3 + 1] * gmx::c_bohr2Nm); + (*x)[j][2] = static_cast(coords[j * 3 + 2] * gmx::c_bohr2Nm); } } void gmx::MimicCommunicator::sendEnergies(real energy) { - double convertedEnergy = energy / (HARTREE2KJ * AVOGADRO); + double convertedEnergy = energy / (gmx::c_hartree2Kj * gmx::c_avogadro); MCL_send(&convertedEnergy, 1, TYPE_DOUBLE, 0); } @@ -273,9 +273,9 @@ void gmx::MimicCommunicator::sendForces(gmx::ArrayRef forces, int nat std::vector convertedForce; for (int j = 0; j < natoms; ++j) { - convertedForce.push_back(static_cast(forces[j][0]) / HARTREE_BOHR2MD); - convertedForce.push_back(static_cast(forces[j][1]) / HARTREE_BOHR2MD); - convertedForce.push_back(static_cast(forces[j][2]) / HARTREE_BOHR2MD); + convertedForce.push_back(static_cast(forces[j][0]) / gmx::c_hartreeBohr2Md); + convertedForce.push_back(static_cast(forces[j][1]) / gmx::c_hartreeBohr2Md); + convertedForce.push_back(static_cast(forces[j][2]) / gmx::c_hartreeBohr2Md); } MCL_send(&*convertedForce.begin(), convertedForce.size(), TYPE_DOUBLE, 0); } diff --git a/src/gromacs/modularsimulator/parrinellorahmanbarostat.cpp b/src/gromacs/modularsimulator/parrinellorahmanbarostat.cpp index 81085d0cda..9b35825a60 100644 --- a/src/gromacs/modularsimulator/parrinellorahmanbarostat.cpp +++ b/src/gromacs/modularsimulator/parrinellorahmanbarostat.cpp @@ -219,7 +219,7 @@ real ParrinelloRahmanBarostat::conservedEnergyContribution() const { for (int j = 0; j <= i; j++) { - real invMass = PRESFAC * (4 * M_PI * M_PI * inputrec_->compress[i][j]) + real invMass = c_presfac * (4 * M_PI * M_PI * inputrec_->compress[i][j]) / (3 * inputrec_->tau_p * inputrec_->tau_p * maxBoxLength); if (invMass > 0) { @@ -235,7 +235,7 @@ real ParrinelloRahmanBarostat::conservedEnergyContribution() const * track of unwrapped box diagonal elements. This case is * excluded in integratorHasConservedEnergyQuantity(). */ - energy += volume * trace(inputrec_->ref_p) / (DIM * PRESFAC); + energy += volume * trace(inputrec_->ref_p) / (DIM * c_presfac); return energy; } diff --git a/src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.cpp b/src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.cpp index 95158d5034..402b13a954 100644 --- a/src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.cpp +++ b/src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.cpp @@ -138,7 +138,7 @@ public: } const real referenceKineticEnergy = - 0.5 * temperatureCouplingData.referenceTemperature[temperatureGroup] * BOLTZ + 0.5 * temperatureCouplingData.referenceTemperature[temperatureGroup] * gmx::c_boltz * temperatureCouplingData.numDegreesOfFreedom[temperatureGroup]; const real newKineticEnergy = diff --git a/src/gromacs/pbcutil/pbc.cpp b/src/gromacs/pbcutil/pbc.cpp index dc66f20558..3e6fbb6623 100644 --- a/src/gromacs/pbcutil/pbc.cpp +++ b/src/gromacs/pbcutil/pbc.cpp @@ -4,7 +4,7 @@ * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -167,7 +167,7 @@ const char* check_box(PbcType pbcType, const matrix box) void matrix_convert(matrix box, const rvec vec, const rvec angleInDegrees) { rvec angle; - svmul(DEG2RAD, angleInDegrees, angle); + svmul(gmx::c_deg2Rad, angleInDegrees, angle); box[XX][XX] = vec[XX]; box[YY][XX] = vec[YY] * cos(angle[ZZ]); box[YY][YY] = vec[YY] * sin(angle[ZZ]); diff --git a/src/gromacs/pulling/pull.cpp b/src/gromacs/pulling/pull.cpp index 41b5bbce40..b0da21552c 100644 --- a/src/gromacs/pulling/pull.cpp +++ b/src/gromacs/pulling/pull.cpp @@ -59,6 +59,7 @@ #include "gromacs/gmxlib/network.h" #include "gromacs/math/functions.h" #include "gromacs/math/units.h" +#include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/math/vectypes.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" @@ -153,7 +154,7 @@ double pull_conversion_factor_userinput2internal(const t_pull_coord* pcrd) { if (pull_coordinate_is_angletype(pcrd)) { - return DEG2RAD; + return gmx::c_deg2Rad; } else { @@ -165,7 +166,7 @@ double pull_conversion_factor_internal2userinput(const t_pull_coord* pcrd) { if (pull_coordinate_is_angletype(pcrd)) { - return RAD2DEG; + return gmx::c_rad2Deg; } else { diff --git a/src/gromacs/selection/params.cpp b/src/gromacs/selection/params.cpp index 7bf73880fb..a5d23a9c48 100644 --- a/src/gromacs/selection/params.cpp +++ b/src/gromacs/selection/params.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 2009,2010,2011,2012,2013 by the GROMACS development team. * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -48,6 +48,7 @@ #include #include "gromacs/math/units.h" +#include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/utility/cstringutil.h" #include "gromacs/utility/exceptions.h" diff --git a/src/gromacs/selection/sm_insolidangle.cpp b/src/gromacs/selection/sm_insolidangle.cpp index bbcffa350e..3dd91a6ab5 100644 --- a/src/gromacs/selection/sm_insolidangle.cpp +++ b/src/gromacs/selection/sm_insolidangle.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 2009,2010,2011,2012,2013 by the GROMACS development team. * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -397,12 +397,12 @@ static void init_insolidangle(const gmx_mtop_t* /* top */, GMX_THROW(gmx::InvalidInputError("Angle cutoff should be > 0")); } - surf->angcut *= DEG2RAD; + surf->angcut *= gmx::c_deg2Rad; surf->distccut = -std::cos(surf->angcut); surf->targetbinsize = surf->angcut / 2; surf->ntbins = static_cast(M_PI / surf->targetbinsize); - surf->tbinsize = (180.0 / surf->ntbins) * DEG2RAD; + surf->tbinsize = (180.0 / surf->ntbins) * gmx::c_deg2Rad; snew(surf->tbin, static_cast(M_PI / surf->tbinsize) + 1); surf->maxbins = 0; diff --git a/src/gromacs/tables/forcetable.cpp b/src/gromacs/tables/forcetable.cpp index 1006e07f1c..1059028d73 100644 --- a/src/gromacs/tables/forcetable.cpp +++ b/src/gromacs/tables/forcetable.cpp @@ -1420,10 +1420,10 @@ bondedtable_t make_bonded_table(FILE* fplog, const char* fn, int angle) /* Convert the table from degrees to radians */ for (i = 0; i < td.nx; i++) { - td.x[i] *= DEG2RAD; - td.f[i] *= RAD2DEG; + td.x[i] *= gmx::c_deg2Rad; + td.f[i] *= gmx::c_rad2Deg; } - td.tabscale *= RAD2DEG; + td.tabscale *= gmx::c_rad2Deg; } tab.n = td.nx; tab.scale = td.tabscale; diff --git a/src/gromacs/tools/check.cpp b/src/gromacs/tools/check.cpp index 856fcee292..d44370c0c6 100644 --- a/src/gromacs/tools/check.cpp +++ b/src/gromacs/tools/check.cpp @@ -494,8 +494,8 @@ static void chk_tps(const char* fn, real vdw_fac, real bon_lo, real bon_hi) ekin += 0.5 * atoms->atom[i].m * v[i][j] * v[i][j]; } } - temp1 = (2.0 * ekin) / (natom * DIM * BOLTZ); - temp2 = (2.0 * ekin) / (natom * (DIM - 1) * BOLTZ); + temp1 = (2.0 * ekin) / (natom * DIM * gmx::c_boltz); + temp2 = (2.0 * ekin) / (natom * (DIM - 1) * gmx::c_boltz); fprintf(stderr, "Kinetic energy: %g (kJ/mol)\n", ekin); fprintf(stderr, "Assuming the number of degrees of freedom to be " diff --git a/src/gromacs/tools/pme_error.cpp b/src/gromacs/tools/pme_error.cpp index 51bebada58..15ce914844 100644 --- a/src/gromacs/tools/pme_error.cpp +++ b/src/gromacs/tools/pme_error.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 2010,2011,2012,2013,2014 by the GROMACS development team. * Copyright (c) 2015,2016,2017,2018,2019 by the GROMACS development team. - * Copyright (c) 2019,2020, by the GROMACS development team, led by + * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -182,7 +182,7 @@ static real estimate_direct(t_inputinfo* info) e_dir = 2.0 * info->q2all * gmx::invsqrt(info->q2allnr * r_coulomb * info->volume); e_dir *= std::exp(-beta * beta * r_coulomb * r_coulomb); - return ONE_4PI_EPS0 * e_dir; + return gmx::c_one4PiEps0 * e_dir; } #define SUMORDER 6 @@ -776,7 +776,7 @@ static real estimate_reciprocal(t_inputinfo* info, e_rec = std::sqrt(e_rec1 + e_rec2 + e_rec3); - return ONE_4PI_EPS0 * e_rec; + return gmx::c_one4PiEps0 * e_rec; } diff --git a/src/gromacs/trajectoryanalysis/modules/angle.cpp b/src/gromacs/trajectoryanalysis/modules/angle.cpp index 80138696cc..53e1010831 100644 --- a/src/gromacs/trajectoryanalysis/modules/angle.cpp +++ b/src/gromacs/trajectoryanalysis/modules/angle.cpp @@ -809,7 +809,7 @@ void Angle::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryA break; default: GMX_THROW(InternalError("invalid -g1 value")); } - dh.setPoint(n, angle * RAD2DEG, bPresent); + dh.setPoint(n, angle * gmx::c_rad2Deg, bPresent); } } dh.finishFrame(); diff --git a/src/gromacs/trajectoryanalysis/modules/freevolume.cpp b/src/gromacs/trajectoryanalysis/modules/freevolume.cpp index f67303c92f..e9f7cbd827 100644 --- a/src/gromacs/trajectoryanalysis/modules/freevolume.cpp +++ b/src/gromacs/trajectoryanalysis/modules/freevolume.cpp @@ -388,7 +388,7 @@ void FreeVolume::writeOutput() printf("Total volume %.2f +/- %.2f nm^3\n", Vaver, Verror); printf("Number of molecules %d total mass %.2f Dalton\n", nmol_, mtot_); - double RhoAver = mtot_ / (Vaver * 1e-24 * AVOGADRO); + double RhoAver = mtot_ / (Vaver * 1e-24 * gmx::c_avogadro); double RhoError = gmx::square(RhoAver / Vaver) * Verror; printf("Average molar mass: %.2f Dalton\n", mtot_ / nmol_); diff --git a/src/gromacs/trajectoryanalysis/modules/sasa.cpp b/src/gromacs/trajectoryanalysis/modules/sasa.cpp index 9ccc4c8f01..00466715da 100644 --- a/src/gromacs/trajectoryanalysis/modules/sasa.cpp +++ b/src/gromacs/trajectoryanalysis/modules/sasa.cpp @@ -1076,7 +1076,7 @@ void Sasa::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAn { totmass += surfaceSel.position(i).mass(); } - const real density = totmass * AMU / (totvolume * NANO * NANO * NANO); + const real density = totmass * gmx::c_amu / (totvolume * gmx::c_nano * gmx::c_nano * gmx::c_nano); vh.startFrame(frnr, fr.time); vh.setPoint(0, totvolume); vh.setPoint(1, density); -- 2.22.0