1#ifndef RPC_AM_ITERATOR_BASIS_TPP
2#define RPC_AM_ITERATOR_BASIS_TPP
11#include "AmIteratorBasis.h"
12#include <rpc/System.h>
13#include <pscf/inter/Interaction.h>
14#include <pscf/iterator/NanException.h>
27 imposedFields_(system)
47 interaction_.setNMonomer(system().mixture().nMonomer());
53 int np = system().domain().unitCell().nParameter();
59 readOptional(in,
"isFlexible", isFlexible_);
64 flexibleParams_.clear();
65 for (
int i = 0; i < np; i++) {
66 flexibleParams_.append(
true);
69 readOptionalFSArray(in,
"flexibleParams", flexibleParams_, np);
70 if (nFlexibleParams() == 0) isFlexible_ =
false;
72 flexibleParams_.clear();
73 for (
int i = 0; i < np; i++) {
74 flexibleParams_.append(
false);
79 readOptional(in,
"scaleStress", scaleStress_);
82 readParamCompositeOptional(in, imposedFields_);
91 out <<
"Iterator times contributions:\n";
101 if (imposedFields_.isActive()) {
102 imposedFields_.setup();
106 interaction_.update(system().interaction());
125 for (
int i = 0; i < n; i++) {
127 if (std::isnan(a[i]) || std::isnan(b[i])) {
128 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
143 for (
int i = 0; i < n; i++) {
145 if (std::isnan(value)) {
146 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
149 if (fabs(value) > max) {
165 const int n = hists[0].capacity();
170 for (
int i = 0; i < n; i++) {
171 newbasis[i] = hists[0][i] - hists[1][i];
173 basis.append(newbasis);
185 for (
int i = 0; i < nHist; i++) {
186 for (
int j = 0; j < n; j++) {
188 trial[j] += coeffs[i] * -1 * basis[i][j];
195 void AmIteratorBasis<D>::addPredictedError(
DArray<double>& fieldTrial,
200 for (
int i = 0; i < n; i++) {
201 fieldTrial[i] += lambda * resTrial[i];
209 bool AmIteratorBasis<D>::hasInitialGuess()
210 {
return system().w().hasData(); }
214 int AmIteratorBasis<D>::nElements()
216 const int nMonomer = system().mixture().nMonomer();
217 const int nBasis = system().domain().basis().nBasis();
219 int nEle = nMonomer*nBasis;
221 nEle += nFlexibleParams();
231 const int nMonomer = system().mixture().nMonomer();
232 const int nBasis = system().domain().basis().nBasis();
236 for (
int i = 0; i < nMonomer; i++) {
237 for (
int k = 0; k < nBasis; k++) {
238 curr[i*nBasis+k] = currSys[i][k];
244 const int nParam = system().domain().unitCell().nParameter();
246 = system().domain().unitCell().parameters();
248 for (
int i = 0; i < nParam; i++) {
249 if (flexibleParams_[i]) {
250 curr[nMonomer*nBasis + counter] = scaleStress_*currParam[i];
261 void AmIteratorBasis<D>::evaluate()
265 system().compute(isFlexible_);
272 const int n = nElements();
273 const int nMonomer = system().mixture().nMonomer();
274 const int nBasis = system().domain().basis().nBasis();
277 for (
int i = 0 ; i < n; ++i) {
282 for (
int i = 0; i < nMonomer; ++i) {
283 for (
int j = 0; j < nMonomer; ++j) {
284 double chi = interaction_.chi(i,j);
285 double p = interaction_.p(i,j);
288 for (
int k = 0; k < nBasis; ++k) {
289 int idx = i*nBasis + k;
290 resid[idx] += chi*c[k] - p*w[k];
296 if (system().hasMask()) {
298 double sumChiInv = interaction_.sumChiInverse();
299 for (
int i = 0; i < nMonomer; ++i) {
300 for (
int k = 0; k < nBasis; ++k) {
301 int idx = i*nBasis + k;
302 resid[idx] -= mask[k] / sumChiInv;
309 if (system().hasExternalFields()) {
310 for (
int i = 0; i < nMonomer; ++i) {
311 for (
int j = 0; j < nMonomer; ++j) {
312 double p = interaction_.p(i,j);
314 for (
int k = 0; k < nBasis; ++k) {
315 int idx = i*nBasis + k;
316 resid[idx] += p * h[k];
323 if (!system().mixture().isCanonical()) {
324 if (!system().hasMask()) {
325 for (
int i = 0; i < nMonomer; ++i) {
326 resid[i*nBasis] -= 1.0 / interaction_.sumChiInverse();
331 for (
int i = 0; i < nMonomer; ++i) {
332 resid[i*nBasis] = 0.0;
338 const int nParam = system().domain().unitCell().nParameter();
350 for (
int i = 0; i < nParam ; i++) {
351 if (flexibleParams_[i]) {
352 double stress = system().mixture().stress(i);
355 if (imposedFields_.isActive()) {
356 stress = imposedFields_.correctedStress(i,stress);
359 resid[nMonomer*nBasis + counter] = -1 * scaleStress_ * stress;
373 const int nMonomer = system().mixture().nMonomer();
374 const int nBasis = system().domain().basis().nBasis();
380 for (
int i = 0; i < nMonomer; i++) {
382 for (
int k = 0; k < nBasis; k++)
384 wField[i][k] = newGuess[i*nBasis + k];
388 if (system().mixture().isCanonical()) {
390 for (
int i = 0; i < nMonomer; ++i) {
392 for (
int j = 0; j < nMonomer; ++j) {
393 chi = interaction_.chi(i,j);
394 wField[i][0] += chi * system().c().basis(j)[0];
398 if (system().hasExternalFields()) {
399 for (
int i = 0; i < nMonomer; ++i) {
400 wField[i][0] += system().h().basis(i)[0];
404 system().setWBasis(wField);
407 const int nParam = system().domain().unitCell().nParameter();
408 const int begin = nMonomer*nBasis;
411 parameters = system().domain().unitCell().parameters();
413 double coeff = 1.0 / scaleStress_;
415 for (
int i = 0; i < nParam; i++) {
416 if (flexibleParams_[i]) {
417 parameters[i] = coeff * newGuess[begin + counter];
423 system().setUnitCell(parameters);
427 if (imposedFields_.isActive()) {
428 imposedFields_.update();
434 void AmIteratorBasis<D>::outputToLog()
436 if (isFlexible() && verbose() > 1) {
437 const int nParam = system().domain().unitCell().nParameter();
438 const int nMonomer = system().mixture().nMonomer();
439 const int nBasis = system().domain().basis().nBasis();
441 for (
int i = 0; i < nParam; i++) {
442 if (flexibleParams_[i]) {
443 double stress = residual()[nMonomer*nBasis + counter] /
444 (-1.0 * scaleStress_);
446 <<
" Cell Param " << i <<
" = "
447 <<
Dbl(system().domain().unitCell().parameters()[i], 15)
462 if (imposedFields_.isActive()) {
463 arr = imposedFields_.getParameterTypes();
471 double value,
bool& success)
473 if (imposedFields_.isActive()) {
474 imposedFields_.setParameter(name, ids, value, success);
486 if (imposedFields_.isActive()) {
487 return imposedFields_.getParameter(name, ids, success);
Template for Anderson mixing iterator algorithm.
Exception thrown when not-a-number (NaN) is encountered.
Base template for UnitCell<D> classes, D=1, 2 or 3.
Rpc implementation of the Anderson Mixing iterator with symmetry.
void setup(bool isContinuation)
Setup iterator just before entering iteration loop.
GArray< ParameterType > getParameterTypes()
Return specialized sweep parameter types to add to the Sweep object.
double getParameter(std::string name, DArray< int > ids, bool &success) const
Get the value of a specialized sweep parameter.
void setClassName(const char *className)
Set class name string.
AmIteratorBasis(System< D > &system)
Constructor.
~AmIteratorBasis()
Destructor.
void readParameters(std::istream &in)
Read all parameters and initialize.
void setParameter(std::string name, DArray< int > ids, double value, bool &success)
Set the value of a specialized sweep parameter.
void outputTimers(std::ostream &out)
Output timing results to log file.
Base class for iterative solvers for SCF equations.
Main class for SCFT or PS-FTS simulation of one 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.
An automatically growable array, analogous to a std::vector.
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.
cudaReal max(DeviceArray< cudaReal > const &in)
Get maximum of array elements (GPU kernel wrapper).
PSCF package top-level namespace.
Utility classes for scientific computation.
float product(float a, float b)
Product for float Data.