8#include "BinaryRelaxIterator.h"
10#include <pscf/inter/Interaction.h>
11#include <util/misc/Log.h>
12#include <util/misc/Timer.h>
39 read(in,
"epsilon", epsilon_);
40 read(in,
"maxIter", maxItr_);
41 read(in,
"lambdaPlus", lambdaPlus_);
42 read(in,
"lambdaMinus", lambdaMinus_);
48 void BinaryRelaxIterator::allocate()
50 if (isAllocated_)
return;
59 wFieldsNew_.allocate(nm);
60 cFieldsNew_.allocate(nm);
61 for (
int i = 0; i < nm; ++i) {
62 wFieldsNew_[i].allocate(nx);
63 cFieldsNew_[i].allocate(nx);
65 dWNew_[i].allocate(nx);
70 void BinaryRelaxIterator::computeDW(
Array<WField> const & wOld,
76 double dWm, dWp, c0, c1, w0, w1, wm;
82 for (
int i = 0; i < nx; ++i) {
89 dWm = 0.5*( c0 - c1 - wm/chi);
97 dWNorm = (dWpNorm + dWmNorm)/
double(nx);
98 dWNorm = sqrt(dWNorm);
101 void BinaryRelaxIterator::updateWFields(
Array<WField> const & wOld,
113 for (j = 0; j < nx; j++){
116 wNew[0][j] = w0 + dW[0][j];
117 wNew[1][j] = w1 + dW[1][j];
122 double shift = wNew[nm-1][nx-1];
123 for (i = 0; i < nm; ++i) {
124 for (j = 0; j < nx; ++j) {
148 for (
int i = 0; i < np; ++i) {
150 if (ensemble == Species::Unknown) {
153 if (ensemble == Species::Open) {
154 isCanonical_ =
false;
161 double shift =
wFields()[nm-1][nx-1];
163 for (i = 0; i < nm; ++i) {
164 for (j = 0; j < nx; ++j) {
176 for (i = 0; i < maxItr_; ++i) {
178 <<
" , error = " << dWNorm_
181 if (dWNorm_ < epsilon_) {
185 Log::file() <<
"The epsilon is " << epsilon_<< std::endl;
192 << timerTotal.
time() <<
" s " << std::endl;
193 Log::file() <<
"Average time cost of each iteration: "
194 << timerTotal.
time()/i <<
" s " << std::endl;
202 computeDW(wFieldsNew_, cFieldsNew_, dWNew_, dWNormNew_);
206 while (dWNormNew_ > dWNorm_ && j < 3) {
209 << dWNormNew_ <<
", decreasing increment" << std::endl;
214 << lambdaPlus_ << std::endl;
216 << lambdaMinus_<< std::endl;
225 if (dWNormNew_ < dWNorm_) {
227 for (j = 0; j < nm; ++j) {
228 for (k = 0; k < nx; ++k) {
231 dW_[j][k] = dWNew_[j][k];
233 dWNorm_ = dWNormNew_;
236 Log::file() <<
"Iteration failed, norm = "
237 << dWNormNew_ << std::endl;
DMatrix< double > const & chi() const
Return the chi matrix by const reference.
Polymer & polymer(int id)
Get a polymer object.
int nMonomer() const
Get number of monomer types.
int nPolymer() const
Get number of polymer species.
virtual ~BinaryRelaxIterator()
Destructor.
int solve(bool isContinuation=false)
Iterate self-consistent field equations to solution.
void readParameters(std::istream &in)
Read all parameters and initialize.
BinaryRelaxIterator(System &system)
Constructor.
int nx() const
Get number of spatial grid points, including both endpoints.
Base class for iterative solvers for SCF equations.
void compute(DArray< WField > const &wFields, DArray< CField > &cFields)
Compute concentrations.
const Domain & domain() const
Get spatial domain (including grid info) by reference.
DArray< System::CField > & cFields()
Get array of all chemical potential fields.
DArray< System::WField > & wFields()
Get array of all chemical potential fields.
const System & system() const
Get parent System by reference.
const Mixture & mixture() const
Get Mixture by reference.
Main class in SCFT simulation of one system.
Interaction & interaction()
Get interaction (i.e., excess free energy) by reference.
CField & cField(int monomerId)
Get chemical potential field for a specific monomer type.
void computeFreeEnergy()
Compute free energy density and pressure for current fields.
WField & wField(int monomerId)
Get chemical potential field for a specific monomer type.
Ensemble
Statistical ensemble for number of molecules.
Array container class template.
void allocate(int capacity)
Allocate the underlying C array.
static std::ostream & file()
Get log ostream by reference.
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
void setClassName(const char *className)
Set class name string.
void start(TimePoint begin)
Start timing from an externally supplied time.
double time()
Return the accumulated time, in seconds.
void stop(TimePoint end)
Stop the clock at an externally supplied time.
#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.
PSCF package top-level namespace.
Utility classes for scientific computation.