PSCF v1.3
rpc/fts/montecarlo/ForceBiasMove.tpp
1#ifndef RPC_FORCE_BIAS_MOVE_TPP
2#define RPC_FORCE_BIAS_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 "ForceBiasMove.h"
12#include "McMove.h"
13#include <rpc/fts/montecarlo/McSimulator.h>
14#include <rpc/fts/compressor/Compressor.h>
15#include <rpc/system/System.h>
16#include <util/param/ParamComposite.h>
17#include <util/random/Random.h>
18
19
20namespace Pscf {
21namespace Rpc {
22
23 using namespace Util;
24
25 /*
26 * Constructor.
27 */
28 template <int D>
30 : McMove<D>(simulator),
31 w_(),
32 dwc_(),
33 mobility_(0.0)
34 { setClassName("ForceBiasMove"); }
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 ForceBiasMove<D>::readParameters(std::istream &in)
48 {
49 //Read the probability
51
52 // Read Brownian dynamics mobility parameter
53 read(in, "mobility", mobility_);
54
55 // Allocate memory for private containers
56 int nMonomer = system().mixture().nMonomer();
57 IntVec<D> meshDimensions = system().domain().mesh().dimensions();
58 w_.allocate(nMonomer);
59 for (int i=0; i < nMonomer; ++i) {
60 w_[i].allocate(meshDimensions);
61 }
62 dc_.allocate(nMonomer-1);
63 dwc_.allocate(nMonomer-1);
64 for (int i=0; i < nMonomer - 1; ++i) {
65 dc_[i].allocate(meshDimensions);
66 dwc_[i].allocate(meshDimensions);
67 }
68
69 }
70
71 template <int D>
73 {
74 // Check array capacities
75 int meshSize = system().domain().mesh().size();
76 int nMonomer = system().mixture().nMonomer();
77 UTIL_CHECK(w_.capacity() == nMonomer);
78 for (int i=0; i < nMonomer; ++i) {
79 UTIL_CHECK(w_[i].capacity() == meshSize);
80 }
81 UTIL_CHECK(dc_.capacity() == nMonomer-1);
82 UTIL_CHECK(dwc_.capacity() == nMonomer-1);
83 for (int i=0; i < nMonomer - 1; ++i) {
84 UTIL_CHECK(dc_[i].capacity() == meshSize);
85 UTIL_CHECK(dwc_[i].capacity() == meshSize);
86 }
87
89 }
90
91 /*
92 * Attempt unconstrained move
93 */
94 template <int D>
96 {
97 totalTimer_.start();
99
100 // Preconditions
101 UTIL_CHECK(simulator().hasWc());
102 UTIL_CHECK(simulator().hasDc());
103 UTIL_CHECK(simulator().hasHamiltonian());
104
105 // Array sizes and indices
106 const int nMonomer = system().mixture().nMonomer();
107 const int meshSize = system().domain().mesh().size();
108 int i, j, k;
109
110 // Get current Hamiltonian
111 double oldHamiltonian = simulator().hamiltonian();
112
113 // Save current state
114 simulator().saveState();
115
116 // Clear both eigen-components of the fields and hamiltonian
117 simulator().clearData();
118
119 attemptMoveTimer_.start();
120
121 // Copy current W fields from parent system into wc_
122 for (i = 0; i < nMonomer; ++i) {
123 w_[i] = system().w().rgrid(i);
124 }
125
126 // Copy current derivative fields from into member variable dc_
127 for (i = 0; i < nMonomer - 1; ++i) {
128 dc_[i] = simulator().dc(i);
129 }
130
131 // Constants for dynamics
132 const double vSystem = system().domain().unitCell().volume();
133 const double vNode = vSystem/double(meshSize);
134 const double a = -1.0*mobility_;
135 const double b = sqrt(2.0*mobility_/vNode);
136
137 // Modify local variables dwc_ and wc_
138 // Loop over eigenvectors of projected chi matrix
139 double dwd, dwr, evec;
140 for (j = 0; j < nMonomer - 1; ++j) {
141 RField<D> const & dc = dc_[j];
142 RField<D> & dwc = dwc_[j];
143 for (k = 0; k < meshSize; ++k) {
144 dwd = a*dc[k];
145 dwr = b*random().gaussian();
146 dwc[k] = dwd + dwr;
147 }
148
149 // Loop over monomer types
150 for (i = 0; i < nMonomer; ++i) {
151 RField<D> const & dwc = dwc_[j];
152 RField<D> & wn = w_[i];
153 evec = simulator().chiEvecs(j,i);
154 for (k = 0; k < meshSize; ++k) {
155 wn[k] += evec*dwc[k];
156 }
157 }
158
159 }
160
161 // Set modified fields in parent system
162 system().w().setRGrid(w_);
163 simulator().clearData();
164
165 attemptMoveTimer_.stop();
166
167 // Call compressor
168 compressorTimer_.start();
169 int compress = simulator().compressor().compress();
170 compressorTimer_.stop();
171
172 bool isConverged = false;
173 if (compress != 0){
175 simulator().restoreState();
176 } else {
177 isConverged = true;
178
179 // Compute eigenvector components of current fields
180 componentTimer_.start();
181 simulator().computeWc();
182 simulator().computeCc();
183 simulator().computeDc();
184 componentTimer_.stop();
185
186 // Evaluate new Hamiltonian
187 hamiltonianTimer_.start();
188 simulator().computeHamiltonian();
189 double newHamiltonian = simulator().hamiltonian();
190 double dH = newHamiltonian - oldHamiltonian;
191
192 // Compute force bias
193 double dp, dm;
194 double bias = 0.0;
195 for (j = 0; j < nMonomer - 1; ++j) {
196 RField<D> const & di = dc_[j];
197 RField<D> const & df = simulator().dc(j);
198 RField<D> const & dwc = dwc_[j];
199 for (k=0; k < meshSize; ++k) {
200 dp = 0.5*(di[k] + df[k]);
201 dm = 0.5*(di[k] - df[k]);
202 bias += dp*( dwc[k] + mobility_*dm );
203 }
204 }
205 bias *= vNode;
206 hamiltonianTimer_.stop();
207
208 // Accept or reject move
209 bool accept = false;
210 decisionTimer_.start();
211 double weight = exp(bias - dH);
212 accept = random().metropolis(weight);
213 if (accept) {
215 simulator().clearState();
216 } else {
217 simulator().restoreState();
218 }
219 decisionTimer_.stop();
220
221 }
222
223 totalTimer_.stop();
224 return isConverged;
225 }
226
227 /*
228 * Trivial default implementation - do nothing
229 */
230 template <int D>
233
234 template<int D>
235 void ForceBiasMove<D>::outputTimers(std::ostream& out)
236 {
237 out << "\n";
238 out << "ForceBiasMove time contributions:\n";
240 }
241
242}
243}
244#endif
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
Field of real double precision values on an FFT mesh.
Definition cpu/RField.h:29
ForceBiasMove(McSimulator< D > &simulator)
Constructor.
void readParameters(std::istream &in)
Read required parameters from file.
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
bool move()
Attempt and accept or reject force bias Monte-Carlo move.
void setClassName(const char *className)
Set class name string.
void outputTimers(std::ostream &out)
Return real move times contributions.
void setup()
Setup before the beginning of each simulation run.
void output()
Output statistics for this move (at the end of simulation)
void incrementNAccept()
Increment the number of accepted moves.
Timer attemptMoveTimer_
Timers for McMove.
Random & random()
Get Random number generator of parent System.
McMove(McSimulator< D > &simulator)
Constructor.
System< D > & system()
Get parent System object.
McSimulator< D > & simulator()
Get parent McSimulator object.
void incrementNAttempt()
Increment the number of attempted moves.
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.
Monte-Carlo simulation coordinator.
#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