1#ifndef RPC_EINSTEIN_CRYSTAL_PERTURBATION_TPP
2#define RPC_EINSTEIN_CRYSTAL_PERTURBATION_TPP
4#include "EinsteinCrystalPerturbation.h"
5#include <rpc/fts/simulator/Simulator.h>
6#include <rpc/system/System.h>
7#include <rpc/solvers/Mixture.h>
8#include <rpc/field/Domain.h>
9#include <prdc/cpu/RField.h>
10#include <util/containers/DArray.h>
26 unperturbedHamiltonian_(0.0),
27 stateEcHamiltonian_(0.0),
28 stateUnperturbedHamiltonian_(0.0),
48 const int nMonomer =
system().mixture().nMonomer();
49 epsilon_.allocate(nMonomer - 1);
50 for (
int i = 0; i < nMonomer - 1 ; ++i) {
58 read(in,
"referenceFieldFileName", referenceFieldFileName_);
67 const int nMonomer =
system().mixture().nMonomer();
69 meshDimensions =
system().domain().mesh().dimensions();
72 w0_.allocate(nMonomer);
73 wc0_.allocate(nMonomer);
74 for (
int i = 0; i < nMonomer; ++i) {
75 w0_[i].allocate(meshDimensions);
76 wc0_[i].allocate(meshDimensions);
85 for (
int i = 0; i < nMonomer - 1 ; ++i) {
86 epsilon_[i] = -1.0 *
simulator().chiEval(i);
91 for (
int i = 0; i < nMonomer - 1 ; ++i) {
102 computeWcReference();
113 const int nMonomer =
system().mixture().nMonomer();
114 const int meshSize =
system().domain().mesh().size();
115 const double vSystem =
system().domain().unitCell().volume();
116 const double vMonomer =
system().mixture().vMonomer();
117 const double nMonomerSystem = vSystem / vMonomer;
119 ecHamiltonian_ = 0.0;
121 for (
int j = 0; j < nMonomer - 1; ++j) {
123 prefactor = double(nMonomer)/(2.0 * epsilon_[j]);
124 for (
int i = 0; i < meshSize; ++i) {
125 w = Wc[i] - wc0_[j][i];
126 ecHamiltonian_ += prefactor*w*w;
131 ecHamiltonian_ /= double(meshSize);
134 ecHamiltonian_ *= nMonomerSystem;
137 unperturbedHamiltonian_ = unperturbedHamiltonian;
139 return (1.0 -
lambda_)*(ecHamiltonian_ - unperturbedHamiltonian_);
149 const int meshSize =
system().domain().mesh().size();
150 const int nMonomer =
system().mixture().nMonomer();
151 const double vMonomer =
system().mixture().vMonomer();
152 double DcEC, prefactor;
155 for (
int i = 0; i < nMonomer - 1; ++i) {
158 prefactor = double(nMonomer) / (epsilon_[i] * vMonomer);
161 for (
int k = 0; k < meshSize; ++k) {
162 DcEC = prefactor * (Wc[k] - wc0_[i][k]);
163 Dc[k] += (1.0 -
lambda_) * (DcEC - Dc[k]);
174 return unperturbedHamiltonian_ - ecHamiltonian_;
182 stateEcHamiltonian_ = ecHamiltonian_;
183 stateUnperturbedHamiltonian_ = unperturbedHamiltonian_;
191 ecHamiltonian_ = stateEcHamiltonian_;
192 unperturbedHamiltonian_ = stateUnperturbedHamiltonian_;
200 void EinsteinCrystalPerturbation<D>::computeWcReference()
202 const int nMonomer = system().mixture().nMonomer();
203 const int meshSize = system().domain().mesh().size();
207 for (j = 0; j < nMonomer; ++j) {
211 for (i = 0; i < meshSize; ++i) {
216 for (k = 0; k < nMonomer; ++k) {
217 double vec = simulator().chiEvecs(j, k)/double(nMonomer);
220 RField<D>
const & Wr = w0_[k];
221 for (i = 0; i < meshSize; ++i) {
228 Log::file() <<
"wc " << wc0_.capacity() <<
"\n";
229 for (i = 0; i < 10; ++i) {
230 Log::file() <<
"wc_1 " << wc0_[0][i] <<
"\n";
231 Log::file() <<
"wc_2 " << wc0_[1][i] <<
"\n";
An IntVec<D, T> is a D-component vector of elements of integer type T.
Field of real double precision values on an FFT mesh.
Base template for UnitCell<D> classes, D=1, 2 or 3.
virtual void restoreState()
Restore any required internal state variables.
DArrayParam< Type > & readOptionalDArray(std::istream &in, const char *label, DArray< Type > &array, int n)
Add and read an optional DArray < Type > parameter.
virtual void readParameters(std::istream &in)
Read parameters from archive.
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
virtual void saveState()
Save any required internal state variables.
void setClassName(const char *className)
Set class name string.
EinsteinCrystalPerturbation(Simulator< D > &simulator)
Constructor.
virtual double df()
Compute and return derivative of free energy.
virtual void incrementDc(DArray< RField< D > > &dc)
Modify the generalized forces to include perturbation.
virtual void setup()
Complete any required initialization.
virtual double hamiltonian(double unperturbedHamiltonian)
Compute and return the perturbation to the Hamiltonian.
virtual ~EinsteinCrystalPerturbation()
Destructor.
File input/output operations and format conversions for fields.
bool readFieldsRGrid(std::istream &in, DArray< RField< D > > &fields, UnitCell< D > &unitCell) const override
Read array of RField objects (r-grid fields) from a stream.
Perturbation(Simulator< D > &simulator)
Constructor.
Simulator< D > const & simulator() const
Get parent Simulator<D> by const reference.
double lambda_
Strength of the perturbation.
System< D > const & system() const
Get parent System<D> by const reference.
Field theoretic simulator (base class).
Dynamically allocatable contiguous array template.
static std::ostream & file()
Get log ostream by reference.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Real periodic fields, SCFT and PS-FTS (CPU).
PSCF package top-level namespace.