1#ifndef RPG_AM_ITERATOR_BASIS_TPP
2#define RPG_AM_ITERATOR_BASIS_TPP
11#include "AmIteratorBasis.h"
12#include <rpg/system/System.h>
14#include <prdc/crystal/UnitCell.h>
15#include <prdc/crystal/Basis.h>
17#include <pscf/inter/Interaction.h>
18#include <pscf/iterator/NanException.h>
52 interaction_.setNMonomer(system().mixture().nMonomer());
58 int np = system().domain().unitCell().nParameter();
69 flexibleParams_.clear();
70 for (
int i = 0; i < np; i++) {
71 flexibleParams_.append(
true);
75 if (nFlexibleParams() == 0) isFlexible_ =
false;
77 flexibleParams_.clear();
78 for (
int i = 0; i < np; i++) {
79 flexibleParams_.append(
false);
93 out <<
"Iterator times contributions:\n";
107 interaction_.update(system().interaction());
127 for (
int i=0; i < n; ++i) {
129 if (std::isnan(a[i]) || std::isnan(b[i])) {
130 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
145 for (
int i = 0; i < n; i++) {
147 if (std::isnan(value)) {
148 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
151 if (fabs(value) > max)
166 const int n = hists[0].capacity();
171 for (
int i = 0; i < n; i++) {
172 newbasis[i] = hists[0][i] - hists[1][i];
175 basis.append(newbasis);
187 for (
int i = 0; i < nHist; i++) {
188 for (
int j = 0; j < n; j++) {
190 trial[j] += coeffs[i] * -1 * basis[i][j];
197 void AmIteratorBasis<D>::addPredictedError(
DArray<double>& fieldTrial,
202 for (
int i = 0; i < n; i++) {
203 fieldTrial[i] += lambda * resTrial[i];
209 bool AmIteratorBasis<D>::hasInitialGuess()
211 return system().w().hasData();
216 int AmIteratorBasis<D>::nElements()
218 const int nMonomer = system().mixture().nMonomer();
219 const int nBasis = system().domain().basis().nBasis();
221 int nEle = nMonomer*nBasis;
224 nEle += nFlexibleParams();
236 const int nMonomer = system().mixture().nMonomer();
237 const int nBasis = system().domain().basis().nBasis();
240 for (
int i = 0; i < nMonomer; i++) {
241 for (
int k = 0; k < nBasis; k++) {
242 curr[i*nBasis+k] = currSys[i][k];
247 const int begin = nMonomer*nBasis;
248 const int nParam = system().domain().unitCell().nParameter();
250 = system().domain().unitCell().parameters();
252 for (
int i = 0; i < nParam; i++) {
253 if (flexibleParams_[i]) {
254 curr[begin + counter] = scaleStress_ * parameters[i];
265 void AmIteratorBasis<D>::evaluate()
269 system().compute(isFlexible_);
276 UTIL_CHECK(system().domain().basis().isInitialized());
277 const int nMonomer = system().mixture().nMonomer();
278 const int nBasis = system().domain().basis().nBasis();
279 const int n = nElements();
282 for (
int i = 0 ; i < n; ++i) {
287 for (
int i = 0; i < nMonomer; ++i) {
288 for (
int j = 0; j < nMonomer; ++j) {
289 double chi = interaction_.chi(i,j);
290 double p = interaction_.p(i,j);
293 for (
int k = 0; k < nBasis; ++k) {
294 int idx = i*nBasis + k;
295 resid[idx] += chi*c[k] - p*w[k];
301 if (system().mask().hasData()) {
303 double sumChiInv = interaction_.sumChiInverse();
304 for (
int i = 0; i < nMonomer; ++i) {
305 for (
int k = 0; k < nBasis; ++k) {
306 int idx = i*nBasis + k;
307 resid[idx] -= mask[k] / sumChiInv;
314 if (system().h().hasData()) {
315 for (
int i = 0; i < nMonomer; ++i) {
316 for (
int j = 0; j < nMonomer; ++j) {
317 double p = interaction_.p(i,j);
319 for (
int k = 0; k < nBasis; ++k) {
320 int idx = i*nBasis + k;
321 resid[idx] += p * h[k];
328 if (!system().mixture().isCanonical()) {
329 for (
int i = 0; i < nMonomer; ++i) {
330 resid[i*nBasis] -= 1.0 / interaction_.sumChiInverse();
334 for (
int i = 0; i < nMonomer; ++i) {
335 resid[i*nBasis] = 0.0;
341 const int nParam = system().domain().unitCell().nParameter();
354 for (
int i = 0; i < nParam ; i++) {
355 if (flexibleParams_[i]) {
356 double str = stress(i);
358 resid[nMonomer*nBasis + counter] = -1 * scaleStress_ * str;
371 UTIL_CHECK(system().domain().basis().isInitialized());
372 const int nMonomer = system().mixture().nMonomer();
373 const int nBasis = system().domain().basis().nBasis();
379 for (
int i = 0; i < nMonomer; i++) {
381 for (
int k = 0; k < nBasis; k++) {
382 wField[i][k] = newGuess[i*nBasis + k];
387 if (system().mixture().isCanonical()) {
389 for (
int i = 0; i < nMonomer; ++i) {
391 for (
int j = 0; j < nMonomer; ++j) {
392 chi = interaction_.chi(i,j);
393 wField[i][0] += chi * system().c().basis(j)[0];
397 if (system().h().hasData()) {
398 for (
int i = 0; i < nMonomer; ++i) {
399 wField[i][0] += system().h().basis(i)[0];
403 system().w().setBasis(wField);
406 const int nParam = system().domain().unitCell().nParameter();
407 const int begin = nMonomer*nBasis;
410 parameters = system().domain().unitCell().parameters();
412 const double coeff = 1.0 / scaleStress_;
414 for (
int i = 0; i < nParam; i++) {
415 if (flexibleParams_[i]) {
416 parameters[i] = coeff * newGuess[begin + counter];
422 system().setUnitCell(parameters);
428 void AmIteratorBasis<D>::outputToLog()
430 if (isFlexible_ && verbose() > 1) {
431 const int nParam = system().domain().unitCell().nParameter();
432 const int nMonomer = system().mixture().nMonomer();
433 const int nBasis = system().domain().basis().nBasis();
435 for (
int i = 0; i < nParam; i++) {
436 if (flexibleParams_[i]) {
437 double str = residual()[nMonomer*nBasis + counter] /
438 (-1.0 * scaleStress_);
440 <<
" Cell Param " << i <<
" = "
441 <<
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 setClassName(const char *className)
Set class name string.
~AmIteratorBasis()
Destructor.
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.
void readParameters(std::istream &in)
Read all parameters and initialize.
void outputTimers(std::ostream &out) const
Output timing results to log file.
void setup(bool isContinuation)
Setup iterator just before entering iteration loop.
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.
Periodic fields and crystallography.
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.
float product(float a, float b)
Product for float Data.