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 <pscf/inter/Interaction.h>
14#include <pscf/iterator/NanException.h>
46 interaction_.setNMonomer(system().mixture().nMonomer());
52 int np = system().domain().unitCell().nParameter();
63 flexibleParams_.clear();
64 for (
int i = 0; i < np; i++) {
65 flexibleParams_.append(
true);
69 if (nFlexibleParams() == 0) isFlexible_ =
false;
71 flexibleParams_.clear();
72 for (
int i = 0; i < np; i++) {
73 flexibleParams_.append(
false);
87 out <<
"Iterator times contributions:\n";
98 interaction_.update(system().interaction());
117 for (
int i = 0; i < n; i++) {
119 if (std::isnan(a[i]) || std::isnan(b[i])) {
120 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
135 for (
int i = 0; i < n; i++) {
137 if (std::isnan(value)) {
138 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
141 if (fabs(value) > max) {
157 const int n = hists[0].capacity();
162 for (
int i = 0; i < n; i++) {
163 newbasis[i] = hists[0][i] - hists[1][i];
165 basis.append(newbasis);
177 for (
int i = 0; i < nHist; i++) {
178 for (
int j = 0; j < n; j++) {
180 trial[j] += coeffs[i] * -1 * basis[i][j];
187 void AmIteratorBasis<D>::addPredictedError(
DArray<double>& fieldTrial,
192 for (
int i = 0; i < n; i++) {
193 fieldTrial[i] += lambda * resTrial[i];
201 bool AmIteratorBasis<D>::hasInitialGuess()
202 {
return system().w().hasData(); }
206 int AmIteratorBasis<D>::nElements()
208 const int nMonomer = system().mixture().nMonomer();
209 const int nBasis = system().domain().basis().nBasis();
211 int nEle = nMonomer*nBasis;
213 nEle += nFlexibleParams();
223 const int nMonomer = system().mixture().nMonomer();
224 const int nBasis = system().domain().basis().nBasis();
228 for (
int i = 0; i < nMonomer; i++) {
229 for (
int k = 0; k < nBasis; k++) {
230 curr[i*nBasis+k] = currSys[i][k];
236 const int nParam = system().domain().unitCell().nParameter();
238 = system().domain().unitCell().parameters();
240 for (
int i = 0; i < nParam; i++) {
241 if (flexibleParams_[i]) {
242 curr[nMonomer*nBasis + counter] = scaleStress_*currParam[i];
253 void AmIteratorBasis<D>::evaluate()
257 system().compute(isFlexible_);
264 const int n = nElements();
265 const int nMonomer = system().mixture().nMonomer();
266 const int nBasis = system().domain().basis().nBasis();
269 for (
int i = 0 ; i < n; ++i) {
274 for (
int i = 0; i < nMonomer; ++i) {
275 for (
int j = 0; j < nMonomer; ++j) {
276 double chi = interaction_.chi(i,j);
277 double p = interaction_.p(i,j);
280 for (
int k = 0; k < nBasis; ++k) {
281 int idx = i*nBasis + k;
282 resid[idx] += chi*c[k] - p*w[k];
288 if (system().mask().hasData()) {
290 double sumChiInv = interaction_.sumChiInverse();
291 for (
int i = 0; i < nMonomer; ++i) {
292 for (
int k = 0; k < nBasis; ++k) {
293 int idx = i*nBasis + k;
294 resid[idx] -= mask[k] / sumChiInv;
301 if (system().h().hasData()) {
302 for (
int i = 0; i < nMonomer; ++i) {
303 for (
int j = 0; j < nMonomer; ++j) {
304 double p = interaction_.p(i,j);
306 for (
int k = 0; k < nBasis; ++k) {
307 int idx = i*nBasis + k;
308 resid[idx] += p * h[k];
315 if (!system().mixture().isCanonical()) {
316 if (!system().mask().hasData()) {
317 for (
int i = 0; i < nMonomer; ++i) {
318 resid[i*nBasis] -= 1.0 / interaction_.sumChiInverse();
323 for (
int i = 0; i < nMonomer; ++i) {
324 resid[i*nBasis] = 0.0;
330 const int nParam = system().domain().unitCell().nParameter();
342 for (
int i = 0; i < nParam ; i++) {
343 if (flexibleParams_[i]) {
344 double str = stress(i);
346 resid[nMonomer*nBasis + counter] = -1 * scaleStress_ * str;
360 const int nMonomer = system().mixture().nMonomer();
361 const int nBasis = system().domain().basis().nBasis();
367 for (
int i = 0; i < nMonomer; i++) {
369 for (
int k = 0; k < nBasis; k++)
371 wField[i][k] = newGuess[i*nBasis + k];
375 if (system().mixture().isCanonical()) {
377 for (
int i = 0; i < nMonomer; ++i) {
379 for (
int j = 0; j < nMonomer; ++j) {
380 chi = interaction_.chi(i,j);
381 wField[i][0] += chi * system().c().basis(j)[0];
385 if (system().h().hasData()) {
386 for (
int i = 0; i < nMonomer; ++i) {
387 wField[i][0] += system().h().basis(i)[0];
391 system().w().setBasis(wField);
394 const int nParam = system().domain().unitCell().nParameter();
395 const int begin = nMonomer*nBasis;
398 parameters = system().domain().unitCell().parameters();
400 double coeff = 1.0 / scaleStress_;
402 for (
int i = 0; i < nParam; i++) {
403 if (flexibleParams_[i]) {
404 parameters[i] = coeff * newGuess[begin + counter];
410 system().setUnitCell(parameters);
416 void AmIteratorBasis<D>::outputToLog()
418 if (isFlexible() && verbose() > 1) {
419 const int nParam = system().domain().unitCell().nParameter();
420 const int nMonomer = system().mixture().nMonomer();
421 const int nBasis = system().domain().basis().nBasis();
423 for (
int i = 0; i < nParam; i++) {
424 if (flexibleParams_[i]) {
425 double str = residual()[nMonomer*nBasis + counter] /
426 (-1.0 * scaleStress_);
428 <<
" Cell Param " << i <<
" = "
429 <<
Dbl(system().domain().unitCell().parameters()[i], 15)
void readErrorType(std::istream &in)
Exception thrown when not-a-number (NaN) is encountered.
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 setup(bool isContinuation)
Setup iterator just before entering iteration loop.
void outputTimers(std::ostream &out) const
Output timing results to log file.
void setClassName(const char *className)
Set class name string.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
AmIteratorBasis(System< D > &system)
Constructor.
~AmIteratorBasis()
Destructor.
void readParameters(std::istream &in)
Read all parameters and initialize.
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.
Class for storing history of previous values in an array.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Real periodic fields, SCFT and PS-FTS (CPU).
PSCF package top-level namespace.
float product(float a, float b)
Product for float Data.