Thermal-FIST 1.6
Package for hadron resonance gas model applications
Loading...
Searching...
No Matches
ZeroTemperatureQvdW.cpp
Go to the documentation of this file.
1/*
2 * Thermal-FIST package
3 *
4 * Copyright (c) 2026 Volodymyr Vovchenko
5 *
6 * GNU General Public License (GPLv3 or later)
7 */
8#include <cassert>
9#include <cmath>
10#include <cstdio>
11#include <cstdlib>
12#include <fstream>
13#include <iomanip>
14#include <iostream>
15#include <limits>
16#include <memory>
17#include <string>
18#include <vector>
19
20#include "HRGBase.h"
21#include "HRGEV.h"
22#include "HRGRealGas.h"
23#include "HRGVDW.h"
24#include "HRGEV/ExcludedVolumeHelper.h"
25#include "ThermalFISTConfig.h"
26
27using namespace std;
28
29#ifdef ThermalFIST_USENAMESPACE
30using namespace thermalfist;
31#endif
32
33// Computes zero-temperature EoS in the QvdW setup for either:
34// - ThermalModelVDW (useRG=0)
35// - ThermalModelRealGas with CS EV + VDW mean field (useRG=1)
36//
37// Usage:
38// ZeroTemperatureQvdW <useRG> <a> <b> <muBmin> <muBmax> <dmuB> <outputfile>
39//
40// Defaults:
41// useRG = 0 (0: ThermalModelVDW, 1: ThermalModelRealGas)
42// a = 0.329 GeV*fm^3
43// b = 3.42 fm^3
44// muBmin = 0.90 GeV
45// muBmax = 1.70 GeV
46// dmuB = 0.01 GeV
47int main(int argc, char *argv[])
48{
49 bool useRG = false;
50 if (argc > 1) useRG = (atoi(argv[1]) != 0);
51
52 double a = 0.329;
53 if (argc > 2) a = atof(argv[2]);
54
55 double b = 3.42;
56 if (argc > 3) b = atof(argv[3]);
57
58 double muBmin = 0.90;
59 if (argc > 4) muBmin = atof(argv[4]);
60
61 double muBmax = 1.70;
62 if (argc > 5) muBmax = atof(argv[5]);
63
64 double dmuB = 0.01;
65 if (argc > 6) dmuB = atof(argv[6]);
66
67 string outputfile = string("ZeroT-EoS-") + (useRG ? "RG-HRG" : "QvdW-HRG") + ".dat";
68 if (argc > 7) outputfile = argv[7];
69
70 if (dmuB <= 0.0) {
71 cerr << "Error: dmuB must be positive." << endl;
72 return 1;
73 }
74 if (muBmax < muBmin) {
75 cerr << "Error: muBmax must be >= muBmin." << endl;
76 return 1;
77 }
78
79 ThermalParticleSystem tps(string(ThermalFIST_INPUT_FOLDER) + "/list/PDG2025/list.dat");
80
81 unique_ptr<ThermalModelBase> model;
82 if (!useRG) {
83 auto *vdw = new ThermalModelVDW(&tps);
84 SetVDWHRGInteractionParameters(vdw, a, b);
85 model.reset(vdw);
86 cout << "Model: ThermalModelVDW (QvdW-HRG)" << endl;
87 } else {
88 auto *rg = new ThermalModelRealGas(&tps);
89
90 if (b != 0.0) {
91 ExcludedVolumeModelBase *evmod = new ExcludedVolumeModelCS();
92 auto vdWb = GetBaryonBaryonInteractionMatrix(rg->TPS(), b);
93 rg->SetExcludedVolumeModel(new ExcludedVolumeModelCrosstermsGeneralized(evmod, vdWb));
94 }
95 if (a != 0.0) {
96 auto vdWa = GetBaryonBaryonInteractionMatrix(rg->TPS(), a);
97 rg->SetMeanFieldModel(new MeanFieldModelMultiVDW(vdWa));
98 }
99
100 model.reset(rg);
101 cout << "Model: ThermalModelRealGas (CS + VDW mean field)" << endl;
102 }
103
104 model->SetUseWidth(false);
105 model->SetStatistics(true);
106 model->SetCalculationType(IdealGasFunctions::Quadratures);
107
108 // Zero-temperature setup with fixed charge chemical potentials.
109 model->SetTemperature(0.0);
110 model->SetElectricChemicalPotential(0.0);
111 model->SetStrangenessChemicalPotential(0.0);
112 model->SetCharmChemicalPotential(0.0);
113
114 ofstream fout(outputfile.c_str());
115 if (!fout.is_open()) {
116 cerr << "Error: cannot open output file " << outputfile << endl;
117 return 1;
118 }
119
120 const int tab = 20;
121 fout << setw(tab) << "muB[GeV]" << " "
122 << setw(tab) << "P[GeV/fm3]" << " "
123 << setw(tab) << "e[GeV/fm3]" << " "
124 << setw(tab) << "s[fm-3]" << " "
125 << setw(tab) << "rhoB[fm-3]" << " "
126 << setw(tab) << "rhoQ[fm-3]" << " "
127 << setw(tab) << "rhoS[fm-3]" << " "
128 << setw(tab) << "(1/3-p/e)" << " "
129 << setw(tab) << "cT2(B)" << " "
130 << setw(tab) << "cs2(B)" << " "
131 << setw(tab) << "cT2(BQ)" << " "
132 << setw(tab) << "cs2(BQ)" << " "
133 << setw(tab) << "cT2(BS)" << " "
134 << setw(tab) << "cs2(BS)" << " "
135 << setw(tab) << "cT2(BQS)" << " "
136 << setw(tab) << "cs2(BQS)" << " "
137 << setw(tab) << "cV(B)[fm-3]" << " "
138 << setw(tab) << "cV(BQ)[fm-3]" << " "
139 << setw(tab) << "cV(BS)[fm-3]" << " "
140 << setw(tab) << "cV(BQS)[fm-3]" << " "
141 << setw(tab) << "cMu[fm-3]" << " "
142 << setw(tab) << "(ds/dT)_mu_num@1e-4" << " "
143 << setw(tab) << "(ds/dT)_mu_num@1e-3" << " "
144 << setw(tab) << "(ds/dT)_mu_num@1e-2" << " "
145 << setw(tab) << "(ds/dT)_mu_an@1e-4" << " "
146 << setw(tab) << "(ds/dT)_mu_an@1e-3" << " "
147 << setw(tab) << "(ds/dT)_mu_an@1e-2" << " "
148 << setw(tab) << "(ds/dT)_mu_an" << "\n";
149
150 const double wt1 = get_wall_time();
151 int iters = 0;
152
153 for (double muB = muBmin; muB <= muBmax + 0.1 * dmuB; muB += dmuB) {
154 model->SetBaryonChemicalPotential(muB);
155 model->CalculatePrimordialDensities();
156 model->CalculateFluctuations();
157 model->CalculateTemperatureDerivatives();
158
159 const double p = model->Pressure();
160 const double e = model->EnergyDensity();
161 const double s = model->EntropyDensity();
162 const double rhoB = model->BaryonDensity();
163 const double rhoQ = model->ElectricChargeDensity();
164 const double rhoS = model->StrangenessDensity();
165 const double trace = (e != 0.0) ? (1.0 / 3.0 - p / e) : 0.0;
166 const double eps = 1.e-12;
167 const bool hasB = (std::abs(rhoB) > eps);
168 const bool hasQ = (std::abs(rhoQ) > eps);
169 const bool hasS = (std::abs(rhoS) > eps);
170
171 const auto safe_cT2 = [&](bool b, bool q, bool ss, double fallback) {
172 const double v = model->cT2(b, q, ss, false);
173 return std::isfinite(v) ? v : fallback;
174 };
175 const auto safe_cs2 = [&](bool b, bool q, bool ss, double fallback) {
176 const double v = model->cs2(b, q, ss, false);
177 return std::isfinite(v) ? v : fallback;
178 };
179
180 double cT2_B = 0.0, cs2_B = 0.0;
181 double cT2_BQ = 0.0, cs2_BQ = 0.0;
182 double cT2_BS = 0.0, cs2_BS = 0.0;
183 double cT2_BQS = 0.0, cs2_BQS = 0.0;
184
185 if (hasB) {
186 cT2_B = safe_cT2(true, false, false, 0.0);
187 cs2_B = safe_cs2(true, false, false, 0.0);
188
189 // At T=0, if extra conserved charges vanish, constraints reduce to B-only.
190 cT2_BQ = hasQ ? safe_cT2(true, true, false, cT2_B) : cT2_B;
191 cs2_BQ = hasQ ? safe_cs2(true, true, false, cs2_B) : cs2_B;
192 cT2_BS = hasS ? safe_cT2(true, false, true, cT2_B) : cT2_B;
193 cs2_BS = hasS ? safe_cs2(true, false, true, cs2_B) : cs2_B;
194 cT2_BQS = (hasQ && hasS) ? safe_cT2(true, true, true, cT2_B) : (hasS ? cT2_BS : cT2_BQ);
195 cs2_BQS = (hasQ && hasS) ? safe_cs2(true, true, true, cs2_B) : (hasS ? cs2_BS : cs2_BQ);
196 }
197
198 // At exactly T=0 these heat capacities vanish by definition.
199 const double cV_B = 0.0;
200 const double cV_BQ = 0.0;
201 const double cV_BS = 0.0;
202 const double cV_BQS = 0.0;
203
204 const double cMu = model->HeatCapacityMu();
205 double dsdT_mu_an = model->CalculateEntropyDensityDerivativeT();
206
207 // Numerical (ds/dT)_mu cross-check from one-sided derivatives at several dT.
208 // For a degenerate Fermi system, corrections are ~T^2 (Sommerfeld), so fit vs dT^2.
209 const std::vector<double> dTscan = {1.e-4, 1.e-3, 1.e-2};
210 std::vector<double> dsdT_scan, dsdT_scan_an;
211 dsdT_scan.reserve(dTscan.size());
212 for (double dTnum : dTscan) {
213 model->SetTemperature(dTnum);
214 model->CalculatePrimordialDensities();
215 const double sT = model->EntropyDensity();
216 dsdT_scan.push_back((sT - s) / dTnum);
217
218 double cVmu = model->HeatCapacityMu();
219 dsdT_scan_an.push_back(cVmu / dTnum);
220 }
221
222 const double dsdT_mu_num_1e4 = dsdT_scan[0];
223 const double dsdT_mu_num_1e3 = dsdT_scan[1];
224 const double dsdT_mu_num_1e2 = dsdT_scan[2];
225
226 const double dsdT_mu_an_1e4 = dsdT_scan_an[0];
227 const double dsdT_mu_an_1e3 = dsdT_scan_an[1];
228 const double dsdT_mu_an_1e2 = dsdT_scan_an[2];
229
230
231
232 // Restore the baseline zero-temperature state.
233 model->SetTemperature(0.0);
234 model->CalculatePrimordialDensities();
235 model->CalculateFluctuations();
236 model->CalculateTemperatureDerivatives();
237
238 fout << setw(tab) << muB << " "
239 << setw(tab) << p << " "
240 << setw(tab) << e << " "
241 << setw(tab) << s << " "
242 << setw(tab) << rhoB << " "
243 << setw(tab) << rhoQ << " "
244 << setw(tab) << rhoS << " "
245 << setw(tab) << trace << " "
246 << setw(tab) << cT2_B << " "
247 << setw(tab) << cs2_B << " "
248 << setw(tab) << cT2_BQ << " "
249 << setw(tab) << cs2_BQ << " "
250 << setw(tab) << cT2_BS << " "
251 << setw(tab) << cs2_BS << " "
252 << setw(tab) << cT2_BQS << " "
253 << setw(tab) << cs2_BQS << " "
254 << setw(tab) << cV_B << " "
255 << setw(tab) << cV_BQ << " "
256 << setw(tab) << cV_BS << " "
257 << setw(tab) << cV_BQS << " "
258 << setw(tab) << cMu << " "
259 << setw(tab) << dsdT_mu_num_1e4 << " "
260 << setw(tab) << dsdT_mu_num_1e3 << " "
261 << setw(tab) << dsdT_mu_num_1e2 << " "
262 << setw(tab) << dsdT_mu_an_1e4 << " "
263 << setw(tab) << dsdT_mu_an_1e3 << " "
264 << setw(tab) << dsdT_mu_an_1e2 << " "
265 << setw(tab) << dsdT_mu_an << "\n";
266
267 ++iters;
268 }
269
270 fout.close();
271
272 const double wt2 = get_wall_time();
273 cout << "Output file: " << outputfile << endl;
274 cout << setw(30) << "Running time: " << (wt2 - wt1) << " s" << endl;
275 if (iters > 0)
276 cout << setw(30) << "Time per point: " << (wt2 - wt1) / iters << " s" << endl;
277
278 return 0;
279}
280
ExcludedVolumeHelper.h
Contains some functions to deal with excluded volumes.
main
int main()
Definition TestEMMDerivatives.cpp:57
thermalfist::ExcludedVolumeModelBase
Base class implementing the ideal gas.
Definition ExcludedVolumeModels.h:26
thermalfist::ExcludedVolumeModelCS
Class implementing auxiliary functions for the Carnahan-Starling excluded volume model.
Definition ExcludedVolumeModels.h:66
thermalfist::ExcludedVolumeModelCrosstermsGeneralized
Implementation of a crossterms generalized excluded volume model.
Definition ExcludedVolumeModelsMulti.h:528
thermalfist::MeanFieldModelMultiVDW
Implementation of the van der Waals mean field model for multiple components.
Definition MeanFieldModelsMulti.h:139
thermalfist::ThermalModelRealGas
Class implementing the quantum real gas HRG model.
Definition ThermalModelRealGas.h:25
thermalfist::ThermalModelVDW
Class implementing the quantum van der Waals HRG model.
Definition ThermalModelVDW.h:37
thermalfist::ThermalParticleSystem
Class containing the particle list.
Definition ThermalParticleSystem.h:30
thermalfist
The main namespace where all classes and functions of the Thermal-FIST library reside.
Definition CosmicEoS.h:9
thermalfist::GetBaryonBaryonInteractionMatrix
std::vector< std::vector< double > > GetBaryonBaryonInteractionMatrix(const ThermalParticleSystem *TPS, double param)
Returns the matrix of attraction and repulsion parameters for baryon-baryon and antibaryon-antibaryon...
Definition ExcludedVolumeHelper.cpp:45
thermalfist::get_wall_time
double get_wall_time()
Definition Utility.cpp:235
thermalfist::SetVDWHRGInteractionParameters
void SetVDWHRGInteractionParameters(ThermalModelBase *model, double a, double b)
Sets vdW interactions for baryon-baryon and antibaryon-antibaryon pairs as in https://arxiv....
Definition ThermalModelVDW.cpp:1461
Generated by doxygen 1.13.2