1#ifndef RPC_FORCE_BIAS_MOVE_TPP
2#define RPC_FORCE_BIAS_MOVE_TPP
11#include "ForceBiasMove.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>
53 read(in,
"mobility", mobility_);
56 int nMonomer =
system().mixture().nMonomer();
58 w_.allocate(nMonomer);
59 for (
int i=0; i < nMonomer; ++i) {
60 w_[i].allocate(meshDimensions);
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);
75 int meshSize =
system().domain().mesh().size();
76 int nMonomer =
system().mixture().nMonomer();
78 for (
int i=0; i < nMonomer; ++i) {
83 for (
int i=0; i < nMonomer - 1; ++i) {
106 const int nMonomer =
system().mixture().nMonomer();
107 const int meshSize =
system().domain().mesh().size();
111 double oldHamiltonian =
simulator().hamiltonian();
122 for (i = 0; i < nMonomer; ++i) {
123 w_[i] =
system().w().rgrid(i);
127 for (i = 0; i < nMonomer - 1; ++i) {
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);
139 double dwd, dwr, evec;
140 for (j = 0; j < nMonomer - 1; ++j) {
143 for (k = 0; k < meshSize; ++k) {
145 dwr = b*
random().gaussian();
150 for (i = 0; i < nMonomer; ++i) {
154 for (k = 0; k < meshSize; ++k) {
155 wn[k] += evec*dwc[k];
162 system().w().setRGrid(w_);
168 compressorTimer_.start();
169 int compress =
simulator().compressor().compress();
170 compressorTimer_.stop();
172 bool isConverged =
false;
180 componentTimer_.start();
184 componentTimer_.stop();
187 hamiltonianTimer_.start();
189 double newHamiltonian =
simulator().hamiltonian();
190 double dH = newHamiltonian - oldHamiltonian;
195 for (j = 0; j < nMonomer - 1; ++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 );
206 hamiltonianTimer_.stop();
210 decisionTimer_.start();
211 double weight = exp(bias - dH);
212 accept =
random().metropolis(weight);
219 decisionTimer_.stop();
238 out <<
"ForceBiasMove time contributions:\n";
An IntVec<D, T> is a D-component vector of elements of integer type T.
Field of real double precision values on an FFT mesh.
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.
~ForceBiasMove()
Destructor.
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.
Real periodic fields, SCFT and PS-FTS (CPU).
PSCF package top-level namespace.