PSCF v1.2
rpg/fts/compressor/intra/IntraCorrelation.tpp
1#ifndef RPG_INTRACORRELATION_TPP
2#define RPG_INTRACORRELATION_TPP
3
4/*
5* PSCF - Polymer Self-Consistent Field Theory
6*
7* Copyright 2016 - 2022, The Regents of the University of Minnesota
8* Distributed under the terms of the GNU General Public License.
9*/
10
11#include "IntraCorrelation.h"
12#include <rpg/System.h>
13#include <util/global.h>
14#include <pscf/mesh/MeshIterator.h>
15#include <pscf/chem/Debye.h>
16#include <pscf/chem/EdgeIterator.h>
17#include <pscf/iterator/NanException.h>
18#include <pscf/chem/PolymerType.h>
19#include <pscf/cuda/HostDArray.h>
20#include <prdc/crystal/shiftToMinimum.h>
21
22namespace Pscf {
23namespace Rpg{
24
25 using namespace Util;
26
27 // Constructor
28 template <int D>
30 : systemPtr_(&system),
31 kSize_(1)
32 { setClassName("IntraCorrelation"); }
33
34 // Destructor
35 template <int D>
38
39
40 template<int D>
41 void
43 {
44 // Local copies of system properties
45 Mixture<D> const & mixture = system().mixture();
46 UnitCell<D> const & unitCell = system().domain().unitCell();
47 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
48 const int nPolymer = mixture.nPolymer();
49 const double vMonomer = mixture.vMonomer();
50
51 // Compute Fourier space kMeshDimensions_
52 for (int i = 0; i < D; ++i) {
53 if (i < D - 1) {
54 kMeshDimensions_[i] = dimensions[i];
55 kSize_ *= dimensions[i];
56 } else {
57 kMeshDimensions_[i] = dimensions[i]/2 + 1;
58 kSize_ *= (dimensions[i]/2 + 1);
59 }
60 }
61 UTIL_CHECK(correlations.capacity() == kSize_);
62
63 // Allocate Gsq (k-space array of square wavenumber values)
65 Gsq.allocate(kSize_);
66
67 // Compute Gsq
68 IntVec<D> G, Gmin;
69 MeshIterator<D> iter;
70 iter.setDimensions(kMeshDimensions_);
71 for (iter.begin(); !iter.atEnd(); ++iter) {
72 G = iter.position();
73 Gmin = shiftToMinimum(G, dimensions, unitCell);
74 Gsq[iter.rank()] = unitCell.ksq(Gmin);
75 }
76
77 double phi, cPolymer, polymerLength;
78 double length, lengthA, lengthB, ksq;
79 double kuhn, kuhnA, kuhnB, eA, eB;
80 int monomerId, monomerIdA, monomerIdB, rank;
81
82 // Allocate local array
83 HostDArray<cudaReal> correlations_h;
84 correlations_h.allocate(kSize_);
85
86 // Initialize correlations_h to zero
87 for (int i = 0; i < correlations_h.capacity(); ++i){
88 correlations_h[i] = 0.0;
89 }
90
91 // Loop over polymer species
92 for (int i = 0; i < nPolymer; i++){
93
94 // Local copies of polymer properties
95 Polymer<D> const & polymer = mixture.polymer(i);
96 const int nBlock = polymer.nBlock();
97 phi = polymer.phi();
98
99 // Compute polymerLength (sum of lengths of all blocks)
100 polymerLength = 0.0;
101 for (int j = 0; j < nBlock; j++) {
102 length = polymer.block(j).length();
103 polymerLength += length;
104 }
105
106 // Compute cPolymer (polymer number concentration)
107 cPolymer = phi/(polymerLength*vMonomer);
108
109 // Compute diagonal (j = k) contributions
110 for (int j = 0; j < nBlock; j++) {
111
112 monomerId = polymer.block(j).monomerId();
113 kuhn = mixture.monomer(monomerId).kuhn();
114 length = polymer.block(j).length();
115
116 // Loop over ksq to increment correlations_h
117 for (iter.begin(); !iter.atEnd(); ++iter) {
118 rank = iter.rank();
119 ksq = Gsq[rank];
120 correlations_h[rank] += cPolymer * Debye::d(ksq, length, kuhn);
121 }
122
123 }
124
125 // Compute off-diagonal contributions
126 if (nBlock > 1) {
127 EdgeIterator EdgeItr(polymer);
128
129 // Outer loop over blocks
130 for (int ia = 1; ia < nBlock; ++ia) {
131
132 // Block A properties
133 Block<D> const & blockA = polymer.block(ia);
134 lengthA = blockA.length();
135 monomerIdA = blockA.monomerId();
136 kuhnA = mixture.monomer(monomerIdA).kuhn();
137
138 // Inner loop over blocks
139 for (int ib = 0; ib < ia; ++ib) {
140
141 // Block B properties
142 Block<D> const & blockB = polymer.block(ib);
143 lengthB = blockB.length();
144 monomerIdB = blockB.monomerId();
145 kuhnB = mixture.monomer(monomerIdB).kuhn();
146
147 // Mean-square end-to-end length of segment between blocks
148 double d = 0.0;
149 int edgeId;
150 EdgeItr.begin(ia, ib);
151 while (EdgeItr.notEnd()) {
152 edgeId = EdgeItr.currentEdgeId();
153 if (edgeId != ia && edgeId != ib){
154 Block<D> const & blockK = polymer.block(edgeId);
155 double lengthK = blockK.length();
156 double monomerIdK = blockK.monomerId();
157 double kuhnK = mixture.monomer(monomerIdK).kuhn();
158 d += lengthK * kuhnK * kuhnK;
159 }
160 ++EdgeItr;
161 }
162
163 // Loop over ksq to increment correlations_h
164 double x;
165 double prefactor = 2.0*cPolymer;
166 for (iter.begin(); !iter.atEnd(); ++iter) {
167 rank = iter.rank();
168 ksq = Gsq[rank];
169 x = std::exp( -d * ksq / 6.0);
170 eA = Debye::e(ksq, lengthA, kuhnA);
171 eB = Debye::e(ksq, lengthB, kuhnB);
172 correlations_h[rank] += prefactor * x * eA * eB;
173 }
174
175 } // loop: block ib
176 } // loop: block ia
177 } // if (nBlock > 1)
178
179 } // loop over polymers
180
181 // Copy local array to device (gpu) memory
182 correlations = correlations_h;
183
184 }
185
186}
187}
188#endif
Edge iterator for graph associated with a polymer.
int currentEdgeId() const
Get index of the current edge.
bool notEnd() const
Return true iff currentId != targetId.
void begin(int sourceId, int targetId)
Initialize iterator.
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.
IntVec< D > position() const
Get current position in the grid, as integer vector.
Polymer & polymer(int id)
Get a polymer object.
double vMonomer() const
Get monomer reference volume (set to 1.0 by default).
Monomer const & monomer(int id) const
Get a Monomer type descriptor (const reference).
int nPolymer() const
Get number of polymer species.
double kuhn() const
Statistical segment length (random walk step size).
Definition Monomer.h:131
Block & block(int id)
Get a specified Block (block solver and descriptor)
double phi() const
Get the overall volume fraction for this species.
Definition Species.h:90
int nBlock() const
Number of blocks.
Field of real double precision values on an FFT mesh.
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
Block within a branched polymer.
int monomerId() const
Get the monomer type id.
Definition Edge.h:235
double length() const
Get the length (number of monomers) in this block.
Definition Edge.h:253
void computeIntraCorrelations(RField< D > &correlations)
Compute and return intramolecular correlations.
Solver for a mixture of polymers and solvents.
Descriptor and solver for a branched polymer species.
Main class for calculations that represent one system.
Definition rpg/System.h:107
int capacity() const
Return allocated size.
Definition Array.h:159
Dynamically allocatable contiguous array template.
void allocate(int capacity)
Allocate the underlying C array.
Definition DArray.h:199
void setClassName(const char *className)
Set class name string.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
double d(double ksq, double length, double kuhn)
Compute and return a homopolymer correlation function.
Definition Debye.cpp:17
double e(double ksq, double length, double kuhn)
Compute and return one-sided correlation factor for one block.
Definition Debye.cpp:32
PSCF package top-level namespace.
Definition param_pc.dox:1
Utility classes for scientific computation.