PSCF v1.2
rpc/fts/compressor/intra/IntraCorrelation.tpp
1#ifndef RPC_INTRACORRELATION_TPP
2#define RPC_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 <rpc/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 <prdc/crystal/shiftToMinimum.h>
18
19namespace Pscf {
20namespace Rpc{
21
22 using namespace Util;
23
24 // Constructor
25 template <int D>
27 : systemPtr_(&system),
28 kSize_(1)
29 { setClassName("IntraCorrelation"); }
30
31 // Destructor
32 template <int D>
35
36 template<int D>
37 void
39 {
40 // Local copies of system properties
41 Mixture<D> const & mixture = system().mixture();
42 UnitCell<D> const & unitCell = system().domain().unitCell();
43 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
44 const int nPolymer = mixture.nPolymer();
45 const double vMonomer = mixture.vMonomer();
46
47 // Compute Fourier space kMeshDimensions_
48 for (int i = 0; i < D; ++i) {
49 if (i < D - 1) {
50 kMeshDimensions_[i] = dimensions[i];
51 kSize_ *= dimensions[i];
52 } else {
53 kMeshDimensions_[i] = dimensions[i]/2 + 1;
54 kSize_ *= (dimensions[i]/2 + 1);
55 }
56 }
57 UTIL_CHECK(correlations.capacity() == kSize_);
58
59 // Allocate Gsq (k-space array of square wavenumber values)
61 Gsq.allocate(kSize_);
62
63 // Compute Gsq
64 IntVec<D> G, Gmin;
65 MeshIterator<D> iter;
66 iter.setDimensions(kMeshDimensions_);
67 for (iter.begin(); !iter.atEnd(); ++iter) {
68 G = iter.position();
69 Gmin = shiftToMinimum(G, dimensions, unitCell);
70 Gsq[iter.rank()] = unitCell.ksq(Gmin);
71 }
72
73 double phi, cPolymer, polymerLength;
74 double length, lengthA, lengthB, ksq;
75 double kuhn, kuhnA, kuhnB, eA, eB;
76 int monomerId, monomerIdA, monomerIdB, rank;
77
78 // Initialize correlations to zero
79 for (int i = 0; i < correlations.capacity(); ++i){
80 correlations[i] = 0.0;
81 }
82
83 // Loop over polymer species
84 for (int i = 0; i < nPolymer; i++){
85
86 // Local copies of polymer properties
87 Polymer<D> const & polymer = mixture.polymer(i);
88 const int nBlock = polymer.nBlock();
89 phi = polymer.phi();
90
91 // Compute polymerLength (sum of lengths of all blocks)
92 polymerLength = 0.0;
93 for (int j = 0; j < nBlock; j++) {
94 length = polymer.block(j).length();
95 polymerLength += length;
96 }
97
98 // Compute cPolymer (polymer number concentration)
99 cPolymer = phi/(polymerLength*vMonomer);
100
101 // Compute diagonal (j = k) contributions
102 for (int j = 0; j < nBlock; j++) {
103
104 monomerId = polymer.block(j).monomerId();
105 kuhn = mixture.monomer(monomerId).kuhn();
106 length = polymer.block(j).length();
107
108 // Loop over ksq to increment correlations
109 for (iter.begin(); !iter.atEnd(); ++iter) {
110 rank = iter.rank();
111 ksq = Gsq[rank];
112 correlations[rank] += cPolymer * Debye::d(ksq, length, kuhn);
113 }
114
115 }
116
117 // Compute off-diagonal contributions
118 if (nBlock > 1) {
119 EdgeIterator EdgeItr(polymer);
120
121 // Outer loop over blocks
122 for (int ia = 1; ia < nBlock; ++ia) {
123
124 // Block A properties
125 Block<D> const & blockA = polymer.block(ia);
126 lengthA = blockA.length();
127 monomerIdA = blockA.monomerId();
128 kuhnA = mixture.monomer(monomerIdA).kuhn();
129
130 // Inner loop over blocks
131 for (int ib = 0; ib < ia; ++ib) {
132
133 // Block B properties
134 Block<D> const & blockB = polymer.block(ib);
135 lengthB = blockB.length();
136 monomerIdB = blockB.monomerId();
137 kuhnB = mixture.monomer(monomerIdB).kuhn();
138
139 // Mean-square end-to-end length of segment between blocks
140 double d = 0.0;
141 int edgeId;
142 EdgeItr.begin(ia, ib);
143 while (EdgeItr.notEnd()) {
144 edgeId = EdgeItr.currentEdgeId();
145 if (edgeId != ia && edgeId != ib){
146 Block<D> const & blockK = polymer.block(edgeId);
147 double lengthK = blockK.length();
148 double monomerIdK = blockK.monomerId();
149 double kuhnK = mixture.monomer(monomerIdK).kuhn();
150 d += lengthK * kuhnK * kuhnK;
151 }
152 ++EdgeItr;
153 }
154
155 // Loop over ksq to increment correlations
156 double x;
157 double prefactor = 2.0*cPolymer;
158 for (iter.begin(); !iter.atEnd(); ++iter) {
159 rank = iter.rank();
160 ksq = Gsq[rank];
161 x = std::exp( -d * ksq / 6.0);
162 eA = Debye::e(ksq, lengthA, kuhnA);
163 eB = Debye::e(ksq, lengthB, kuhnB);
164 correlations[rank] += prefactor * x * eA * eB;
165 }
166
167 } // loop: block ib
168 } // loop: block ia
169 } // if (nBlock > 1)
170
171 } // loop over polymers
172
173 }
174
175}
176}
177#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.
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 modify intramolecular correlations.
Solver for a mixture of polymers and solvents.
Descriptor and solver for one polymer species.
Main class for SCFT or PS-FTS simulation of one system.
Definition rpc/System.h:100
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.