1#ifndef RPG_EINSTEIN_CRYSTAL_PERTURBATION_TPP
2#define RPG_EINSTEIN_CRYSTAL_PERTURBATION_TPP
4#include "EinsteinCrystalPerturbation.h"
5#include <rpg/fts/simulator/Simulator.h>
6#include <rpg/fts/VecOpFts.h>
8#include <prdc/cuda/resources.h>
15 using namespace Prdc::Cuda;
25 unperturbedHamiltonian_(0.0),
26 stateEcHamiltonian_(0.0),
27 stateUnperturbedHamiltonian_(0.0),
44 read(in,
"lambda", lambda_);
47 const int nMonomer = system().mixture().nMonomer();
48 epsilon_.allocate(nMonomer - 1);
49 for (
int i = 0; i < nMonomer - 1 ; ++i) {
55 readOptionalDArray(in,
"epsilon", epsilon_, nMonomer-1).isActive();
57 read(in,
"referenceFieldFileName", referenceFieldFileName_);
66 const int nMonomer = system().mixture().nMonomer();
68 meshDimensions = system().domain().mesh().dimensions();
71 w0_.allocate(nMonomer);
72 wc0_.allocate(nMonomer);
73 for (
int i = 0; i < nMonomer; ++i) {
74 w0_[i].allocate(meshDimensions);
75 wc0_[i].allocate(meshDimensions);
83 for (
int i = 0; i < nMonomer - 1 ; ++i) {
84 epsilon_[i] = -1.0 * simulator().chiEval(i);
89 for (
int i = 0; i < nMonomer - 1 ; ++i) {
94 FieldIo<D> const & fieldIo = system().domain().fieldIo();
96 w0_, system().domain().unitCell());
110 const int nMonomer = system().mixture().nMonomer();
111 const int meshSize = system().domain().mesh().size();
112 const IntVec<D> meshDimensions = system().domain().mesh().dimensions();
113 const double vSystem = system().domain().unitCell().volume();
114 const double vMonomer = system().mixture().vMonomer();
115 const double nMonomerSystem = vSystem / vMonomer;
119 ecHamiltonian_ = 0.0;
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]);
127 ecHamiltonian_ += prefactor * wSqure;
131 ecHamiltonian_ /= double(meshSize);
134 ecHamiltonian_ *= nMonomerSystem;
137 unperturbedHamiltonian_ = unperturbedHamiltonian;
139 return (1.0 - lambda_)*(ecHamiltonian_ - unperturbedHamiltonian_);
149 const IntVec<D> meshDimensions = system().domain().mesh().dimensions();
150 const int nMonomer = system().mixture().nMonomer();
151 const double vMonomer = system().mixture().vMonomer();
157 for (
int i = 0; i < nMonomer - 1; ++i) {
159 RField<D> const & Wc = simulator().wc(i);
160 prefactor = double(nMonomer) / (epsilon_[i] * vMonomer);
178 return unperturbedHamiltonian_ - ecHamiltonian_;
186 stateEcHamiltonian_ = ecHamiltonian_;
187 stateUnperturbedHamiltonian_ = unperturbedHamiltonian_;
195 ecHamiltonian_ = stateEcHamiltonian_;
196 unperturbedHamiltonian_ = stateUnperturbedHamiltonian_;
206 const int nMonomer = system().mixture().nMonomer();
210 for (i = 0; i < nMonomer; ++i) {
218 for (j = 0; j < nMonomer; ++j) {
220 vec = (cudaReal) simulator().chiEvecs(i, j)/nMonomer;
230 std::string filename =
"wc";
231 system().fieldIo().writeFieldsRGrid(filename, wc0_,
232 system().domain().unitCell(),
false);
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.
void allocate(IntVec< D > const &meshDimensions)
Allocate the underlying C array for an FFT grid.
Perturbation for Einstein crystal thermodynamic integration method.
virtual double df()
Compute and return derivative of free energy.
void setClassName(const char *className)
Set class name string.
EinsteinCrystalPerturbation(Simulator< D > &simulator)
Constructor.
virtual double hamiltonian(double unperturbedHamiltonian)
Compute and return the perturbation to the Hamiltonian.
virtual void setup()
Complete any required initialization.
virtual void incrementDc(DArray< RField< D > > &dc)
Modify the generalized forces to include perturbation.
virtual ~EinsteinCrystalPerturbation()
Destructor.
virtual void restoreState()
Restore any required internal state variables.
virtual void saveState()
Save any required internal state variables.
virtual void readParameters(std::istream &in)
Read parameters from archive.
File input/output operations and format conversions for fields.
void readFieldsRGrid(std::istream &in, DArray< RField< D > > &fields, UnitCell< D > &unitCell) const
Read array of RField objects (r-grid fields) from a stream.
Base class for additive perturbations of standard FTS Hamiltonian.
Field theoretic simulator (base class).
Dynamically allocatable contiguous array template.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
cudaReal innerProduct(DeviceArray< cudaReal > const &a, DeviceArray< cudaReal > const &b)
Compute inner product of two real arrays (GPU kernel wrapper).
void subVV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, DeviceArray< cudaReal > const &c, const int beginIdA, const int beginIdB, const int beginIdC, const int n)
Vector subtraction, a[i] = b[i] - c[i], kernel wrapper (cudaReal).
void addVcVc(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c, DeviceArray< cudaReal > const &d, cudaReal const e)
Vector addition w/ coefficient, a[i] = (b[i]*c) + (d[i]*e), kernel wrapper.
void mulEqS(DeviceArray< cudaReal > &a, cudaReal const b, const int beginIdA, const int n)
Vector multiplication in-place, a[i] *= b, kernel wrapper (cudaReal).
void eqS(DeviceArray< cudaReal > &a, cudaReal const b, const int beginIdA, const int n)
Vector assignment, a[i] = b, kernel wrapper (cudaReal).
void addEqVc(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c)
Vector addition in-place w/ coefficient, a[i] += b[i] * c, kernel wrapper.
PSCF package top-level namespace.
Utility classes for scientific computation.