PSCF v1.3
rpg/fts/analyzer/BinaryStructureFactorGrid.tpp
1#ifndef RPG_BINARY_STRUCTURE_FACTOR_GRID_TPP
2#define RPG_BINARY_STRUCTURE_FACTOR_GRID_TPP
3
4#include "BinaryStructureFactorGrid.h"
5
6#include <rpg/fts/simulator/Simulator.h>
7#include <rpg/system/System.h>
8
9#include <prdc/crystal/shiftToMinimum.h>
10#include <prdc/cuda/FFT.h>
11#include <prdc/cuda/resources.h>
12#include <prdc/cuda/complex.h>
13
14#include <pscf/inter/Interaction.h>
15#include <pscf/mesh/MeshIterator.h>
16
17#include <util/misc/FileMaster.h>
18#include <util/misc/ioUtil.h>
19#include <util/format/Dbl.h>
20#include <util/format/Int.h>
21
22#include <fftw3.h>
23
24#include <vector>
25#include <unordered_map>
26#include <algorithm>
27
28namespace Pscf {
29namespace Rpg {
30
31 using namespace Util;
32 using namespace Pscf::Prdc;
33 using namespace Pscf::Prdc::Cuda;
34
35 /*
36 * Constructor.
37 */
38 template <int D>
40 : Analyzer<D>(),
43 isInitialized_(false),
44 nSamplePerBlock_(1)
45 { setClassName("BinaryStructureFactorGrid"); }
46
47 /*
48 * Read parameters from file, and allocate memory.
49 */
50 template <int D>
52 {
53 // Precondition: Require that the system has two monomer types
54 UTIL_CHECK(system().mixture().nMonomer() == 2);
55
56 readInterval(in);
58 readOptional(in,"nSamplePerBlock", nSamplePerBlock_);
59 }
60
61 /*
62 * BinaryStructureFactorGrid setup
63 */
64 template <int D>
66 {
67 // Copies of system properties
68 const int nMonomer = system().mixture().nMonomer();
69 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
70
71 FFT<D>::computeKMesh(dimensions, kMeshDimensions_, kSize_);
72
73 #if 0
74 // Compute Fourier space kMeshDimensions_
75 for (int i = 0; i < D; ++i) {
76 if (i < D - 1) {
77 kMeshDimensions_[i] = dimensions[i];
78 } else {
79 kMeshDimensions_[i] = dimensions[i]/2 + 1;
80 }
81 }
82 kSize_ = 1;
83 for(int i = 0; i < D; ++i) {
84 kSize_ *= kMeshDimensions_[i];
85 }
86 #endif
87
88 nWave_ = kSize_;
89
90 if (!isInitialized_){
91 wm_.allocate(dimensions);
92 wk_.allocate(dimensions);
93 }
94
95 structureFactors_.allocate(nWave_);
96 accumulators_.allocate(nWave_);
97
98 // Convert real grid to KGrid format
99 qList_.resize(kSize_);
100 RField<D> const & kSq = system().domain().waveList().kSq();
101 HostDArray<cudaReal> kSqHost(kSize_);
102 kSqHost = kSq;
103 MeshIterator<D> itr;
104 itr.setDimensions(kMeshDimensions_);
105 for (itr.begin(); !itr.atEnd(); ++itr){
106 qList_[itr.rank()] = sqrt(double(kSqHost[itr.rank()]));
107 }
108
109 // Clear accumulators
110 for (int i = 0; i < nWave_; ++i) {
111 structureFactors_[i] = 0.0;
112 accumulators_[i].setNSamplePerBlock(nSamplePerBlock_);
113 accumulators_[i].clear();
114 }
115
116 isInitialized_ = true;
117
118 }
119
120 /*
121 * Increment structure factors for all wavevectors and modes.
122 */
123 template <int D>
125 {
126 UTIL_CHECK(system().w().hasData());
127
128 if (isAtInterval(iStep)) {
129 //IntVec<D> const & dimensions
130 // = system().domain().mesh().dimensions();
131
132 // Compute W: (rgrid(0) - rgrid(1)) / 2
133 VecOp::addVcVc(wm_, system().w().rgrid(0), 0.5,
134 system().w().rgrid(1), -0.5);
135
136 // Convert real grid to KGrid format
137 system().domain().fft().forwardTransform(wm_, wk_);
138
139 HostDArray<cudaComplex> wkCpu(kSize_);
140 wkCpu = wk_; // copy from device to host
141 for (int k = 0; k < wk_.capacity(); k++) {
142 accumulators_[k].sample(absSq<cudaComplex, cudaReal>(wkCpu[k]));
143 }
144 }
145
146 }
147
148 template <int D>
150 {
151 const double vSystem = system().domain().unitCell().volume();
152 const double vMonomer = system().mixture().vMonomer();
153 double n = vSystem / vMonomer;
154 double chi= system().interaction().chi(0,1);
155 MeshIterator<D> itr;
156 itr.setDimensions(kMeshDimensions_);
157 for (itr.begin(); !itr.atEnd(); ++itr) {
158 // Compute vS(q)
159 structureFactors_[itr.rank()] = n / (chi * chi) * accumulators_[itr.rank()].average() - 1.0/(2.0*chi);
160
161 // Compute S(q)
162 structureFactors_[itr.rank()] /= vMonomer;
163 }
164 }
165
166 // Average S(k) over k of equal magnitude
167 template <int D>
169 {
170 UTIL_CHECK(qList_.capacity() == structureFactors_.capacity());
171
172 std::map<double, double> SMap;
173 for (int i = 0; i < qList_.capacity(); ++i) {
174 double q = qList_[i];
175 double s = structureFactors_[i];
176 SMap[q] += s;
177 }
178
179 // Average structure factor with same magnitude value of q
180 for (auto& i : SMap) {
181 double q = i.first;
182 double sum = i.second;
183 int count = std::count(qList_.begin(), qList_.end(), q);
184 i.second = sum / count;
185 }
186
187 // Average structure factor for q in range of +-tolerance
188 double tolerance = 1e-5;
189 auto it = SMap.begin();
190 double currentq = it->first;
191 double sumS = it->second;
192 int count = 1;
193 for (++it; it != SMap.end(); ++it) {
194 double key = it->first;
195 double s = it->second;
196 if (std::abs(key - currentq) <= tolerance) {
197 sumS += s;
198 count++;
199 } else {
200 averageSMap_[currentq] = sumS / count;
201 currentq = key;
202 sumS = s;
203 count = 1;
204 }
205 }
206
207 // Add last batch of data (last number within tolerance)
208 averageSMap_[currentq] = sumS / count;
209 }
210
211
212
213
214 /*
215 * Output final results to output file.
216 */
217 template <int D>
219 {
222
223 // Output structure factors to one file
224 system().fileMaster().openOutputFile(outputFileName(), outputFile_);
225 outputFile_ << "\t" << "q";
226 outputFile_ << "\t" <<"\t" <<"S(q)";
227 outputFile_<< std::endl;
228 for (const auto& i : averageSMap_) {
229 double q = i.first;
230 double averageS = i.second;
231 outputFile_ << Dbl(q, 18, 8);
232 outputFile_ << Dbl(averageS, 18, 8);
233 outputFile_ << std::endl;
234 }
235 outputFile_.close();
236 }
237
238}
239}
240#endif
Template for dynamic array stored in host CPU memory.
Definition HostDArray.h:43
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.
static void computeKMesh(IntVec< D > const &rMeshDimensions, IntVec< D > &kMeshDimensions, int &kSize)
Compute dimensions and size of k-space mesh for DFT of real data.
Definition cpu/FFT.tpp:262
Field of real double precision values on an FFT mesh.
Definition cpu/RField.h:29
const std::string & outputFileName() const
Return outputFileName string.
Analyzer()
Default constructor.
void readOutputFileName(std::istream &in)
Read outputFileName from file.
void readInterval(std::istream &in)
Read interval from file, with error checking.
bool isAtInterval(long counter) const
Return true iff counter is a multiple of the interval.
void sample(long iStep)
Add particles to BinaryStructureFactor accumulators.
void setClassName(const char *className)
Set class name string.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
void averageStructureFactor()
Compute average S(k) over k of equal magnitude.
System< D > * systemPtr_
Pointer to the parent system.
void output()
Output results to predefined output file.
Simulator< D > & simulator()
Return reference to parent Simulator.
void readParameters(std::istream &in)
Read parameters from file.
System< D > & system()
Return reference to parent system.
Simulator< D > * simulatorPtr_
Pointer to parent Simulator.
BinaryStructureFactorGrid(Simulator< D > &simulator, System< D > &system)
Constructor.
Field theoretic simulator (base class).
Main class, representing one complete system.
Wrapper for a double precision number, for formatted ostream output.
Definition Dbl.h:40
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
RT absSq(CT const &a)
Return square of absolute magnitude of a complex number.
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.
Definition VecOpMisc.cu:304
Fields, FFTs, and utilities for periodic boundary conditions (CUDA)
Definition Reduce.cpp:14
Periodic fields and crystallography.
Definition CField.cpp:11
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.
Definition param_pc.dox:1