1#ifndef RPC_AM_ITERATOR_BASIS_TPP
2#define RPC_AM_ITERATOR_BASIS_TPP
11#include "AmIteratorBasis.h"
12#include <rpc/system/System.h>
13#include <rpc/solvers/Mixture.h>
14#include <rpc/field/Domain.h>
15#include <prdc/crystal/UnitCell.h>
16#include <pscf/inter/Interaction.h>
17#include <pscf/iterator/NanException.h>
18#include <util/containers/DArray.h>
58 UnitCell<D> const & unitCell = system().domain().unitCell();
71 flexibleParams_.clear();
72 for (
int i = 0; i < np; i++) {
73 flexibleParams_.append(
true);
77 if (nFlexibleParams() == 0) isFlexible_ =
false;
79 flexibleParams_.clear();
80 for (
int i = 0; i < np; i++) {
81 flexibleParams_.append(
false);
93 const int nMonomer = system().mixture().nMonomer();
94 interaction_.setNMonomer(nMonomer);
105 out <<
"Iterator times contributions:\n";
116 interaction_.update(system().interaction());
125 int AmIteratorBasis<D>::nElements()
127 const int nMonomer = system().mixture().nMonomer();
128 const int nBasis = system().domain().basis().nBasis();
130 int nEle = nMonomer*nBasis;
132 nEle += nFlexibleParams();
141 bool AmIteratorBasis<D>::hasInitialGuess()
142 {
return system().w().hasData(); }
150 const int nMonomer = system().mixture().nMonomer();
151 const int nBasis = system().domain().basis().nBasis();
152 const int nEle = nElements();
157 for (
int i = 0; i < nMonomer; i++) {
160 for (
int k = 0; k < nBasis; k++) {
161 state[begin + k] = field[k];
168 UnitCell<D> const & unitCell = system().domain().unitCell();
170 const int nParam = unitCell.nParameter();
171 begin = nMonomer*nBasis;
173 for (
int i = 0; i < nParam; i++) {
174 if (flexibleParams_[i]) {
175 state[begin + counter] = scaleStress_*parameters[i];
188 void AmIteratorBasis<D>::evaluate()
189 { system().compute(isFlexible_); }
197 const int nMonomer = system().mixture().nMonomer();
198 const int nBasis = system().domain().basis().nBasis();
199 const int n = nElements();
207 for (i = 0 ; i < n; ++i) {
212 for (i = 0; i < nMonomer; ++i) {
214 for (j = 0; j < nMonomer; ++j) {
215 chi = interaction_.chi(i,j);
216 p = interaction_.p(i,j);
219 if (system().h().hasData()) {
221 for (k = 0; k < nBasis; ++k) {
222 resid[begin + k] += chi*c[k] + p*(h[k] - w[k]);
225 for (k = 0; k < nBasis; ++k) {
226 resid[begin + k] += chi*c[k] - p*w[k];
233 double shift = -1.0 / interaction_.sumChiInverse();
234 if (system().mask().hasData()) {
236 for (i = 0; i < nMonomer; ++i) {
238 for (k = 0; k < nBasis; ++k) {
239 resid[begin + k] += shift*mask[k];
243 for (i = 0; i < nMonomer; ++i) {
244 resid[i*nBasis] += shift;
249 if (system().mixture().isCanonical()) {
250 for (i = 0; i < nMonomer; ++i) {
251 resid[i*nBasis] = 0.0;
267 double coeff = -1.0 * scaleStress_;
268 const int nParam = system().domain().unitCell().nParameter();
269 begin = nMonomer*nBasis;
271 for (i = 0; i < nParam ; i++) {
272 if (flexibleParams_[i]) {
273 resid[begin + counter] = coeff * stress(i);
289 const int nMonomer = system().mixture().nMonomer();
290 const int nBasis = system().domain().basis().nBasis();
295 for (
int i = 0; i < nMonomer; i++) {
301 for (
int i = 0; i < nMonomer; i++) {
303 for (
int k = 0; k < nBasis; k++) {
304 wFields[i][k] = newState[begin + k];
309 if (system().mixture().isCanonical()) {
312 for (
int i = 0; i < nMonomer; ++i) {
317 double chi, wAve, cAve;
318 for (
int i = 0; i < nMonomer; ++i) {
320 for (
int j = 0; j < nMonomer; ++j) {
321 chi = interaction_.chi(i,j);
322 cAve = system().c().basis(j)[0];
325 wFields[i][0] = wAve;
329 if (system().h().hasData()) {
330 for (
int i = 0; i < nMonomer; ++i) {
331 wFields[i][0] += system().h().basis(i)[0];
337 system().w().setBasis(wFields);
344 parameters = system().domain().unitCell().parameters();
347 const double coeff = 1.0 / scaleStress_;
348 const int nParam = system().domain().unitCell().nParameter();
349 const int begin = nMonomer*nBasis;
351 for (
int i = 0; i < nParam; i++) {
352 if (flexibleParams_[i]) {
353 parameters[i] = coeff * newState[begin + counter];
360 system().setUnitCell(parameters);
368 void AmIteratorBasis<D>::outputToLog()
370 if (isFlexible() && verbose() > 1) {
372 UnitCell<D> const & unitCell = system().domain().unitCell();
374 const int nMonomer = system().mixture().nMonomer();
375 const int nBasis = system().domain().basis().nBasis();
376 const int begin = nMonomer*nBasis;
378 for (
int i = 0; i < nParam; i++) {
379 if (flexibleParams_[i]) {
380 str = -1.0 * residual()[begin + counter] / scaleStress_;
382 <<
" Cell Param " << i <<
" = "
383 <<
Dbl(unitCell.parameters()[i], 15)
virtual void setup(bool isContinuation)
void readMixingParameters(std::istream &in, bool useLambdaRamp=true)
void outputTimers(std::ostream &out) const
void readErrorType(std::istream &in)
virtual void readParameters(std::istream &in)
int nParameter() const
Get the number of parameters in the unit cell.
Base template for UnitCell<D> classes, D=1, 2 or 3.
FSArrayParam< Type, N > & readOptionalFSArray(std::istream &in, const char *label, FSArray< Type, N > &array, int size)
Add and read an optional FSArray < Type, N > array parameter.
void readParameters(std::istream &in) override
Read all parameters and initialize.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
void outputTimers(std::ostream &out) const override
Output timing results to log file.
AmIteratorBasis(System< D > &system)
Constructor.
~AmIteratorBasis()
Destructor.
void setup(bool isContinuation) override
Setup iterator just before entering iteration loop.
Base class for iterative solvers for SCF equations.
Main class, representing a complete physical system.
int capacity() const
Return allocated size.
Dynamically allocatable contiguous array template.
void allocate(int capacity)
Allocate the underlying C array.
Wrapper for a double precision number, for formatted ostream output.
A fixed capacity (static) contiguous array with a variable logical size.
static std::ostream & file()
Get log ostream by reference.
void setClassName(const char *className)
Set class name string.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Periodic fields and crystallography.
Real periodic fields, SCFT and PS-FTS (CPU).
PSCF package top-level namespace.