13#include <prdc/cpu/WaveList.h>
14#include <prdc/cpu/FFT.h>
15#include <prdc/crystal/UnitCell.h>
16#include <prdc/crystal/shiftToMinimum.h>
18#include <pscf/chem/Edge.h>
19#include <pscf/mesh/Mesh.h>
20#include <pscf/mesh/MeshIterator.h>
21#include <pscf/math/IntVec.h>
23#include <util/containers/DMatrix.h>
24#include <util/containers/DArray.h>
25#include <util/containers/FSArray.h>
42 unitCellPtr_(nullptr),
43 waveListPtr_(nullptr),
79 waveListPtr_ = &wavelist;
96 UTIL_CHECK(mesh().dimensions() == fft().meshDimensions());
103 expKsq_.allocate(kMeshDimensions_);
104 expKsq2_.allocate(kMeshDimensions_);
105 expW_.allocate(mesh().dimensions());
106 qr_.allocate(mesh().dimensions());
107 qk_.allocate(mesh().dimensions());
108 qr2_.allocate(mesh().dimensions());
109 qk2_.allocate(mesh().dimensions());
111 expW2_.allocate(mesh().dimensions());
113 expWInv_.allocate(mesh().dimensions());
117 cField().allocate(mesh().dimensions());
126 tempNs = floor(
length() / (2.0 *
ds) + 0.5);
131 ds_ =
length()/double(ns_-1);
162 tempNs = floor(
length()/(2.0 *dsTarget_) + 0.5 );
167 ds_ =
length()/double(ns_-1);
203 void Block<D>::computeExpKsq()
209 double bSqFactor = -1.0*kuhn()*kuhn() / 6.0;
215 if (!waveListPtr_->hasKSq()) {
216 waveListPtr_->computeKSq();
218 RField<D> const & kSq = waveListPtr_->kSq();
226 arg = kSq[i]*bSqFactor;
227 expKsq_[i] = exp(arg);
228 expKsq2_[i] = exp(0.5*arg);
242 int nx = mesh().size();
250 for (
int i = 0; i < nx; ++i) {
253 expW2_[i] = exp(0.5*arg);
257 for (
int i = 0; i < nx; ++i) {
260 expWInv_[i] = 1.0/expW_[i];
280 int nx = mesh().size();
283 UTIL_CHECK(mesh().dimensions() == fft().meshDimensions());
287 int nk = qk_.capacity();
306 for (i = 0; i < nx; ++i) {
307 qr_[i] = q[i]*expW_[i];
308 qr2_[i] = q[i]*expW2_[i];
310 fft().forwardTransform(qr_, qk_);
311 fft().forwardTransform(qr2_, qk2_);
312 for (i = 0; i < nk; ++i) {
313 qk_[i][0] *= expKsq_[i];
314 qk_[i][1] *= expKsq_[i];
315 qk2_[i][0] *= expKsq2_[i];
316 qk2_[i][1] *= expKsq2_[i];
318 fft().inverseTransformUnsafe(qk_, qr_);
319 fft().inverseTransformUnsafe(qk2_, qr2_);
320 for (i = 0; i < nx; ++i) {
321 qr_[i] = qr_[i]*expW_[i];
322 qr2_[i] = qr2_[i]*expW_[i];
330 fft().forwardTransform(qr2_, qk2_);
331 for (i = 0; i < nk; ++i) {
332 qk2_[i][0] *= expKsq2_[i];
333 qk2_[i][1] *= expKsq2_[i];
335 fft().inverseTransformUnsafe(qk2_, qr2_);
336 for (i = 0; i < nx; ++i) {
337 qr2_[i] = qr2_[i]*expW2_[i];
341 for (i = 0; i < nx; ++i) {
342 qout[i] = (4.0*qr2_[i] - qr_[i])/3.0;
367 int nx = mesh().size();
371 int nk = qk_.capacity();
376 fft().forwardTransform(q, qk_);
377 for (
int i = 0; i < nk; ++i) {
378 qk_[i][0] *= expKsq_[i];
379 qk_[i][1] *= expKsq_[i];
381 fft().inverseTransformUnsafe(qk_, qout);
393 int nx = mesh().size();
397 int nk = qk_.capacity();
402 fft().forwardTransform(q, qk_);
403 for (
int i = 0; i < nk; ++i) {
404 qk_[i][0] *= expKsq2_[i];
405 qk_[i][1] *= expKsq2_[i];
407 fft().inverseTransformUnsafe(qk_, qout);
417 int nx = mesh().size();
423 for (
int i = 0; i < nx; ++i) {
436 int nx = mesh().size();
445 for (i = 0; i < nx; ++i) {
456 for (i = 0; i < nx; ++i) {
457 cField()[i] += p0.
q(0)[i]*p1.
q(ns_ - 1)[i];
458 cField()[i] += p0.
q(ns_ -1)[i]*p1.
q(0)[i];
463 for (j = 1; j < (ns_ -1); j += 2) {
466 for (i = 0; i < nx; ++i) {
468 cField()[i] += qf[i] * qr[i] * 4.0;
473 for (j = 2; j < (ns_ -2); j += 2) {
476 for (i = 0; i < nx; ++i) {
478 cField()[i] += qf[i] * qr[i] * 2.0;
483 prefactor *= ds_ / 3.0;
484 for (i = 0; i < nx; ++i) {
498 int nx = mesh().size();
507 for (i = 0; i < nx; ++i) {
517 for (j = 1; j < (ns_ -1); ++j) {
520 for (i = 0; i < nx; ++i) {
521 cField()[i] += qf[i] * qr[i] * expWInv_[i];
526 for (i = 0; i < nx; ++i) {
539 int nx = mesh().size();
542 UTIL_CHECK(mesh().dimensions() == fft().meshDimensions());
554 if (!waveListPtr_->hasdKSq()) {
555 waveListPtr_->computedKSq();
560 const int nParam = unitCell().nParameter();
561 for (
int i = 0; i < nParam; ++i) {
566 const double bSq =
kuhn()*
kuhn()/6.0;
567 double dels, prod, increment;
571 for (
int j = 0; j < ns_ ; ++j) {
574 fft().forwardTransform(qr_, qk_);
576 qr2_ = p1.
q(ns_ - 1 - j);
577 fft().forwardTransform(qr2_, qk2_);
581 if (j != 0 && j != ns_ - 1) {
590 for (n = 0; n < nParam ; ++n) {
594 for (m = 0; m < kSize_ ; ++m) {
595 prod = (qk2_[m][0] * qk_[m][0]) + (qk2_[m][1] * qk_[m][1]);
599 increment *= bSq * dels;
600 dQ[n] = dQ[n] - increment;
607 for (
int n = 0; n < nParam; ++n) {
608 stress_.append(-1.0*prefactor*dQ[n]);
621 int nx = mesh().size();
624 UTIL_CHECK(mesh().dimensions() == fft().meshDimensions());
636 if (!waveListPtr_->hasdKSq()) {
637 waveListPtr_->computedKSq();
642 int nParam = unitCell().nParameter();
643 for (
int i = 0; i < nParam; ++i) {
647 const double bSq =
kuhn()*
kuhn()/6.0;
648 double increment, prod;
655 fft().forwardTransform(qr_, qk_);
658 qr2_ = p1.
q(ns_ - 2);
659 fft().forwardTransform(qr2_, qk2_);
662 for (
int n = 0; n < nParam ; ++n) {
666 for (
int m = 0; m < kSize_ ; ++m) {
667 prod = (qk2_[m][0] * qk_[m][0]) + (qk2_[m][1] * qk_[m][1]);
668 prod *= dKSq[m]*expKsq2_[m];
671 increment *= 0.5*bSq;
672 dQ[n] = dQ[n] - increment;
677 for (
int j = 1; j < ns_ - 2 ; ++j) {
681 fft().forwardTransform(qr_, qk_);
684 qr2_ = p1.
q(ns_ - 2 - j);
685 fft().forwardTransform(qr2_, qk2_);
688 for (
int n = 0; n < nParam ; ++n) {
692 for (
int m = 0; m < kSize_ ; ++m) {
693 prod = (qk2_[m][0] * qk_[m][0]) + (qk2_[m][1] * qk_[m][1]);
694 prod *= dKSq[m]*expKsq_[m];
698 dQ[n] = dQ[n] - increment;
708 fft().forwardTransform(qr_, qk_);
712 fft().forwardTransform(qr2_, qk2_);
715 for (
int n = 0; n < nParam ; ++n) {
719 for (
int m = 0; m < kSize_ ; ++m) {
720 prod = (qk2_[m][0] * qk_[m][0]) + (qk2_[m][1] * qk_[m][1]);
721 prod *= dKSq[m]*expKsq2_[m];
724 increment *= 0.5*bSq;
725 dQ[n] = dQ[n] - increment;
732 for (
int i = 0; i < nParam; ++i) {
733 stress_.append(-1.0*prefactor*dQ[i]);
virtual void setLength(double length)
Set the length of this block (only valid for thread model).
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.
Description of a regular grid of points in a periodic domain.
Fourier transform wrapper.
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.
bool isAllocated() const
Return true if the FftwDArray has been allocated, false otherwise.
Field of real double precision values on an FFT mesh.
Class to calculate and store properties of wavevectors.
Base template for UnitCell<D> classes, D=1, 2 or 3.
void stepBondBead(RField< D > const &qin, RField< D > &qout) const
Apply a bond operator for the bead model.
void setLength(double newLength)
Set or reset block length (only used in thread model).
void computeConcentrationBead(double prefactor)
Compute the concentration for this block, using the bead model.
void setKuhn(double kuhn)
Set or reset monomer statistical segment length.
void clearUnitCellData()
Clear all internal data that depends on the unit cell parameters.
void setupSolver(RField< D > const &w)
Set up the MDE solver for this block.
void allocate(double ds)
Allocate memory and set contour step size.
void computeStressBead(double prefactor)
Compute stress contribution for this block, using bead model.
int nBead() const
Get the number of beads in this block, in the bead model.
Propagator< D > & propagator(int directionId)
Get a Propagator for a specified direction.
void associate(Mesh< D > const &mesh, FFT< D > const &fft, UnitCell< D > const &cell, WaveList< D > &wavelist)
Create permanent associations with related objects.
double kuhn() const
Get monomer statistical segment length.
Propagator< D >::FieldT & cField()
Get the associated monomer concentration field.
void computeConcentrationThread(double prefactor)
Compute the concentration for this block, for the thread model.
void stepThread(RField< D > const &qin, RField< D > &qout) const
Compute one step of solution of MDE for the thread model.
void computeStressThread(double prefactor)
Compute stress contribution for this block, using thread model.
double ds() const
Get contour length step size.
void stepFieldBead(RField< D > &q) const
Apply the exponential field operator for the bead model.
void stepHalfBondBead(RField< D > const &qin, RField< D > &qout) const
Apply a half-bond operator for the bead model.
void stepBead(RField< D > const &qin, RField< D > &qout) const
Compute one step of solution of MDE for the bead model.
double length() const
Get the length of this block, in the thread model.
MDE solver for one direction of one block.
bool isHeadEnd() const
Is the head vertex a chain end?
bool isSolved() const
Has the modified diffusion equation been solved?
bool isTailEnd() const
Is the tail vertex a chain end?
bool isAllocated() const
Has memory been allocated for this propagator?
const FieldT & q(int i) const
Return q-field at specified step.
int capacity() const
Return allocated size.
A fixed capacity (static) contiguous array with a variable logical size.
void append(Data const &data)
Append data to the end of the array.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
bool isThread()
Is the thread model in use ?
bool isBead()
Is the bead model in use ?
Fields and FFTs for periodic boundary conditions (CPU)
Periodic fields and crystallography.
Real periodic fields, SCFT and PS-FTS (CPU).
PSCF package top-level namespace.