PSCF v1.3
rpc/fts/montecarlo/McMove.tpp
1#ifndef RPC_MC_MOVE_TPP
2#define RPC_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 <rpc/fts/montecarlo/McSimulator.h>
14#include <rpc/fts/compressor/Compressor.h>
15#include <rpc/system/System.h>
16#include <util/archives/Serializable_includes.h>
17
18namespace Pscf {
19namespace Rpc {
20
21 using namespace Util;
22
23 /*
24 * Constructor.
25 */
26 template <int D>
28 : simulatorPtr_(&simulator),
29 systemPtr_(&(simulator.system())),
30 randomPtr_(&(simulator.random())),
31 nAttempt_(0),
32 nAccept_(0),
33 nFail_(0)
34 {}
35
36 /*
37 * Destructor, empty default implementation.
38 */
39 template <int D>
42
43 /*
44 * ReadParameters, empty default implementation.
45 */
46 template <int D>
47 void McMove<D>::readParameters(std::istream &in)
48 {}
49
50 /*
51 * Read the probability from file.
52 */
53 template <int D>
54 void McMove<D>::readProbability(std::istream &in)
55 { read<double>(in, "probability", probability_); }
56
57 /*
58 * Setup at beginning of loop.
59 *
60 * Trivial default implementation - initializes counters.
61 */
62 template <int D>
64 {
65 nAttempt_ = 0;
66 nAccept_ = 0;
67 nFail_ = 0;
69 simulator().computeWc();
70
71 if (simulator().needsCc() || simulator().needsDc()){
72 system().compute();
73 simulator().computeCc();
74 }
75
76 if (simulator().needsDc()){
77 simulator().computeDc();
78 }
79 }
80
81 /*
82 * Trivial default implementation - always returns false.
83 */
84 template <int D>
86 {
87 // Start timers
88 totalTimer_.start();
89 attemptMoveTimer_.start();
91
92 // Get current Hamiltonian
93 double oldHamiltonian = simulator().hamiltonian();
94
95 // Save current state
96 simulator().saveState();
97
98 // Clear both eigen-components of the fields and hamiltonian
99 simulator().clearData();
100
101 // Attempt modification of exchange fields
102 attemptMove();
103 attemptMoveTimer_.stop();
104
105 // Call compressor to modify pressure-like fields
106 compressorTimer_.start();
107 int compress = simulator().compressor().compress();
108 UTIL_CHECK(system().c().hasData());
109 compressorTimer_.stop();
110
111 bool isConverged = false;
112 if (compress != 0){
114 simulator().restoreState();
115 } else {
116 isConverged = true;
117
118 // Compute eigenvector components of the current fields
119 componentTimer_.start();
120 simulator().computeWc();
121 // Compute cc fields if any move require cc fields
122 if (simulator().needsCc() || simulator().needsDc()){
123 UTIL_CHECK(system().c().hasData());
124 //system().compute();
125 simulator().computeCc();
126 }
127 // Compute dc fields if any move require dc fields
128 if (simulator().needsDc()){
129 simulator().computeDc();
130 }
131 componentTimer_.stop();
132
133 // Evaluate new Hamiltonian
134 hamiltonianTimer_.start();
135 simulator().computeHamiltonian();
136 double newHamiltonian = simulator().hamiltonian();
137 hamiltonianTimer_.stop();
138
139 // Accept or reject move
140 decisionTimer_.start();
141 bool accept = false;
142 double weight = exp(-(newHamiltonian - oldHamiltonian));
143 accept = random().metropolis(weight);
144 if (accept) {
146 simulator().clearState();
147 } else {
148 simulator().restoreState();
149 }
150 decisionTimer_.stop();
151
152 }
153 totalTimer_.stop();
154
155 return isConverged;
156 }
157
158 /*
159 * Trivial default implementation - do nothing
160 */
161 template <int D>
163 {}
164
165 template<int D>
166 void McMove<D>::outputTimers(std::ostream& out)
167 {
168 // Output timing results, if requested.
169 double total = totalTimer_.time();
170 out << " "
171 << "Total" << std::setw(17) << "Per Move" << std::setw(14) << "Fraction" << "\n";
172 out << "Attempt Move: "
173 << Dbl(attemptMoveTimer_.time(), 9, 3) << " s, "
174 << Dbl(attemptMoveTimer_.time()/nAttempt_, 9, 3) << " s, "
175 << Dbl(attemptMoveTimer_.time()/total, 9, 3) << "\n";
176 out << "Compressor: "
177 << Dbl(compressorTimer_.time(), 9, 3) << " s, "
178 << Dbl(compressorTimer_.time()/nAttempt_, 9, 3) << " s, "
179 << Dbl(compressorTimer_.time()/total, 9, 3) << "\n";
180 out << "Compute eigen-components: "
181 << Dbl(componentTimer_.time(), 9, 3) << " s, "
182 << Dbl(componentTimer_.time()/nAttempt_, 9, 3) << " s, "
183 << Dbl(componentTimer_.time()/total, 9, 3) << "\n";
184 out << "Compute Hamiltonian: "
185 << Dbl(hamiltonianTimer_.time(), 9, 3) << " s, "
186 << Dbl(hamiltonianTimer_.time()/nAttempt_, 9, 3) << " s, "
187 << Dbl(hamiltonianTimer_.time()/total, 9, 3) << "\n";
188 out << "Accept or Reject: "
189 << Dbl(decisionTimer_.time(), 9, 3) << " s, "
190 << Dbl(decisionTimer_.time()/nAttempt_, 9, 3) << " s, "
191 << Dbl(decisionTimer_.time()/total, 9, 3) << "\n";
192 out << "total time: "
193 << Dbl(total, 9, 3) << " s, "
194 << Dbl(total/nAttempt_, 9, 3) << " s \n";
195 out << "\n";
196 }
197
198 template<int D>
200 {
201 attemptMoveTimer_.clear();
202 compressorTimer_.clear();
203 componentTimer_.clear();
204 hamiltonianTimer_.clear();
205 decisionTimer_.clear();
206 totalTimer_.clear();
207 }
208
209}
210}
211#endif
void incrementNAccept()
Increment the number of accepted moves.
virtual void clearTimers()
Clear timers.
Timer attemptMoveTimer_
Timers for McMove.
virtual bool needsDc()
Decide whether dc fields need to be saved for move The default implementation is false.
Random & random()
Get Random number generator of parent System.
McMove(McSimulator< D > &simulator)
Constructor.
virtual bool needsCc()
Decide whether cc fields need to be saved for move The default implementation is false.
virtual void readParameters(std::istream &in)
Read required parameters from file.
System< D > & system()
Get parent System object.
McSimulator< D > & simulator()
Get parent McSimulator object.
void incrementNAttempt()
Increment the number of attempted moves.
virtual void attemptMove()
Attempt unconstrained move.
virtual ~McMove()
Destructor.
virtual bool move()
Generate, attempt, and accept or reject a Monte Carlo move.
virtual void outputTimers(std::ostream &out)
Log output timing results.
void readProbability(std::istream &in)
Read the probability from file.
void incrementNFail()
Increment the number of failed moves.
virtual void setup()
Setup before the beginning of each simulation run.
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
Real periodic fields, SCFT and PS-FTS (CPU).
Definition param_pc.dox:2
PSCF package top-level namespace.
Definition param_pc.dox:1