1#ifndef RPG_MC_SIMULATOR_TPP
2#define RPG_MC_SIMULATOR_TPP
11#include "McSimulator.h"
13#include <rpg/System.h>
14#include <rpg/fts/montecarlo/McMoveFactory.h>
15#include <rpg/fts/analyzer/AnalyzerFactory.h>
16#include <rpg/fts/trajectory/TrajectoryReader.h>
17#include <rpg/fts/trajectory/TrajectoryReaderFactory.h>
18#include <rpg/fts/compressor/Compressor.h>
19#include <rpg/fts/perturbation/Perturbation.h>
20#include <rpg/fts/perturbation/PerturbationFactory.h>
21#include <rpg/fts/ramp/RampFactory.h>
22#include <rpg/fts/ramp/Ramp.h>
24#include <util/random/Random.h>
25#include <pscf/cuda/CudaRandom.h>
26#include <util/misc/Timer.h>
29#include <gsl/gsl_eigen.h>
43 mcMoveManager_(*this, system),
44 analyzerManager_(*this, system),
45 trajectoryReaderFactoryPtr_(0)
48 trajectoryReaderFactoryPtr_
58 if (trajectoryReaderFactoryPtr_) {
59 delete trajectoryReaderFactoryPtr_;
73 readParamCompositeOptional(in, mcMoveManager_);
77 readCompressor(in, isEnd);
78 readPerturbation(in, isEnd);
88 readParamCompositeOptional(in, analyzerManager_);
91 state_.needsCc =
false;
92 state_.needsDc =
false;
93 state_.needsHamiltonian =
true;
94 if (mcMoveManager_.needsCc()){
95 state_.needsCc =
true;
97 if (mcMoveManager_.needsDc()){
98 state_.needsDc =
true;
116 if (hasPerturbation()) {
117 perturbation().setup();
128 if (hasCompressor()) {
129 compressor().compress();
130 compressor().clearTimers();
135 if (state_.needsCc || state_.needsDc) {
138 if (state_.needsDc) {
141 computeHamiltonian();
145 mcMoveManager_.setup();
149 if (analyzerManager_.size() > 0){
150 analyzerManager_.setup();
175 ramp().setParameters(iStep_);
179 analyzerTimer.
start();
180 analyzerManager_.sample(iStep_);
181 analyzerTimer.
stop();
183 for (iTotalStep_ = 0; iTotalStep_ < nStep; ++iTotalStep_) {
186 converged = mcMoveManager_.chooseMove().move();
192 ramp().setParameters(iStep_);
196 analyzerTimer.
start();
199 if (analyzerManager_.size() > 0) {
200 analyzerManager_.sample(iStep_);
204 analyzerTimer.
stop();
207 Log::file() <<
"Step: "<< iTotalStep_<<
" fail to converge" <<
"\n";
212 double time = timer.
time();
213 double analyzerTime = analyzerTimer.
time();
216 mcMoveManager_.output();
218 analyzerManager_.output();
228 Log::file() <<
"nStep " << nStep << std::endl;
229 if (iStep_ != nStep){
230 Log::file() <<
"nFail Step " << (nStep - iStep_) << std::endl;
233 <<
" sec" << std::endl;
234 double rStep = double(nStep);
235 Log::file() <<
"time / nStep " << time / rStep
236 <<
" sec" << std::endl;
237 Log::file() <<
"Analyzer run time " << analyzerTime
238 <<
" sec" << std::endl;
249 Log::file() <<
"Move Statistics:" << endl << endl;
250 Log::file() << setw(20) << left <<
"Move Name"
251 << setw(10) << right <<
"Attempted"
252 << setw(10) << right <<
"Accepted"
253 << setw(13) << right <<
"AcceptRate"
254 << setw(10) << right <<
"Failed"
255 << setw(13) << right <<
"FailRate"
257 int nMove = mcMoveManager_.size();
258 for (
int iMove = 0; iMove < nMove; ++iMove) {
259 attempt = mcMoveManager_[iMove].nAttempt();
260 accept = mcMoveManager_[iMove].nAccept();
261 fail = mcMoveManager_[iMove].nFail();
263 << mcMoveManager_[iMove].className()
264 << setw(10) << right << attempt
265 << setw(10) << accept
266 << setw(13) << fixed << setprecision(5)
267 << ( attempt == 0 ? 0.0 : double(accept)/double(attempt) )
269 << setw(13) << fixed << setprecision(5)
270 << ( attempt == 0 ? 0.0 : double(fail)/double(attempt) )
282 std::string classname,
283 std::string filename)
293 trajectoryReaderPtr = trajectoryReaderFactory().factory(classname);
294 if (!trajectoryReaderPtr) {
296 message =
"Invalid TrajectoryReader class name " + classname;
301 trajectoryReaderPtr->
open(filename);
302 trajectoryReaderPtr->readHeader();
308 hasFrame = trajectoryReaderPtr->
readFrame();
310 for (iStep_ = 0; iStep_ <= max && hasFrame; ++iStep_) {
322 analyzerManager_.sample(iStep_);
326 hasFrame = trajectoryReaderPtr->
readFrame();
329 Log::file() <<
"end main loop" << std::endl;
330 int nFrames = iStep_ - min;
331 trajectoryReaderPtr->
close();
332 delete trajectoryReaderPtr;
335 analyzerManager_.output();
339 Log::file() <<
"# of frames " << nFrames << std::endl;
341 <<
" sec" << std::endl;
342 Log::file() <<
"time / frame " << timer.
time()/double(nFrames)
343 <<
" sec" << std::endl;
354 compressor().outputTimers(out);
356 out <<
"MC move time contributions:\n";
357 mcMoveManager_.outputTimers(out);
365 { mcMoveManager_.clearTimers(); }
Abstract base for periodic output and/or analysis actions.
Monte-Carlo simulation coordinator.
virtual void outputTimers(std::ostream &out)
Output timing results.
McSimulator(System< D > &system)
Constructor.
void setClassName(const char *className)
Set class name string.
virtual void analyze(int min, int max, std::string classname, std::string filename)
Read and analyze a trajectory file.
void simulate(int nStep)
Perform a field theoretic Monte-Carlo simulation.
virtual void clearTimers()
Clear timers.
virtual void readParameters(std::istream &in)
Read parameters for a MC simulation.
~McSimulator()
Destructor.
Field theoretic simulator (base class).
System< D > & system()
Get parent system by reference.
Main class for calculations that represent one system.
Factory for subclasses of TrajectoryReader.
virtual bool readFrame()=0
Read a single frame.
virtual void close()=0
Close the trajectory file.
virtual void open(std::string filename)=0
Open trajectory file and read header, if any.
static std::ostream & file()
Get log ostream by reference.
void start(TimePoint begin)
Start timing from an externally supplied time.
double time()
Return the accumulated time, in seconds.
void stop(TimePoint end)
Stop the clock at an externally supplied time.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Fields, FFTs, and utilities for periodic boundary conditions (CUDA)
PSCF package top-level namespace.
Utility classes for scientific computation.