PSCF v1.3.3
rpc/fts/perturbation/EinsteinCrystalPerturbation.tpp
1#ifndef RPC_EINSTEIN_CRYSTAL_PERTURBATION_TPP
2#define RPC_EINSTEIN_CRYSTAL_PERTURBATION_TPP
3
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>
11#include <util/global.h>
12
13namespace Pscf {
14namespace Rpc {
15
16 using namespace Util;
17
18 /*
19 * Constructor.
20 */
21 template <int D>
25 ecHamiltonian_(0.0),
26 unperturbedHamiltonian_(0.0),
27 stateEcHamiltonian_(0.0),
28 stateUnperturbedHamiltonian_(0.0),
29 hasEpsilon_(false)
30 { setClassName("EinsteinCrystal"); }
31
32 /*
33 * Destructor.
34 */
35 template <int D>
38
39 /*
40 * Read parameters from stream, empty default implementation.
41 */
42 template <int D>
44 {
45 read(in, "lambda", lambda_);
46
47 // Allocate and initialize epsilon_ array
48 const int nMonomer = system().mixture().nMonomer();
49 epsilon_.allocate(nMonomer - 1);
50 for (int i = 0; i < nMonomer - 1 ; ++i) {
51 epsilon_[i] = 0.0;
52 }
53
54 // Optionally read the parameters used in Einstein crystal integration
55 hasEpsilon_ =
56 readOptionalDArray(in, "epsilon", epsilon_, nMonomer-1).isActive();
57
58 read(in, "referenceFieldFileName", referenceFieldFileName_);
59 }
60
61 /*
62 * Setup before simulation.
63 */
64 template <int D>
66 {
67 const int nMonomer = system().mixture().nMonomer();
68 const IntVec<D>
69 meshDimensions = system().domain().mesh().dimensions();
70
71 // Allocate memory for reference field
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);
77 }
78
79 /*
80 * If the user did not input values for epsilon, set
81 * the values of epsilon_ to the nontrivial -1.0 * eigenvalues
82 * of the projected chi matrix by default.
83 */
84 if (!hasEpsilon_){
85 for (int i = 0; i < nMonomer - 1 ; ++i) {
86 epsilon_[i] = -1.0 * simulator().chiEval(i);
87 }
88 }
89
90 // Check that all epsilon values are positive
91 for (int i = 0; i < nMonomer - 1 ; ++i) {
92 UTIL_CHECK(epsilon_[i] > 0.0);
93 }
94
95 // Read in reference field from a file
96 UnitCell<D> tempUnitCell;
97 FieldIo<D> const & fieldIo = system().domain().fieldIo();
98 fieldIo.readFieldsRGrid(referenceFieldFileName_,
99 w0_, tempUnitCell);
100
101 // Compute eigenvector components of the reference field
102 computeWcReference();
103 }
104
105 /*
106 * Compute and return perturbation to Hamiltonian.
107 */
108 template <int D>
109 double
110 EinsteinCrystalPerturbation<D>::hamiltonian(double unperturbedHamiltonian)
111 {
112 // Compute Einstein crystal Hamiltonian
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;
118 double prefactor, w;
119 ecHamiltonian_ = 0.0;
120
121 for (int j = 0; j < nMonomer - 1; ++j) {
122 RField<D> const & Wc = simulator().wc(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;
127 }
128 }
129
130 // Normalize EC hamiltonian to equal a value per monomer
131 ecHamiltonian_ /= double(meshSize);
132
133 // Compute EC hamiltonian of system
134 ecHamiltonian_ *= nMonomerSystem;
135
136 // Set unperturbedHamiltonian_ member variable
137 unperturbedHamiltonian_ = unperturbedHamiltonian;
138
139 return (1.0 - lambda_)*(ecHamiltonian_ - unperturbedHamiltonian_);
140 }
141
142 /*
143 * Modify functional derivatives, empty default implementation.
144 */
145 template <int D>
146 void
148 {
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;
153
154 // Loop over composition eigenvectors (exclude the last)
155 for (int i = 0; i < nMonomer - 1; ++i) {
156 RField<D>& Dc = dc[i];
157 RField<D> const & Wc = simulator().wc(i);
158 prefactor = double(nMonomer) / (epsilon_[i] * vMonomer);
159
160 // Loop over grid points, compute composite derivative
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]);
164 }
165
166 }
167 }
168
169 /*
170 * Compute and return derivative of free energy with respect to lambda.
171 */
172 template <int D>
174 return unperturbedHamiltonian_ - ecHamiltonian_;
175 }
176
177 /*
178 * Save any required internal state variables.
179 */
180 template <int D>
182 stateEcHamiltonian_ = ecHamiltonian_;
183 stateUnperturbedHamiltonian_ = unperturbedHamiltonian_;
184 }
185
186 /*
187 * Save any required internal state variables.
188 */
189 template <int D>
191 ecHamiltonian_ = stateEcHamiltonian_;
192 unperturbedHamiltonian_ = stateUnperturbedHamiltonian_;
193 }
194
195 /*
196 * Compute the eigenvector components of the w fields, using the
197 * eigenvectors chiEvecs of the projected chi matrix as a basis.
198 */
199 template <int D>
200 void EinsteinCrystalPerturbation<D>::computeWcReference()
201 {
202 const int nMonomer = system().mixture().nMonomer();
203 const int meshSize = system().domain().mesh().size();
204 int i, j, k;
205
206 // Loop over eigenvectors (j is an eigenvector index)
207 for (j = 0; j < nMonomer; ++j) {
208
209 // Loop over grid points to zero out field wc_[j]
210 RField<D>& Wc = wc0_[j];
211 for (i = 0; i < meshSize; ++i) {
212 Wc[i] = 0.0;
213 }
214
215 // Loop over monomer types (k is a monomer index)
216 for (k = 0; k < nMonomer; ++k) {
217 double vec = simulator().chiEvecs(j, k)/double(nMonomer);
218
219 // Loop over grid points
220 RField<D> const & Wr = w0_[k];
221 for (i = 0; i < meshSize; ++i) {
222 Wc[i] += vec*Wr[i];
223 }
224 }
225 }
226
227 #if 0
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";
232 }
233 #endif
234 }
235
236}
237}
238#endif
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
Field of real double precision values on an FFT mesh.
Definition cpu/RField.h:29
Base template for UnitCell<D> classes, D=1, 2 or 3.
Definition UnitCell.h:56
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.
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.
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.
Definition DArray.h:32
static std::ostream & file()
Get log ostream by reference.
Definition Log.cpp:59
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
Real periodic fields, SCFT and PS-FTS (CPU).
Definition param_pc.dox:2
PSCF package top-level namespace.