1#ifndef RPC_SIMULATOR_TPP
2#define RPC_SIMULATOR_TPP
13#include <rpc/system/System.h>
14#include <rpc/fts/compressor/Compressor.h>
15#include <rpc/fts/compressor/CompressorFactory.h>
16#include <rpc/fts/perturbation/Perturbation.h>
17#include <rpc/fts/perturbation/PerturbationFactory.h>
18#include <rpc/fts/ramp/Ramp.h>
19#include <rpc/fts/ramp/RampFactory.h>
21#include <pscf/inter/Interaction.h>
23#include <util/misc/Timer.h>
24#include <util/random/Random.h>
28#include <gsl/gsl_eigen.h>
53 compressorFactoryPtr_(nullptr),
54 compressorPtr_(nullptr),
55 perturbationFactoryPtr_(nullptr),
56 perturbationPtr_(nullptr),
57 rampFactoryPtr_(nullptr),
73 if (compressorFactoryPtr_) {
74 delete compressorFactoryPtr_;
76 if (compressorPtr_ ) {
77 delete compressorPtr_;
79 if (perturbationFactoryPtr_) {
80 delete perturbationFactoryPtr_;
82 if (perturbationPtr_) {
83 delete perturbationPtr_;
85 if (rampFactoryPtr_) {
86 delete rampFactoryPtr_;
101 const int nMonomer =
system().mixture().nMonomer();
104 chiP_.allocate(nMonomer, nMonomer);
105 chiEvals_.allocate(nMonomer);
106 chiEvecs_.allocate(nMonomer, nMonomer);
107 sc_.allocate(nMonomer);
110 wc_.allocate(nMonomer);
111 cc_.allocate(nMonomer);
113 for (
int i = 0; i < nMonomer; ++i) {
114 wc_[i].allocate(dimensions);
115 cc_[i].allocate(dimensions);
119 dc_.allocate(nMonomer-1);
120 for (
int i = 0; i < nMonomer - 1; ++i) {
121 dc_[i].allocate(dimensions);
125 if (!
state_.isAllocated) {
127 state_.allocate(nMonomer, dimensions);
160 {
UTIL_THROW(
"Error: Unimplemented function Simulator<D>::simulate"); }
167 std::string classname,
168 std::string filename)
169 {
UTIL_THROW(
"Error: Unimplemented function Simulator<D>::analyze"); }
198 const int nMonomer = mixture.
nMonomer();
199 const int meshSize = domain.
mesh().size();
210 for (
int i = 0; i < np; ++i) {
211 polymerPtr = &mixture.
polymer(i);
212 phi = polymerPtr->
phi();
213 mu = polymerPtr->
mu();
215 length = polymerPtr->
length();
217 length = (double)polymerPtr->
nBead();
221 lnQ += phi*( -mu + 1.0 )/length;
230 for (
int i = 0; i < ns; ++i) {
231 solventPtr = &mixture.
solvent(i);
232 phi = solventPtr->
phi();
233 mu = solventPtr->
mu();
234 size = solventPtr->
size();
237 lnQ += phi*( -mu + 1.0 )/size;
244 for (
int i = 0; i < meshSize; ++i) {
245 lnQ += Wc[i]/double(meshSize);
253 double prefactor, w, s;
255 for (j = 0; j < nMonomer - 1; ++j) {
257 prefactor = -0.5*double(nMonomer)/chiEvals_[j];
259 for (i = 0; i < meshSize; ++i) {
266 HW /= double(meshSize);
269 HW += 0.5*sc_[nMonomer - 1];
272 const double vSystem = domain.
unitCell().volume();
273 const double vMonomer = mixture.
vMonomer();
274 const double nMonomerSystem = vSystem / vMonomer;
296 const int nMonomer =
system().mixture().nMonomer();
298 double d = 1.0/double(nMonomer);
304 for (i = 0; i < nMonomer; ++i) {
305 for (j = 0; j < nMonomer; ++j) {
314 for (i = 0; i < nMonomer; ++i) {
315 for (j = 0; j < nMonomer; ++j) {
317 for (k = 0; k < nMonomer; ++k) {
318 T(i,j) += chi(i,k)*P(k,j);
326 for (i = 0; i < nMonomer; ++i) {
327 for (j = 0; j < nMonomer; ++j) {
329 for (k = 0; k < nMonomer; ++k) {
330 chiP(i,j) += P(i,k)*T(k,j);
339 gsl_matrix* A = gsl_matrix_alloc(nMonomer, nMonomer);
342 for (i = 0; i < nMonomer; ++i) {
343 for (j = 0; j < nMonomer; ++j) {
344 gsl_matrix_set(A, i, j, chiP(i, j));
349 gsl_eigen_symmv_workspace* work = gsl_eigen_symmv_alloc(nMonomer);
350 gsl_vector* Avals = gsl_vector_alloc(nMonomer);
351 gsl_matrix* Avecs = gsl_matrix_alloc(nMonomer, nMonomer);
353 error = gsl_eigen_symmv(A, Avals, Avecs, work);
365 for (i = 0; i < nMonomer; ++i) {
366 val = gsl_vector_get(Avals, i);
367 if (std::abs(val) < 1.0E-8) {
374 for (j = 0; j < nMonomer; ++j) {
375 chiEvecs_(k, j) = gsl_matrix_get(Avecs, j, i);
377 if (chiEvecs_(k, 0) < 0.0) {
378 for (j = 0; j < nMonomer; ++j) {
379 chiEvecs_(k, j) = -chiEvecs_(k, j);
391 for (j = 0; j < nMonomer; ++j) {
392 chiEvecs_(i, j) = gsl_matrix_get(Avecs, j, iNull);
394 if (chiEvecs_(i, 0) < 0) {
395 for (j = 0; j < nMonomer; ++j) {
396 chiEvecs_(i, j) = -chiEvecs_(i, j);
401 double vec, norm, prefactor;
402 for (i = 0; i < nMonomer; ++i) {
404 for (j = 0; j < nMonomer; ++j) {
405 vec = chiEvecs_(i, j);
408 prefactor = sqrt(
double(nMonomer)/norm );
409 for (j = 0; j < nMonomer; ++j) {
410 chiEvecs_(i, j) *= prefactor;
415 for (j = 0; j < nMonomer; ++j) {
416 UTIL_CHECK(std::abs(chiEvecs_(nMonomer-1, j) - 1.0) < 1.0E-8);
422 for (i = 0; i < nMonomer; ++i) {
424 for (j = 0; j < nMonomer; ++j) {
427 s[i] = s[i]/double(nMonomer);
431 for (i = 0; i < nMonomer; ++i) {
433 for (j = 0; j < nMonomer; ++j) {
434 sc_[i] += chiEvecs_(i,j)*s[j];
436 sc_[i] = sc_[i]/double(nMonomer);
441 for (i = 0; i < nMonomer; ++i) {
442 Log::file() <<
"Eigenpair " << i <<
"\n";
443 Log::file() <<
"value = " << chiEvals_[i] <<
"\n";
445 for (j = 0; j < nMonomer; ++j) {
449 Log::file() <<
" sc[i] = " << sc_[i] << std::endl;
464 const int nMonomer =
system().mixture().nMonomer();
465 const int meshSize =
system().domain().mesh().size();
469 for (j = 0; j < nMonomer; ++j) {
473 for (i = 0; i < meshSize; ++i) {
478 for (k = 0; k < nMonomer; ++k) {
479 double vec = chiEvecs_(j, k)/double(nMonomer);
483 for (i = 0; i < meshSize; ++i) {
505 const int nMonomer =
system().mixture().nMonomer();
506 const int meshSize =
system().domain().mesh().size();
510 for (i = 0; i < nMonomer; ++i) {
514 for (k = 0; k < meshSize; ++k) {
519 for (j = 0; j < nMonomer; ++j) {
521 double vec = chiEvecs_(i, j);
524 for (k = 0; k < meshSize; ++k) {
546 const int meshSize =
system().domain().mesh().size();
547 const int nMonomer =
system().mixture().nMonomer();
548 const double vMonomer =
system().mixture().vMonomer();
549 const double a = 1.0/vMonomer;
555 for (i = 0; i < nMonomer - 1; ++i) {
559 b = -1.0*double(nMonomer)/chiEvals_[i];
562 for (k = 0; k < meshSize; ++k) {
563 Dc[k] = a*( b*(Wc[k] - s) + Cc[k] );
589 int nMonomer =
system().mixture().nMonomer();
591 for (
int i = 0; i < nMonomer; ++i) {
600 for (
int i = 0; i < nMonomer; ++i) {
609 for (
int i = 0; i < nMonomer - 1; ++i) {
615 if (
state_.needsHamiltonian){
641 const int nMonomer =
system().mixture().nMonomer();
647 if (
state_.needsHamiltonian){
655 for (
int i = 0; i < nMonomer; ++i) {
661 for (
int i = 0; i < nMonomer; ++i) {
668 for (
int i = 0; i < nMonomer - 1; ++i) {
688 {
state_.hasData =
false; }
699 compressorPtr_->outputTimers(out);
709 out <<
"MDE counter "
710 << compressorPtr_->mdeCounter() << std::endl;
721 compressorPtr_->clearTimers();
786 perturbationPtr_ = ptr;
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.
SolventT & solvent(int id)
Get a solvent solver object.
PolymerT & polymer(int id)
Get a polymer solver object by non-const reference.
Factory for subclasses of Compressor.
Spatial domain for a periodic structure with real fields, on a CPU.
Mesh< D > & mesh()
Get the Mesh by non-const reference.
UnitCell< D > & unitCell()
Get the UnitCell by non-const reference.
Solver and descriptor for a mixture of polymers and solvents.
int nPolymer() const
Get number of polymer species.
int nMonomer() const
Get number of monomer types.
int nSolvent() const
Get number of solvent (point particle) species.
double vMonomer() const
Get monomer reference volume (set to 1.0 by default).
Factory for subclasses of Perturbation.
Base class for additive perturbations of standard FTS Hamiltonian.
Descriptor and solver for one polymer species.
double phi() const
Get the overall volume fraction for this species.
int nBead() const
Total number of beads in the polymer (bead model).
double mu() const
Get the chemical potential for this species (units kT=1).
double length() const
Sum of the lengths of all blocks in the polymer (thread model).
Factory for subclasses of Ramp.
Class that varies parameters during a simulation (abstract).
bool hasCc_
Have eigen-components of the current c fields been computed ?
void readPerturbation(std::istream &in, bool &isEnd)
Optionally read an associated perturbation.
Perturbation< D > & perturbation()
Get the perturbation factory by non-const reference.
System< D > & system()
Get parent system by reference.
bool hasWc_
Have eigen-components of the current w fields been computed ?
void computeDc()
Compute functional derivatives of the Hamiltonian.
SimState< D > state_
Previous state saved during at the beginning of a step.
bool hasHamiltonian_
Has the Hamiltonian been computed for the current w and c fields?
virtual void analyze(int min, int max, std::string classname, std::string filename)
Read and analyze a trajectory file.
void setClassName(const char *className)
Set class name string.
double perturbationHamiltonian() const
Get the perturbation to the standard Hamiltonian (if any).
void clearState()
Clear the saved copy of the fts state.
PerturbationFactory< D > & perturbationFactory()
Get the perturbation factory by reference.
void computeWc()
Compute eigenvector components of the current w fields.
DArray< RField< D > > const & cc() const
Get all eigenvector components of the current c fields.
void setPerturbation(Perturbation< D > *ptr)
Set the associated perturbation.
void clearData()
Clear field eigen-components and hamiltonian components.
bool hasDc() const
Are the current d fields valid ?
void analyzeChi()
Perform eigenvalue analysis of projected chi matrix.
double fieldHamiltonian_
Field contribution (H_W) to Hamiltonian.
virtual void outputMdeCounter(std::ostream &out) const
Output MDE counter.
long seed_
Random number generator seed.
virtual void readParameters(std::istream &in)
Read parameters for a simulation.
void allocate()
Allocate required memory.
double idealHamiltonian_
Ideal gas contribution (-lnQ) to Hamiltonian H[W].
virtual void simulate(int nStep)
Perform a field theoretic Monte-Carlo simulation.
void setRamp(Ramp< D > *ptr)
Set the associated ramp.
DArray< RField< D > > wc_
Eigenvector components of w fields on a real space grid.
DArray< RField< D > > dc_
Components of d fields on a real space grid.
bool hasCc() const
Are eigen-components of current c fields valid ?
RampFactory< D > & rampFactory()
Get the ramp factory by reference.
double hamiltonian_
Total field theoretic Hamiltonian H[W] (extensive value).
double perturbationHamiltonian_
Perturbation to the standard Hamiltonian (if any).
virtual void clearTimers()
Clear timers.
bool hasWc() const
Are eigen-components of current w fields valid ?
CompressorFactory< D > & compressorFactory()
Get the compressor factory by reference.
bool hasDc_
Have functional derivatives of H[W] been computed ?
void restoreState()
Restore the saved copy of the fts move state.
void saveState()
Save a copy of the fts move state.
long iTotalStep_
Step counter - total number of attempted BD or MC steps.
double hamiltonian() const
Get the Hamiltonian used in PS-FTS.
bool hasCompressor() const
Does this Simulator have a Compressor?
void readRandomSeed(std::istream &in)
Optionally read a random number generator seed.
bool hasPerturbation() const
Does this Simulator have a Perturbation?
DArray< RField< D > > cc_
Eigenvector components of c fields on a real space grid.
virtual void outputTimers(std::ostream &out) const
Output timing results.
void readCompressor(std::istream &in, bool &isEnd)
Read the compressor block of the parameter file.
double fieldHamiltonian() const
Get the quadratic field contribution to the Hamiltonian.
long iStep_
Step counter - attempted steps for which compressor converges.
void readRamp(std::istream &in, bool &isEnd)
Optionally read an associated ramp.
bool hasHamiltonian() const
Has the Hamiltonian been computed for current w and c fields?
DArray< RField< D > > const & dc() const
Get all of the current d fields.
Random random_
Random number generator.
DArray< RField< D > > const & wc() const
Get all eigenvector components of the current w fields.
void computeHamiltonian()
Compute the Hamiltonian used in PS-FTS.
double idealHamiltonian() const
Get ideal gas contribution to the Hamiltonian.
Random & random()
Get random number generator by reference.
void computeCc()
Compute eigenvector components of the current c fields.
Solver and descriptor for a solvent species.
double phi() const
Get the overall volume fraction for this species.
double mu() const
Get the chemical potential for this species (units kT=1).
double size() const
Get the size (number of monomers) in this solvent.
Main class, representing one complete system.
Dynamically allocatable contiguous array template.
void allocate(int capacity)
Allocate the underlying C array.
Dynamically allocated Matrix.
void allocate(int capacity1, int capacity2)
Allocate memory for a matrix.
static std::ostream & file()
Get log ostream by reference.
std::string indent() const
Return indent string for this object (string of spaces).
static bool echo()
Get echo parameter.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
std::string className() const
Get class name string.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
bool isThread()
Is the thread model in use ?
Real periodic fields, SCFT and PS-FTS (CPU).
PSCF package top-level namespace.