1#ifndef RPC_BINARY_STRUCTURE_FACTOR_GRID_TPP
2#define RPC_BINARY_STRUCTURE_FACTOR_GRID_TPP
4#include "BinaryStructureFactorGrid.h"
6#include <rpc/fts/simulator/Simulator.h>
9#include <prdc/cpu/RField.h>
10#include <prdc/crystal/shiftToMinimum.h>
12#include <pscf/mesh/MeshIterator.h>
13#include <pscf/math/IntVec.h>
14#include <pscf/math/RealVec.h>
16#include <util/param/ParamComposite.h>
17#include <util/misc/FileMaster.h>
18#include <util/misc/ioUtil.h>
19#include <util/format/Int.h>
20#include <util/format/Dbl.h>
46 simulatorPtr_(&simulator),
47 systemPtr_(&(simulator.system())),
48 isInitialized_(false),
62 UTIL_CHECK(system().mixture().nMonomer() == 2);
65 readOutputFileName(in);
66 readOptional(in,
"nSamplePerBlock", nSamplePerBlock_);
76 const int nMonomer = system().mixture().nMonomer();
77 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
81 for (
int i = 0; i < D; ++i) {
83 kMeshDimensions_[i] = dimensions[i];
85 kMeshDimensions_[i] = dimensions[i]/2 + 1;
87 kSize_ *= kMeshDimensions_[i];
91 wKGrid_.allocate(nMonomer);
92 wm_.allocate(dimensions);
93 wk_.allocate(dimensions);
94 for (
int i = 0; i < nMonomer; ++i) {
95 wKGrid_[i].allocate(dimensions);
99 qList_.resize(kSize_);
102 UnitCell<D> const & unitCell = system().domain().unitCell();
108 Gmin = shiftToMinimum(G, dimensions, unitCell);
109 Gsq = unitCell.
ksq(Gmin);
110 qList_[itr.
rank()] = sqrt(Gsq);
116 structureFactors_.allocate(nWave_);
117 accumulators_.allocate(nWave_);
118 isInitialized_ =
true;
121 for (
int i = 0; i < nWave_; ++i) {
122 structureFactors_[i] = 0.0;
123 accumulators_[i].setNSamplePerBlock(nSamplePerBlock_);
124 accumulators_[i].clear();
126 if (!isInitialized_) {
127 UTIL_THROW(
"Error: object is not initialized");
139 if (isAtInterval(iStep)) {
142 for (
int i = 0; i < wm_.capacity(); ++i) {
143 wm_[i] = (system().w().rgrid(0)[i] - system().w().rgrid(1)[i])/2;
147 system().domain().fft().forwardTransform(wm_, wk_);
150 for (
int k=0; k< wk_.capacity(); k++) {
151 std::complex<double> wmKGrid(wk_[k][0], wk_[k][1]);
152 double squared_magnitude = std::norm(wmKGrid);
153 accumulators_[k].sample(squared_magnitude);
162 const double vSystem = system().domain().unitCell().volume();
163 const double vMonomer = system().mixture().vMonomer();
164 double n = vSystem / vMonomer;
165 double chi= system().interaction().chi(0,1);
170 structureFactors_[itr.
rank()] = n / (chi * chi) * accumulators_[itr.
rank()].average() - 1.0/(2.0*chi);
173 structureFactors_[itr.
rank()] /= vMonomer;
182 const int nMonomer = system().mixture().nMonomer();
183 for (
int i = 0; i < nMonomer; ++i) {
184 system().domain().fft().forwardTransform(system().w().rgrid()[i],
188 std::map<double, double> SMap;
190 int qListcapacity = (int)qList_.capacity();
192 UTIL_CHECK(qListcapacity == structureFactors_.capacity());
193 for (
int i = 0; i < qListcapacity; ++i) {
195 s = structureFactors_[i];
201 for (
auto& i : SMap) {
203 double sum = i.second;
204 int count = std::count(qList_.begin(), qList_.end(), q);
205 i.second = sum / count;
209 double tolerance = 1e-5;
210 auto it = SMap.begin();
211 double currentq = it->first;
212 double sumS = it->second;
214 for (++it; it != SMap.end(); ++it) {
215 double key = it->first;
216 double s = it->second;
217 if (std::abs(key - currentq) <= tolerance) {
221 averageSMap_[currentq] = sumS / count;
229 averageSMap_[currentq] = sumS / count;
238 computeStructureFactor();
239 averageStructureFactor();
242 system().fileMaster().openOutputFile(outputFileName(), outputFile_);
243 outputFile_ <<
"\t" <<
"q";
244 outputFile_ <<
"\t" <<
"\t" <<
"S(q)/(\u03C1 N)";
245 outputFile_<< std::endl;
246 for (
const auto& i : averageSMap_) {
248 double averageS = i.second;
249 outputFile_ <<
Dbl(q, 18, 8);
250 outputFile_ <<
Dbl(averageS, 18, 8);
251 outputFile_ << std::endl;
An IntVec<D, T> is a D-component vector of elements of integer type T.
Iterator over points in a Mesh<D>.
int rank() const
Get the rank of current element.
void begin()
Set iterator to the first point in the mesh.
bool atEnd() const
Is this the end (i.e., one past the last point)?
void setDimensions(const IntVec< D > &dimensions)
Set the grid dimensions in all directions.
IntVec< D > position() const
Get current position in the grid, as integer vector.
virtual double ksq(IntVec< D > const &k) const
Compute square magnitude of reciprocal lattice vector.
Base template for UnitCell<D> classes, D=1, 2 or 3.
Abstract base for periodic output and/or analysis actions.
void setup()
Clear accumulators.
void setClassName(const char *className)
Set class name string.
void readParameters(std::istream &in)
Read parameters from file.
BinaryStructureFactorGrid(Simulator< D > &simulator, System< D > &system)
Constructor.
void sample(long iStep)
Add particles to BinaryStructureFactor accumulators.
void averageStructureFactor()
Compute average S(k) over k of equal magnitude.
void computeStructureFactor()
Compute structure factor.
void output()
Output results to predefined output file.
Field theoretic simulator (base class).
Main class for SCFT or PS-FTS simulation of one system.
Wrapper for a double precision number, for formatted ostream output.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Periodic fields and crystallography.
PSCF package top-level namespace.
Utility classes for scientific computation.