PSCF v1.3.3
rpg/fts/montecarlo/RealMove.tpp
1#ifndef RPG_REAL_MOVE_TPP
2#define RPG_REAL_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 "RealMove.h"
12#include "McMove.h"
13#include <rpg/fts/montecarlo/McSimulator.h>
14#include <rpg/system/System.h>
15#include <rpg/solvers/Mixture.h>
16#include <rpg/field/Domain.h>
17#include <rpg/fts/VecOpFts.h>
18#include <prdc/cuda/VecOp.h>
19#include <pscf/math/IntVec.h>
20#include <pscf/cuda/CudaRandom.h>
21#include <util/param/ParamComposite.h>
22
23namespace Pscf {
24namespace Rpg {
25
26 using namespace Util;
27
28 /*
29 * Constructor.
30 */
31 template <int D>
33 : McMove<D>(simulator),
34 w_(),
35 dwc_(),
36 sigma_(0.0),
37 isAllocated_(false)
38 { setClassName("RealMove"); }
39
40 /*
41 * Destructor, empty default implementation.
42 */
43 template <int D>
46
47 /*
48 * ReadParameters, empty default implementation.
49 */
50 template <int D>
51 void RealMove<D>::readParameters(std::istream &in)
52 {
53 // Read the probability
55
56 // The standard deviation of the Gaussian distribution
57 read(in, "sigma", sigma_);
58
59 }
60
61
62 template <int D>
64 {
66 const int nMonomer = system().mixture().nMonomer();
67 IntVec<D> meshDimensions = system().domain().mesh().dimensions();
68
69 if (!isAllocated_){
70 w_.allocate(nMonomer);
71 for (int i=0; i < nMonomer; ++i) {
72 w_[i].allocate(meshDimensions);
73 }
74 dwc_.allocate(meshDimensions);
75 gaussianField_.allocate(meshDimensions);
76 isAllocated_ = true;
77 }
78 }
79
80 /*
81 * Attempt unconstrained move
82 */
83 template <int D>
85 {
86 const int nMonomer = system().mixture().nMonomer();
87 const int meshSize = system().domain().mesh().size();
88 int i, j;
89
90 // Copy current W fields from parent system into w_
91 for (i = 0; i < nMonomer; ++i) {
92 VecOp::eqV(w_[i], system().w().rgrid(i));
93 }
94
95 double evec;
96 double mean = 0.0;
97
98 // Loop over composition eigenvectors of projected chi matrix
99 for (j = 0; j < nMonomer - 1; ++j) {
100
101 // Generate Gaussian distributed random numbers
102 cudaRandom().normal(gaussianField_, sigma_, mean);
103 VecOp::eqV(dwc_, gaussianField_);
104
105 // Loop over monomer types
106 for (i = 0; i < nMonomer; ++i) {
107 RField<D> & w = w_[i];
108 evec = simulator().chiEvecs(j,i);
109 VecOp::addEqVc(w, dwc_, evec);
110 }
111
112 }
113
114 // Set w-fields of parent System
115 system().w().setRGrid(w_);
116
117 }
118
119
120 /*
121 * Trivial default implementation - do nothing
122 */
123 template <int D>
125 {}
126
127 template<int D>
128 void RealMove<D>::outputTimers(std::ostream& out)
129 {
130 out << "\n";
131 out << "RealMove time contributions:\n";
133 }
134
135}
136}
137#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
CudaRandom & cudaRandom()
Get cuda random number generator by reference.
void readProbability(std::istream &in)
Read the probability from file.
McSimulator< D > & simulator()
Get parent McSimulator object.
System< D > & system()
Get parent System object.
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.
Monte-Carlo simulation coordinator.
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
void readParameters(std::istream &in)
Read required parameters from file.
void setClassName(const char *className)
Set class name string.
void attemptMove()
Attempt unconstrained move.
RealMove(McSimulator< D > &simulator)
Constructor.
void setup()
Setup before the beginning of each simulation run.
void outputTimers(std::ostream &out)
Return real move times contributions.
void output()
Output statistics for this move (at the end of simulation)
void eqV(Array< double > &a, Array< double > const &b)
Vector assignment, a[i] = b[i].
Definition VecOp.cpp:23
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.
Definition VecOpMisc.cu:343
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.