PSCF v1.3
rpg/fts/montecarlo/McMove.tpp
1#ifndef RPG_MC_MOVE_TPP
2#define RPG_MC_MOVE_TPP
3
4/*
5* PSCF - Polymer Self-Consistent Field
6*
7* Copyright 2015 - 2025, The Regents of the University of Minnesota
8* Distributed under the terms of the GNU General Public License.
9*/
10
11#include "McMove.h"
12
13#include <rpg/fts/montecarlo/McSimulator.h>
14#include <rpg/fts/compressor/Compressor.h>
15#include <rpg/system/System.h>
16
17namespace Pscf {
18namespace Rpg {
19
20 using namespace Util;
21
22 /*
23 * Constructor.
24 */
25 template <int D>
27 : simulatorPtr_(&simulator),
28 systemPtr_(&(simulator.system())),
29 randomPtr_(&(simulator.random())),
30 cudaRandomPtr_(&(simulator.cudaRandom())),
31 nAttempt_(0),
32 nAccept_(0)
33 {}
34
35 /*
36 * Destructor, empty default implementation.
37 */
38 template <int D>
41
42 /*
43 * ReadParameters, empty default implementation.
44 */
45 template <int D>
46 void McMove<D>::readParameters(std::istream &in)
47 {}
48
49 /*
50 * Read the probability from file.
51 */
52 template <int D>
53 void McMove<D>::readProbability(std::istream &in)
54 { read<double>(in, "probability", probability_); }
55
56 /*
57 * Setup at beginning of loop.
58 *
59 * Trivial default implementation - initializes counters.
60 */
61 template <int D>
63 {
64 nAttempt_ = 0;
65 nAccept_ = 0;
66 nFail_ = 0;
68 simulator().computeWc();
69 if (simulator().needsCc() || simulator().needsDc()){
70 system().compute();
71 simulator().computeCc();
72 }
73 if (simulator().needsDc()){
74 simulator().computeDc();
75 }
76 }
77
78 /*
79 * Trivial default implementation - always returns false.
80 */
81 template <int D>
83 {
84 totalTimer_.start();
86
87 // Get current Hamiltonian
88 double oldHamiltonian = simulator().hamiltonian();
89
90 // Save current state
91 simulator().saveState();
92
93 // Clear both eigen-components of the fields and hamiltonian
94 simulator().clearData();
95
96 // Attempt modification
97 attemptMoveTimer_.start();
99 attemptMoveTimer_.stop();
100
101 // Call compressor
102 compressorTimer_.start();
103 int compress = simulator().compressor().compress();
104 UTIL_CHECK(system().c().hasData());
105 compressorTimer_.stop();
106
107 bool isConverged = false;
108 if (compress != 0){
110 simulator().restoreState();
111 } else {
112 isConverged = true;
113
114 // Compute eigenvector components of the current w fields
115 componentTimer_.start();
116 simulator().computeWc();
117 // Compute cc fields if any move require cc fields
118 if (simulator().needsCc() || simulator().needsDc()){
119 //system().compute();
120 UTIL_CHECK(system().c().hasData());
121 simulator().computeCc();
122 }
123 // Compute dc fields if any move require dc fields
124 if (simulator().needsDc()){
125 simulator().computeDc();
126 }
127 componentTimer_.stop();
128
129 // Evaluate new Hamiltonian
130 hamiltonianTimer_.start();
131 simulator().computeHamiltonian();
132 double newHamiltonian = simulator().hamiltonian();
133 hamiltonianTimer_.stop();
134
135 // Accept or reject move
136 bool accept = false;
137 decisionTimer_.start();
138 double weight = exp(-(newHamiltonian - oldHamiltonian));
139 accept = random().metropolis(weight);
140 if (accept) {
142 simulator().clearState();
143 } else {
144 simulator().restoreState();
145 }
146 decisionTimer_.stop();
147
148 }
149 totalTimer_.stop();
150
151 return isConverged;
152 }
153
154 /*
155 * Trivial default implementation - do nothing
156 */
157 template <int D>
159 {}
160
161 template<int D>
162 void McMove<D>::outputTimers(std::ostream& out)
163 {
164 // Output timing results, if requested.
165 double total = totalTimer_.time();
166 out << " "
167 << "Total" << std::setw(17) << "Per Move" << std::setw(14) << "Fraction" << "\n";
168 out << "Attempt Move: "
169 << Dbl(attemptMoveTimer_.time(), 9, 3) << " s, "
170 << Dbl(attemptMoveTimer_.time()/nAttempt_, 9, 3) << " s, "
171 << Dbl(attemptMoveTimer_.time()/total, 9, 3) << "\n";
172 out << "Compressor: "
173 << Dbl(compressorTimer_.time(), 9, 3) << " s, "
174 << Dbl(compressorTimer_.time()/nAttempt_, 9, 3) << " s, "
175 << Dbl(compressorTimer_.time()/total, 9, 3) << "\n";
176 out << "Compute eigen-components: "
177 << Dbl(componentTimer_.time(), 9, 3) << " s, "
178 << Dbl(componentTimer_.time()/nAttempt_, 9, 3) << " s, "
179 << Dbl(componentTimer_.time()/total, 9, 3) << "\n";
180 out << "Compute Hamiltonian: "
181 << Dbl(hamiltonianTimer_.time(), 9, 3) << " s, "
182 << Dbl(hamiltonianTimer_.time()/nAttempt_, 9, 3) << " s, "
183 << Dbl(hamiltonianTimer_.time()/total, 9, 3) << "\n";
184 out << "Accept or Reject: "
185 << Dbl(decisionTimer_.time(), 9, 3) << " s, "
186 << Dbl(decisionTimer_.time()/nAttempt_, 9, 3) << " s, "
187 << Dbl(decisionTimer_.time()/total, 9, 3) << "\n";
188 out << "total time: "
189 << Dbl(total, 9, 3) << " s, "
190 << Dbl(total/nAttempt_, 9, 3) << " s \n";
191 out << "\n";
192 }
193
194 template<int D>
196 {
197 attemptMoveTimer_.clear();
198 compressorTimer_.clear();
199 componentTimer_.clear();
200 hamiltonianTimer_.clear();
201 decisionTimer_.clear();
202 totalTimer_.clear();
203 }
204
205
206}
207}
208#endif
CudaRandom & cudaRandom()
Get cuda random number generator by reference.
virtual ~McMove()
Destructor.
void incrementNAttempt()
Increment the number of attempted moves.
Random & random()
Get Random number generator of parent System.
virtual void attemptMove()
Attempt unconstrained move.
virtual bool move()
Generate, attempt, and accept or reject a Monte Carlo move.
void readProbability(std::istream &in)
Read the probability from file.
virtual bool needsCc()
Decide whether cc fields need to be saved for move The default implementation is false.
McSimulator< D > & simulator()
Get parent McSimulator object.
virtual void readParameters(std::istream &in)
Read required parameters from file.
System< D > & system()
Get parent System object.
void incrementNAccept()
Increment the number of accepted moves.
void incrementNFail()
Increment the number of failed moves.
virtual void outputTimers(std::ostream &out)
Log output timing results.
virtual void setup()
Setup before the beginning of each simulation run.
McMove(McSimulator< D > &simulator)
Constructor.
Timer attemptMoveTimer_
Timers for McMove.
virtual bool needsDc()
Decide whether dc fields need to be saved for move The default implementation is false.
virtual void clearTimers()
Clear timers.
virtual void output()
Output statistics for this move (at the end of simulation)
Monte-Carlo simulation coordinator.
Wrapper for a double precision number, for formatted ostream output.
Definition Dbl.h:40
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.
Definition param_pc.dox:1