Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
hmc_General.cpp
Go to the documentation of this file.
1 
14 #include "hmc_General.h"
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 //- parameter entries
21 namespace {
22  void append_entry(Parameters& param)
23  {
24  param.Register_int("Metropolis_test", 0);
25 
26  param.Register_string("verbose_level", "NULL");
27  }
28 
29 
30 #ifdef USE_PARAMETERS_FACTORY
31  bool init_param = ParametersFactory::Register("HMC.General", append_entry);
32 #endif
33 }
34 //- end
35 
36 //- parameters class
37 Parameters_HMC_General::Parameters_HMC_General() { append_entry(*this); }
38 //- end
39 
40 //====================================================================
41 void HMC_General::set_parameters(const Parameters& params)
42 {
43  const string str_vlevel = params.get_string("verbose_level");
44 
45  m_vl = vout.set_verbose_level(str_vlevel);
46 
47  //- fetch and check input parameters
48  int Mtpl_test;
49 
50  int err = 0;
51  err += params.fetch_int("Metropolis_test", Mtpl_test);
52 
53  if (err) {
54  vout.crucial(m_vl, "HMC_General: fetch error, input parameter not found.\n");
55  abort();
56  }
57 
58 
59  set_parameters(Mtpl_test);
60 }
61 
62 
63 //====================================================================
64 void HMC_General::set_parameters(int Mtpl_test)
65 {
66  //- print input parameters
67  vout.general(m_vl, "HMC (general) setup:\n");
68  vout.general(m_vl, " Number of actions: %4d\n", m_action.size());
69  vout.general(m_vl, " Mtlp_test = %4d\n", Mtpl_test);
70 
71  //- range check
72  // NB. Mtpl_test == 0 is allowed.
73 
74  //- store values
75  m_Mtpl_test = Mtpl_test;
76 }
77 
78 
79 //====================================================================
80 double HMC_General::update(Field_G& conf_org)
81 {
82  int Nc = CommonParameters::Nc();
83  int Lvol = CommonParameters::Lvol();
84 
85  Field_G U(conf_org);
86 
87  for (int i = 0; i < m_action.size(); ++i) {
88  m_action[i]->set_config(&U);
89  }
90 
91  int Nin = U.nin();
92  int Nvol = U.nvol();
93  int Nex = U.nex();
94 
95  Field_G iP(Nvol, Nex);
96 
97 
98  vout.general(m_vl, "HMC (general) start.\n");
99 
100  // Langevin step
101  double H_tot0 = langevin(iP, U);
102  vout.general(m_vl, " H_tot(init) = %20.10f\n", H_tot0);
103 
104  double plaq0 = m_staple->plaquette(U);
105  vout.general(m_vl, " plaq(init) = %20.12f\n", plaq0);
106 
107  // initial Hamiltonian
108  // H_tot0 = calc_Hamiltonian(iP,U);
109 
110  // molecular dynamical integration
111  m_integrator->evolve(iP, U);
112 
113  // final Hamiltonian
114  double H_tot1 = calc_Hamiltonian(iP, U);
115  vout.general(m_vl, " H_tot(fin) = %20.10f\n", H_tot1);
116 
117  double plaq1 = m_staple->plaquette(U);
118  vout.general(m_vl, " plaq(fin) = %20.12f\n", plaq1);
119 
120 
121  // Metropolis test
122  vout.general(m_vl, "Metropolis test.\n");
123 
124  double diff_H = H_tot1 - H_tot0;
125  double Arprob = exp(-diff_H);
126  double rn = -1.0;
127  if (m_Mtpl_test != 0) {
128  rn = m_rand->get();
129  }
130  vout.general(m_vl, " H(diff) = %14.8f\n", diff_H);
131  vout.general(m_vl, " Acceptance = %14.8f\n", Arprob);
132  vout.general(m_vl, " Random number = %14.8f\n", rn);
133 
134  if (rn <= Arprob) { // accepted
135  vout.general(m_vl, " Accepted\n");
136  conf_org = U;
137  } else { // rejected
138  vout.general(m_vl, " Rejected\n");
139  }
140 
141  return m_staple->plaquette(U);
142 }
143 
144 
145 //====================================================================
146 double HMC_General::langevin(Field_G& iP, Field_G& U)
147 {
148  int Nc = CommonParameters::Nc();
149  int Lvol = CommonParameters::Lvol();
150  int Ndim = CommonParameters::Ndim();
151  int NcA = Nc * Nc - 1;
152 
153  int size_iP = NcA * Lvol * Ndim;
154 
155  vout.general(m_vl, "Langevin step.\n");
156 
157  for (int i = 0; i < m_action.size(); ++i) {
158  m_action[i]->notify_linkv();
159  }
160 
161  for (int i = 0; i < m_director.size(); ++i) {
162  m_director[i]->notify_linkv();
163  }
164 
165  // kinetic term
166  double H_iP = langevin_P(iP);
167  vout.general(m_vl, " Kinetic term:\n");
168  vout.general(m_vl, " H_kin = %18.8f\n", H_iP);
169  vout.general(m_vl, " H_kin/dof = %18.8f\n", H_iP / size_iP);
170 
171 
172  double H_actions = 0.0;
173  for (int i = 0; i < m_action.size(); ++i) {
174  H_actions += m_action[i]->langevin(m_rand);
175  }
176 
177  double H_tot = H_iP + H_actions;
178  return H_tot;
179 }
180 
181 
182 //====================================================================
183 double HMC_General::langevin_P(Field_G& iP)
184 {
185  int Nin = iP.nin();
186  int Nvol = iP.nvol();
187  int Nex = iP.nex();
188 
189  int Nc = CommonParameters::Nc();
190  int NcA = Nc * Nc - 1;
191 
192  // confirm that now gauge group is SU(3)
193  assert(NcA == 8);
194 
195  Field Hrand(NcA, Nvol, Nex);
196 
197  m_rand->gauss_lex_global(Hrand);
198 
199  double sq3r = 1.0 / sqrt(3.0);
200  double sq3r2 = 2.0 * sq3r;
201 
202  for (int ex = 0; ex < Nex; ++ex) {
203  for (int site = 0; site < Nvol; ++site) {
204  //here SU(3) is explicitly assumed. need to generalize.
205  double hc1 = Hrand.cmp(0, site, ex);
206  double hc2 = Hrand.cmp(1, site, ex);
207  double hc3 = Hrand.cmp(2, site, ex);
208  double hc4 = Hrand.cmp(3, site, ex);
209  double hc5 = Hrand.cmp(4, site, ex);
210  double hc6 = Hrand.cmp(5, site, ex);
211  double hc7 = Hrand.cmp(6, site, ex);
212  double hc8 = Hrand.cmp(7, site, ex);
213  iP.set(0, site, ex, 0.0);
214  iP.set(1, site, ex, hc3 + hc8 * sq3r);
215  iP.set(2, site, ex, hc2);
216  iP.set(3, site, ex, hc1);
217  iP.set(4, site, ex, hc5);
218  iP.set(5, site, ex, hc4);
219  iP.set(6, site, ex, -hc2);
220  iP.set(7, site, ex, hc1);
221  iP.set(8, site, ex, 0.0);
222  iP.set(9, site, ex, -hc3 + hc8 * sq3r);
223  iP.set(10, site, ex, hc7);
224  iP.set(11, site, ex, hc6);
225  iP.set(12, site, ex, -hc5);
226  iP.set(13, site, ex, hc4);
227  iP.set(14, site, ex, -hc7);
228  iP.set(15, site, ex, hc6);
229  iP.set(16, site, ex, 0.0);
230  iP.set(17, site, ex, -hc8 * sq3r2);
231  }
232  }
233 
234  double iP2 = iP.norm();
235  iP2 = 0.5 * iP2 * iP2;
236 
237  return iP2;
238 }
239 
240 
241 //====================================================================
242 double HMC_General::calc_Hamiltonian(Field_G& iP, Field_G& U)
243 {
244  int Nin = U.nin();
245  int Nvol = U.nvol();
246  int Nex = U.nex();
247 
248  int Nc = CommonParameters::Nc();
249  int Nd = CommonParameters::Nd();
250  int Lvol = CommonParameters::Lvol();
251  int NcA = Nc * Nc - 1;
252 
253  int size_iP = NcA * Lvol * Nex;
254  int size_U = Lvol * 6.0;
255  int size_psf = Lvol * Nc * Nd;
256 
257 
258  vout.general(m_vl, "Hamiltonian calculation.\n");
259 
260  // kinetic term
261  double H_iP = calcH_P(iP);
262  vout.general(m_vl, " Kinetic term:\n");
263  vout.general(m_vl, " H_kin = %18.8f\n", H_iP);
264  vout.general(m_vl, " H_kin/dof = %18.8f\n", H_iP / size_iP);
265 
266  double H_actions = 0.0;
267  for (int i = 0; i < m_action.size(); ++i) {
268  H_actions += m_action[i]->calcH();
269  }
270 
271  double H_tot = H_iP + H_actions;
272  return H_tot;
273 }
274 
275 
276 //====================================================================
277 double HMC_General::calcH_P(Field_G& iP)
278 {
279  double hn = iP.norm();
280  double H_iP = 0.5 * hn * hn;
281 
282  return H_iP;
283 }
284 
285 
286 //====================================================================
287 //============================================================END=====