1#ifndef RPG_FORCE_BIAS_MOVE_TPP
2#define RPG_FORCE_BIAS_MOVE_TPP
11#include "ForceBiasMove.h"
13#include <rpg/fts/montecarlo/McSimulator.h>
14#include <rpg/fts/compressor/Compressor.h>
15#include <rpg/fts/VecOpFts.h>
16#include <rpg/system/System.h>
17#include <prdc/cuda/resources.h>
18#include <pscf/cuda/CudaRandom.h>
19#include <util/param/ParamComposite.h>
20#include <util/random/Random.h>
55 read(in,
"mobility", mobility_);
58 int nMonomer =
system().mixture().nMonomer();
60 w_.allocate(nMonomer);
61 for (
int i=0; i < nMonomer; ++i) {
62 w_[i].allocate(meshDimensions);
64 dc_.allocate(nMonomer-1);
65 dwc_.allocate(nMonomer-1);
66 for (
int i=0; i < nMonomer - 1; ++i) {
67 dc_[i].allocate(meshDimensions);
68 dwc_[i].allocate(meshDimensions);
77 int meshSize =
system().domain().mesh().size();
79 int nMonomer =
system().mixture().nMonomer();
81 for (
int i=0; i < nMonomer; ++i) {
86 for (
int i=0; i < nMonomer - 1; ++i) {
92 biasField_.allocate(dimensions);
93 gaussianField_.allocate(dimensions);
111 const int nMonomer =
system().mixture().nMonomer();
112 const int meshSize =
system().domain().mesh().size();
117 double oldHamiltonian =
simulator().hamiltonian();
128 for (i = 0; i < nMonomer; ++i) {
133 for (i = 0; i < nMonomer - 1; ++i) {
138 const double vSystem =
system().domain().unitCell().volume();
139 const double vNode = vSystem/double(meshSize);
140 double a = -1.0*mobility_;
141 double b = sqrt(2.0*mobility_/vNode);
150 for (j = 0; j < nMonomer - 1; ++j) {
154 cudaRandom().normal(gaussianField_, stddev, mean);
160 for (i = 0; i < nMonomer; ++i) {
170 system().w().setRGrid(w_);
176 compressorTimer_.start();
177 int compress =
simulator().compressor().compress();
179 compressorTimer_.stop();
181 bool isConverged =
false;
189 componentTimer_.start();
194 componentTimer_.stop();
197 hamiltonianTimer_.start();
199 double newHamiltonian =
simulator().hamiltonian();
200 double dH = newHamiltonian - oldHamiltonian;
204 for (j = 0; j < nMonomer - 1; ++j) {
210 bias += Reduce::sum(biasField_);
213 hamiltonianTimer_.stop();
216 decisionTimer_.start();
218 double weight = exp(bias - dH);
219 accept =
random().metropolis(weight);
226 decisionTimer_.stop();
245 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.
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
void setClassName(const char *className)
Set class name string.
void setup()
Setup before the beginning of each simulation run.
void output()
Output statistics for this move (at the end of simulation)
ForceBiasMove(McSimulator< D > &simulator)
Constructor.
void readParameters(std::istream &in)
Read required parameters from file.
bool move()
Attempt and accept or reject force bias Monte-Carlo move.
void outputTimers(std::ostream &out)
Return real move times contributions.
~ForceBiasMove()
Destructor.
CudaRandom & cudaRandom()
Get cuda random number generator by reference.
void incrementNAttempt()
Increment the number of attempted moves.
Random & random()
Get Random number generator of parent System.
void readProbability(std::istream &in)
Read the probability from file.
McSimulator< D > & simulator()
Get parent McSimulator object.
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.
Monte-Carlo simulation coordinator.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
void eqV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, const int beginIdA, const int beginIdB, const int n)
Vector assignment, a[i] = b[i], kernel wrapper (cudaReal).
void addVcVc(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c, DeviceArray< cudaReal > const &d, cudaReal const e)
Vector addition w/ coefficient, a[i] = (b[i]*c) + (d[i]*e), kernel wrapper.
void addEqVc(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c)
Vector addition in-place w/ coefficient, a[i] += b[i] * c, kernel wrapper.
void computeForceBias(DeviceArray< cudaReal > &result, DeviceArray< cudaReal > const &di, DeviceArray< cudaReal > const &df, DeviceArray< cudaReal > const &dwc, cudaReal mobility)
Compute force bias.
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.