PSCF v1.2
rpc/fts/analyzer/BinaryStructureFactorGrid.tpp
1#ifndef RPC_BINARY_STRUCTURE_FACTOR_GRID_TPP
2#define RPC_BINARY_STRUCTURE_FACTOR_GRID_TPP
3
4#include "BinaryStructureFactorGrid.h"
5
6#include <rpc/fts/simulator/Simulator.h>
7#include <rpc/System.h>
8
9#include <prdc/cpu/RField.h>
10#include <prdc/crystal/shiftToMinimum.h>
11
12#include <pscf/mesh/MeshIterator.h>
13#include <pscf/math/IntVec.h>
14#include <pscf/math/RealVec.h>
15
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>
21#include <util/global.h>
22
23#include <fftw3.h>
24
25#include <iostream>
26#include <complex>
27#include <vector>
28#include <algorithm>
29
30//#include <unordered_map>
31
32namespace Pscf {
33namespace Rpc {
34
35 using namespace Util;
36 using namespace Pscf::Prdc;
37
38 /*
39 * Constructor.
40 */
41 template <int D>
43 Simulator<D>& simulator,
44 System<D>& system)
45 : Analyzer<D>(),
46 simulatorPtr_(&simulator),
47 systemPtr_(&(simulator.system())),
48 isInitialized_(false),
49 nSamplePerBlock_(1),
50 kMeshDimensions_(0),
51 kSize_(0)
52 { setClassName("BinaryStructureFactorGrid"); }
53
54
55 /*
56 * Read parameters from file, and allocate memory.
57 */
58 template <int D>
60 {
61 // Precondition: Require that the system has two monomer types
62 UTIL_CHECK(system().mixture().nMonomer() == 2);
63
64 readInterval(in);
65 readOutputFileName(in);
66 readOptional(in,"nSamplePerBlock", nSamplePerBlock_);
67 }
68
69 /*
70 * BinaryStructureFactorGrid setup
71 */
72 template <int D>
74 {
75 //Allocate variables
76 const int nMonomer = system().mixture().nMonomer();
77 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
78
79 // Compute Fourier space kMeshDimensions_ and kSize_
80 kSize_ = 1;
81 for (int i = 0; i < D; ++i) {
82 if (i < D - 1) {
83 kMeshDimensions_[i] = dimensions[i];
84 } else {
85 kMeshDimensions_[i] = dimensions[i]/2 + 1;
86 }
87 kSize_ *= kMeshDimensions_[i];
88 }
89
90 if (!isInitialized_){
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);
96 }
97
98 // Compute qList
99 qList_.resize(kSize_);
100 IntVec<D> G;
101 IntVec<D> Gmin;
102 UnitCell<D> const & unitCell = system().domain().unitCell();
103 MeshIterator<D> itr(kMeshDimensions_);
104 double Gsq;
105 for (itr.begin(); !itr.atEnd(); ++itr){
106 // Obtain square magnitude of reciprocal basis vector
107 G = itr.position();
108 Gmin = shiftToMinimum(G, dimensions, unitCell);
109 Gsq = unitCell.ksq(Gmin);
110 qList_[itr.rank()] = sqrt(Gsq);
111 }
112
113 }
114
115 nWave_ = kSize_;
116 structureFactors_.allocate(nWave_);
117 accumulators_.allocate(nWave_);
118 isInitialized_ = true;
119
120 // Clear accumulators
121 for (int i = 0; i < nWave_; ++i) {
122 structureFactors_[i] = 0.0;
123 accumulators_[i].setNSamplePerBlock(nSamplePerBlock_);
124 accumulators_[i].clear();
125 }
126 if (!isInitialized_) {
127 UTIL_THROW("Error: object is not initialized");
128 }
129 }
130
131 /*
132 * Increment structure factors for all wavevectors and modes.
133 */
134 template <int D>
136 {
137 UTIL_CHECK(system().w().hasData());
138
139 if (isAtInterval(iStep)) {
140
141 // Compute W-
142 for (int i = 0; i < wm_.capacity(); ++i) {
143 wm_[i] = (system().w().rgrid(0)[i] - system().w().rgrid(1)[i])/2;
144 }
145
146 // Convert real grid to KGrid format
147 system().domain().fft().forwardTransform(wm_, wk_);
148
149 // Pass square magnitudes of Fourier components to accumulators
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);
154 }
155
156 }
157 }
158
159 template <int D>
161 {
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);
166 MeshIterator<D> itr;
167 itr.setDimensions(wKGrid_[0].dftDimensions());
168 for (itr.begin(); !itr.atEnd(); ++itr) {
169 // Compute vS(q)
170 structureFactors_[itr.rank()] = n / (chi * chi) * accumulators_[itr.rank()].average() - 1.0/(2.0*chi);
171
172 // Compute S(q)
173 structureFactors_[itr.rank()] /= vMonomer;
174 }
175 }
176
177 // Average S(k) over k of equal magnitude
178 template <int D>
180 {
181 // Convert real grid to KGrid format
182 const int nMonomer = system().mixture().nMonomer();
183 for (int i = 0; i < nMonomer; ++i) {
184 system().domain().fft().forwardTransform(system().w().rgrid()[i],
185 wKGrid_[i]);
186 }
187
188 std::map<double, double> SMap;
189 {
190 int qListcapacity = (int)qList_.capacity();
191 double q, s;
192 UTIL_CHECK(qListcapacity == structureFactors_.capacity());
193 for (int i = 0; i < qListcapacity; ++i) {
194 q = qList_[i];
195 s = structureFactors_[i];
196 SMap[q] += s;
197 }
198 }
199
200 // Average structure factor with same magnitude value of q
201 for (auto& i : SMap) {
202 double q = i.first;
203 double sum = i.second;
204 int count = std::count(qList_.begin(), qList_.end(), q);
205 i.second = sum / count;
206 }
207
208 // Average structure factor for q in range of +-tolerance
209 double tolerance = 1e-5;
210 auto it = SMap.begin();
211 double currentq = it->first;
212 double sumS = it->second;
213 int count = 1;
214 for (++it; it != SMap.end(); ++it) {
215 double key = it->first;
216 double s = it->second;
217 if (std::abs(key - currentq) <= tolerance) {
218 sumS += s;
219 count++;
220 } else {
221 averageSMap_[currentq] = sumS / count;
222 currentq = key;
223 sumS = s;
224 count = 1;
225 }
226 }
227
228 // Add last batch of data (last number within tolerance)
229 averageSMap_[currentq] = sumS / count;
230 }
231
232 /*
233 * Output final results to output file.
234 */
235 template <int D>
237 {
238 computeStructureFactor();
239 averageStructureFactor();
240
241 // Output structure factors to one file
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_) {
247 double q = i.first;
248 double averageS = i.second;
249 outputFile_ << Dbl(q, 18, 8);
250 outputFile_ << Dbl(averageS, 18, 8);
251 outputFile_ << std::endl;
252 }
253 outputFile_.close();
254 }
255
256}
257}
258#endif
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
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.
Definition rpg/System.h:34
Abstract base for periodic output and/or analysis actions.
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 output()
Output results to predefined output file.
Field theoretic simulator (base class).
Definition rpc/System.h:38
Main class for SCFT or PS-FTS simulation of one system.
Definition rpc/System.h:100
Wrapper for a double precision number, for formatted ostream output.
Definition Dbl.h:40
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Definition global.h:51
Periodic fields and crystallography.
Definition CField.cpp:11
PSCF package top-level namespace.
Definition param_pc.dox:1
Utility classes for scientific computation.