1#ifndef RPG_INTRACORRELATION_TPP
2#define RPG_INTRACORRELATION_TPP
11#include "IntraCorrelation.h"
12#include <rpg/System.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>
30 : systemPtr_(&system),
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();
52 for (
int i = 0; i < D; ++i) {
54 kMeshDimensions_[i] = dimensions[i];
55 kSize_ *= dimensions[i];
57 kMeshDimensions_[i] = dimensions[i]/2 + 1;
58 kSize_ *= (dimensions[i]/2 + 1);
73 Gmin = shiftToMinimum(G, dimensions, unitCell);
74 Gsq[iter.
rank()] = unitCell.
ksq(Gmin);
77 double phi, cPolymer, polymerLength;
78 double length, lengthA, lengthB, ksq;
79 double kuhn, kuhnA, kuhnB, eA, eB;
80 int monomerId, monomerIdA, monomerIdB, rank;
87 for (
int i = 0; i < correlations_h.
capacity(); ++i){
88 correlations_h[i] = 0.0;
92 for (
int i = 0; i < nPolymer; i++){
96 const int nBlock = polymer.
nBlock();
101 for (
int j = 0; j < nBlock; j++) {
103 polymerLength += length;
107 cPolymer = phi/(polymerLength*vMonomer);
110 for (
int j = 0; j < nBlock; j++) {
120 correlations_h[rank] += cPolymer *
Debye::d(ksq, length, kuhn);
130 for (
int ia = 1; ia < nBlock; ++ia) {
134 lengthA = blockA.
length();
139 for (
int ib = 0; ib < ia; ++ib) {
143 lengthB = blockB.
length();
150 EdgeItr.
begin(ia, ib);
151 while (EdgeItr.
notEnd()) {
153 if (edgeId != ia && edgeId != ib){
155 double lengthK = blockK.
length();
158 d += lengthK * kuhnK * kuhnK;
165 double prefactor = 2.0*cPolymer;
169 x = std::exp( -d * ksq / 6.0);
172 correlations_h[rank] += prefactor * x * eA * eB;
182 correlations = correlations_h;
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.
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.
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).
Block & block(int id)
Get a specified Block (block solver and descriptor)
double phi() const
Get the overall volume fraction for this species.
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.
Block within a branched polymer.
int monomerId() const
Get the monomer type id.
double length() const
Get the length (number of monomers) in this block.
void computeIntraCorrelations(RField< D > &correlations)
Compute and return intramolecular correlations.
IntraCorrelation(System< D > &system)
Constructor.
~IntraCorrelation()
Destructor.
Solver for a mixture of polymers and solvents.
Descriptor and solver for a branched polymer species.
Main class for calculations that represent one system.
int capacity() const
Return allocated size.
Dynamically allocatable contiguous array template.
void allocate(int capacity)
Allocate the underlying C array.
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.
double d(double ksq, double length, double kuhn)
Compute and return a homopolymer correlation function.
double e(double ksq, double length, double kuhn)
Compute and return one-sided correlation factor for one block.
PSCF package top-level namespace.
Utility classes for scientific computation.