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 one complete 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.