1#ifndef RPG_AM_ITERATOR_BASIS_TPP
2#define RPG_AM_ITERATOR_BASIS_TPP
11#include "AmIteratorBasis.h"
12#include <rpg/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>
33 imposedFields_(system)
53 interaction_.setNMonomer(system().mixture().nMonomer());
59 int np = system().domain().unitCell().nParameter();
65 readOptional(in,
"isFlexible", isFlexible_);
70 flexibleParams_.clear();
71 for (
int i = 0; i < np; i++) {
72 flexibleParams_.append(
true);
75 readOptionalFSArray(in,
"flexibleParams", flexibleParams_, np);
76 if (nFlexibleParams() == 0) isFlexible_ =
false;
78 flexibleParams_.clear();
79 for (
int i = 0; i < np; i++) {
80 flexibleParams_.append(
false);
85 readOptional(in,
"scaleStress", scaleStress_);
88 readParamCompositeOptional(in, imposedFields_);
97 out <<
"Iterator times contributions:\n";
107 if (imposedFields_.isActive()) {
108 imposedFields_.setup();
115 interaction_.update(system().interaction());
135 for (
int i=0; i < n; ++i) {
137 if (std::isnan(a[i]) || std::isnan(b[i])) {
138 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
153 for (
int i = 0; i < n; i++) {
155 if (std::isnan(value)) {
156 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
159 if (fabs(value) > max)
174 const int n = hists[0].capacity();
179 for (
int i = 0; i < n; i++) {
180 newbasis[i] = hists[0][i] - hists[1][i];
183 basis.append(newbasis);
195 for (
int i = 0; i < nHist; i++) {
196 for (
int j = 0; j < n; j++) {
198 trial[j] += coeffs[i] * -1 * basis[i][j];
205 void AmIteratorBasis<D>::addPredictedError(
DArray<double>& fieldTrial,
210 for (
int i = 0; i < n; i++) {
211 fieldTrial[i] += lambda * resTrial[i];
217 bool AmIteratorBasis<D>::hasInitialGuess()
219 return system().w().hasData();
224 int AmIteratorBasis<D>::nElements()
226 const int nMonomer = system().mixture().nMonomer();
227 const int nBasis = system().basis().nBasis();
229 int nEle = nMonomer*nBasis;
232 nEle += nFlexibleParams();
244 const int nMonomer = system().mixture().nMonomer();
245 const int nBasis = system().basis().nBasis();
248 for (
int i = 0; i < nMonomer; i++) {
249 for (
int k = 0; k < nBasis; k++) {
250 curr[i*nBasis+k] = currSys[i][k];
255 const int begin = nMonomer*nBasis;
256 const int nParam = system().unitCell().nParameter();
258 = system().unitCell().parameters();
260 for (
int i = 0; i < nParam; i++) {
261 if (flexibleParams_[i]) {
262 curr[begin + counter] = scaleStress_ * parameters[i];
273 void AmIteratorBasis<D>::evaluate()
277 system().compute(isFlexible_);
285 const int nMonomer = system().mixture().nMonomer();
286 const int nBasis = system().basis().nBasis();
287 const int n = nElements();
290 for (
int i = 0 ; i < n; ++i) {
295 for (
int i = 0; i < nMonomer; ++i) {
296 for (
int j = 0; j < nMonomer; ++j) {
297 double chi = interaction_.chi(i,j);
298 double p = interaction_.p(i,j);
301 for (
int k = 0; k < nBasis; ++k) {
302 int idx = i*nBasis + k;
303 resid[idx] += chi*c[k] - p*w[k];
309 if (system().hasMask()) {
311 double sumChiInv = interaction_.sumChiInverse();
312 for (
int i = 0; i < nMonomer; ++i) {
313 for (
int k = 0; k < nBasis; ++k) {
314 int idx = i*nBasis + k;
315 resid[idx] -= mask[k] / sumChiInv;
322 if (system().hasExternalFields()) {
323 for (
int i = 0; i < nMonomer; ++i) {
324 for (
int j = 0; j < nMonomer; ++j) {
325 double p = interaction_.p(i,j);
327 for (
int k = 0; k < nBasis; ++k) {
328 int idx = i*nBasis + k;
329 resid[idx] += p * h[k];
336 if (!system().mixture().isCanonical()) {
337 for (
int i = 0; i < nMonomer; ++i) {
338 resid[i*nBasis] -= 1.0 / interaction_.sumChiInverse();
342 for (
int i = 0; i < nMonomer; ++i) {
343 resid[i*nBasis] = 0.0;
349 const int nParam = system().unitCell().nParameter();
362 for (
int i = 0; i < nParam ; i++) {
363 if (flexibleParams_[i]) {
364 double stress = system().mixture().stress(i);
367 if (imposedFields_.isActive()) {
368 stress = imposedFields_.correctedStress(i,stress);
371 resid[nMonomer*nBasis + counter] = -1 * scaleStress_ * stress;
385 const int nMonomer = system().mixture().nMonomer();
386 const int nBasis = system().basis().nBasis();
392 for (
int i = 0; i < nMonomer; i++) {
394 for (
int k = 0; k < nBasis; k++) {
395 wField[i][k] = newGuess[i*nBasis + k];
400 if (system().mixture().isCanonical()) {
402 for (
int i = 0; i < nMonomer; ++i) {
404 for (
int j = 0; j < nMonomer; ++j) {
405 chi = interaction_.chi(i,j);
406 wField[i][0] += chi * system().c().basis(j)[0];
410 if (system().hasExternalFields()) {
411 for (
int i = 0; i < nMonomer; ++i) {
412 wField[i][0] += system().h().basis(i)[0];
416 system().setWBasis(wField);
419 const int nParam = system().unitCell().nParameter();
420 const int begin = nMonomer*nBasis;
423 parameters = system().domain().unitCell().parameters();
425 const double coeff = 1.0 / scaleStress_;
427 for (
int i = 0; i < nParam; i++) {
428 if (flexibleParams_[i]) {
429 parameters[i] = coeff * newGuess[begin + counter];
435 system().setUnitCell(parameters);
439 if (imposedFields_.isActive()) {
440 imposedFields_.update();
446 void AmIteratorBasis<D>::outputToLog()
448 if (isFlexible_ && verbose() > 1) {
449 const int nParam = system().domain().unitCell().nParameter();
450 const int nMonomer = system().mixture().nMonomer();
451 const int nBasis = system().domain().basis().nBasis();
453 for (
int i = 0; i < nParam; i++) {
454 if (flexibleParams_[i]) {
455 double stress = residual()[nMonomer*nBasis + counter] /
456 (-1.0 * scaleStress_);
458 <<
" Cell Param " << i <<
" = "
459 <<
Dbl(system().domain().unitCell().parameters()[i], 15)
474 if (imposedFields_.isActive()) {
475 arr = imposedFields_.getParameterTypes();
483 double value,
bool& success)
485 if (imposedFields_.isActive()) {
486 imposedFields_.setParameter(name, ids, value, success);
498 if (imposedFields_.isActive()) {
499 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.
Rpg implementation of the Anderson Mixing iterator.
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()
Destructor.
AmIteratorBasis(System< D > &system)
Constructor.
void readParameters(std::istream &in)
Read all parameters and initialize.
void setup(bool isContinuation)
Setup iterator just before entering iteration loop.
void outputTimers(std::ostream &out)
Output timing results to log file.
void setParameter(std::string name, DArray< int > ids, double value, bool &success)
Set the value of a specialized sweep parameter.
GArray< ParameterType > getParameterTypes()
Return specialized sweep parameter types to add to the Sweep object.
Base class for iterative solvers for SCF equations.
Main class for calculations that represent 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.