Thermal-FIST 1.5
Package for hadron resonance gas model applications
Loading...
Searching...
No Matches
ThermodynamicsBQS.cpp
Go to the documentation of this file.
1/*
2 * Thermal-FIST package
3 *
4 * Copyright (c) 2025 Volodymyr Vovchenko
5 *
6 * GNU General Public License (GPLv3 or later)
7 */
8#include <string.h>
9#include <fstream>
10#include <iostream>
11#include <sstream>
12#include <iomanip>
13#include <ctime>
14#include <cstdio>
15#include <cassert>
16
17#include "HRGBase.h"
18#include "HRGEV.h"
19#include "HRGFit.h"
20#include "HRGVDW.h"
21#include "HRGRealGas.h"
22
24
25#include "ThermalFISTConfig.h"
26
27using namespace std;
28
29#ifdef ThermalFIST_USENAMESPACE
30using namespace thermalfist;
31#endif
32
33// The range in T, muB, muQ, muS to scan
34map<string, vector<double>> param_range = {
35 { "T", { 0.100, 0.180, 0.001 } }, // GeV
36 { "muB", { 0.000, 0.000, 0.01 } }, // GeV
37 { "muQ", { 0.000, 0.000, 0.01 } }, // GeV
38 { "muS", { 0.000, 0.000, 0.01 } } // GeV
39};
40
41// Read the parameter range from file
42void ReadParameterRangeFromFile(const std::string &filename)
43{
44 if (filename == "")
45 return;
46 std::ifstream fin(filename);
47 if (!fin.is_open()) {
48 std::cerr << "Error: cannot open file " << filename << " using default parameter range" << std::endl;
49 return;
50 }
51
52 std::string line;
53 while (std::getline(fin, line)) {
54 if (line.empty() || line[0] == '#') { continue; } // skip empty lines and comments
55 std::istringstream iss(line);
56 std::string param;
57 if (!(iss >> param)) continue;
58 if (param[0] == '#') { continue; } // skip comments
59 double min, max, step;
60 if (!(iss >> min >> max >> step)) { break; } // error
61 param_range[param] = { min, max, step };
62 }
63
64 fin.close();
65}
66
67// Calculates the HRG model equation of state properties
68// for a given range of T, muB, muQ, muS
69// Usage: ThermodynamicsBQS <param_range_file> <outputfile> <a> <b>
70// <a> -- parameter a for the QvdW model (GeV fm^3)
71// <b> -- parameter b for the QvdW model (fm^3)
72// useCS -- Use real gas model and Carnahan-Starling EV
73// <param_range_file> -- file with the parameter range
74// <outputfile> -- file to write the results to
75int main(int argc, char *argv[])
76{
77 // van der Waals attraction
78 double a = 0.;
79 if (argc > 1)
80 a = atof(argv[1]);
81
82 // van der Waals repulsion
83 double b = 0.;
84 if (argc > 2)
85 b = atof(argv[2]);
86
87 // Use real gas
88 bool useRG = false;
89 if (argc > 3)
90 useRG = atoi(argv[3]);
91
92 // Parameter range from file
93 string param_range_file = "";
94 if (argc > 4)
95 param_range_file = argv[4];
96 ReadParameterRangeFromFile(param_range_file);
97
98 // Model type
99 // 0 - Ideal HRG, 1 - EV HRG (no attraction), 2 - QvdW HRG (from 1609.03975), 3 - Real gas HRG
100 int ModelType = 0;
101 string ModelPrefix = "Id-HRG";
102 if (a == 0. && b == 0.) {
103 ModelType = 0;
104 ModelPrefix = "Id-HRG";
105 }
106 else if (a == 0.) {
107 ModelType = 1;
108 ModelPrefix = "EV-HRG";
109 }
110 else {
111 ModelType = 2;
112 ModelPrefix = "QvdW-HRG";
113 }
114
115 if (useRG) {
116 ModelType = 3;
117 ModelPrefix = "RG-HRG";
118 }
119
120 std::string outputfile = "Thermodynamics-" + ModelPrefix + "-output.dat";
121 if (argc > 5)
122 outputfile = argv[5];
123
124
125 // Create the hadron list instance and read the list from file
126 ThermalParticleSystem TPS(string(ThermalFIST_INPUT_FOLDER) + "/list/PDG2020/list.dat"); // <-- Default list, no light nuclei
127 //ThermalParticleSystem TPS(string(ThermalFIST_INPUT_FOLDER) + "/list/PDG2020/list-withnuclei.dat"); // <-- Default list, with light nuclei
128
129 // Create the ThermalModel instance
130 // Choose the class which fits the required variant of HRG model
131 ThermalModelBase *model;
132 assert(ModelType >= 0 && ModelType <= 3);
133 if (ModelType == 0) {
134 model = new ThermalModelIdeal(&TPS);
135 cout << "Performing calculation in the Id-HRG model..." << endl;
136 }
137 else if (ModelType == 1) {
138 model = new ThermalModelEVCrossterms(&TPS);
140 cout << "Performing calculation in the EV-HRG model..." << endl;
141 }
142 else if (ModelType == 2) {
143 model = new ThermalModelVDW(&TPS);
145 cout << "Performing calculation in the QvdW-HRG model..." << endl;
146 }
147 else if (ModelType == 3) {
148 model = new ThermalModelRealGas(&TPS);
149
150 // Excluded volume model
151 if (b != 0.) {
152 // ExcludedVolumeModelBase *evmod = new ExcludedVolumeModelVDW();
154 auto vdWb = GetBaryonBaryonInteractionMatrix(model->TPS(), b);
155 static_cast<ThermalModelRealGas *>(model)->SetExcludedVolumeModel(
157 }
158
159 // Mean field model
160 if (a != 0.) {
161 auto vdWa = GetBaryonBaryonInteractionMatrix(model->TPS(), a);
162 static_cast<ThermalModelRealGas *>(model)->SetMeanFieldModel(new MeanFieldModelMultiVDW(vdWa));
163 }
164 cout << "Performing calculation in the RG-HRG model..." << endl;
165 }
166
167 // Use (or not) finite resonance width
168 bool useWidth = false;
169 model->SetUseWidth(useWidth);
170 if (useWidth)
172
173 // Include (or not) quantum statistics
174 bool useQStats = true;
175 model->SetStatistics(useQStats);
176 if (useQStats)
178
179 // Prepare the output
180 ofstream fout(outputfile);
181 const int tabsize = 20;
182 fout << std::setw(tabsize) << "T[GeV]" << " ";
183 fout << std::setw(tabsize) << "muB[GeV]" << " ";
184 fout << std::setw(tabsize) << "muQ[GeV]" << " ";
185 fout << std::setw(tabsize) << "muS[GeV]" << " ";
186 fout << std::setw(tabsize) << "P[GeV/fm3]" << " ";
187 fout << std::setw(tabsize) << "e[GeV/fm3]" << " ";
188 fout << std::setw(tabsize) << "s[fm-3]" << " ";
189 fout << std::setw(tabsize) << "(1/3-p/e)" << " "; // Another defition of the trace anomaly
190 fout << std::setw(tabsize) << "rhoB[fm-3]" << " ";
191 fout << std::setw(tabsize) << "rhoQ[fm-3]" << " ";
192 fout << std::setw(tabsize) << "rhoS[fm-3]" << " ";
193 // Temperature normalized quantities (a la lattice QCD)
194 fout << std::setw(tabsize) << "P/T^4" << " ";
195 fout << std::setw(tabsize) << "e/T^4" << " ";
196 fout << std::setw(tabsize) << "s/T^3" << " ";
197 fout << std::setw(tabsize) << "(e-3p)/T^3" << " ";
198 fout << std::setw(tabsize) << std::endl;
199
200 // Timer
201 double wt1 = get_wall_time();
202 int iters = 0; // number of iterations (phase diagram points)
203
204 // Loop over the thermodynamic variables
205 double Tmin = param_range["T"][0];
206 double Tmax = param_range["T"][1];
207 double dT = param_range["T"][2];
208 for(double T = Tmin; T <= Tmax + 0.1 * dT; T += dT) {
209 model->SetTemperature(T);
210
211 // Output the outer loop
212 cout << "Calculating the temperature " << T << " GeV" << endl;
213
214 double muBmin = param_range["muB"][0];
215 double muBmax = param_range["muB"][1];
216 double dmuB = param_range["muB"][2];
217 for(double muB = muBmin; muB <= muBmax + 0.1 * dmuB; muB += dmuB) {
218 model->SetBaryonChemicalPotential(muB);
219
220 double muQmin = param_range["muQ"][0];
221 double muQmax = param_range["muQ"][1];
222 double dmuQ = param_range["muQ"][2];
223 for(double muQ = muQmin; muQ <= muQmax + 0.1 * dmuQ; muQ += dmuQ) {
225
226 double muSmin = param_range["muS"][0];
227 double muSmax = param_range["muS"][1];
228 double dmuS = param_range["muS"][2];
229 for(double muS = muSmin; muS <= muSmax + 0.1 * dmuS; muS += dmuS) {
231
232 // Perform the calculation
234 iters++;
235
236 // Collect the output
237 double p = model->Pressure(); // Pressure in GeV/fm3
238 double e = model->EnergyDensity(); // Energy density in GeV/fm3
239 double s = model->EntropyDensity(); // Entropy density in fm-3
240 double rhoB = model->BaryonDensity(); // Baryon density in fm-3
241 double rhoQ = model->ElectricChargeDensity(); // Electric charge density in fm-3
242 double rhoS = model->StrangenessDensity(); // Strangeness density in fm-3
243 double trace_anomaly = (1./3. - p/e);
244
245 double pT4 = p / pow(T, 4) / xMath::GeVtoifm3();
246 double eT4 = e / pow(T, 4) / xMath::GeVtoifm3();
247 double sT3 = s / pow(T, 3) / xMath::GeVtoifm3();
248 double IT4 = (e - 3 * p) / pow(T, 4) / xMath::GeVtoifm3();
249
250 // Print to file
251 fout << setw(tabsize) << T << " ";
252 fout << setw(tabsize) << muB << " ";
253 fout << setw(tabsize) << muQ << " ";
254 fout << setw(tabsize) << muS << " ";
255 fout << setw(tabsize) << p << " ";
256 fout << setw(tabsize) << e << " ";
257 fout << setw(tabsize) << s << " ";
258 fout << setw(tabsize) << trace_anomaly << " ";
259 fout << setw(tabsize) << rhoB << " ";
260 fout << setw(tabsize) << rhoQ << " ";
261 fout << setw(tabsize) << rhoS << " ";
262 fout << setw(tabsize) << pT4 << " ";
263 fout << setw(tabsize) << eT4 << " ";
264 fout << setw(tabsize) << sT3 << " ";
265 fout << setw(tabsize) << IT4 << " ";
266 fout << setw(tabsize) << endl;
267 }
268 }
269 }
270 }
271
272 fout.close();
273
274 double wt2 = get_wall_time();
275 cout << setw(30) << "Running time: " << (wt2 - wt1) << " s" << endl;
276 cout << setw(30) << "Time per single calculation: " << (wt2 - wt1) / iters << " s" << endl;
277
278 // Cleanup
279 delete model;
280
281 return 0;
282}
283
Contains some functions to deal with excluded volumes.
map< string, vector< double > > param_range
int main(int argc, char *argv[])
void ReadParameterRangeFromFile(const std::string &filename)
Base class implementing the ideal gas.
Class implementing auxiliary functions for the Carnahan-Starling excluded volume model.
Implementation of a crossterms generalized excluded volume model.
Implementation of the van der Waals mean field model for multiple components.
Abstract base class for an HRG model implementation.
virtual void SetTemperature(double T)
Set the temperature.
void SetUseWidth(bool useWidth)
Sets whether finite resonance widths are used. Deprecated.
virtual void SetBaryonChemicalPotential(double muB)
Set the baryon chemical potential.
virtual void SetElectricChemicalPotential(double muQ)
Set the electric chemical potential.
virtual void CalculatePrimordialDensities()=0
Calculates the primordial densities of all species.
virtual void SetStrangenessChemicalPotential(double muS)
Set the strangeness chemical potential.
virtual void SetCalculationType(IdealGasFunctions::QStatsCalculationType type)
Sets the CalculationType() method to evaluate quantum statistics. Calls the corresponding method in T...
virtual void SetStatistics(bool stats)
Class implementing the crossterms excluded-volume model.
Class implementing the Ideal HRG model.
Class implementing the quantum real gas HRG model.
Class implementing the quantum van der Waals HRG model.
@ eBWconstBR
Energy-dependent Breit-Wigner scheme (eBW) with constant branching ratios when evaluating feeddown.
Class containing the particle list.
constexpr double GeVtoifm3()
A constant to transform GeV into fm .
Definition xMath.h:31
The main namespace where all classes and functions of the Thermal-FIST library reside.
Definition CosmicEoS.h:9
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...
double get_wall_time()
Definition Utility.cpp:191
void SetVDWHRGInteractionParameters(ThermalModelBase *model, double a, double b)
Sets vdW interactions for baryon-baryon and antibaryon-antibaryon pairs as in https://arxiv....
void SetEVHRGInteractionParameters(ThermalModelBase *model, double b)
Sets EV interactions for baryon-baryon and antibaryon-antibaryon pairs as in https://arxiv....