1#ifndef RPG_AM_ITERATOR_BASIS_TPP
2#define RPG_AM_ITERATOR_BASIS_TPP
11#include "AmIteratorBasis.h"
12#include <rpg/system/System.h>
13#include <prdc/crystal/UnitCell.h>
14#include <prdc/crystal/Basis.h>
15#include <pscf/inter/Interaction.h>
16#include <pscf/iterator/NanException.h>
61 int np = system().domain().unitCell().nParameter();
64 UTIL_CHECK(system().domain().unitCell().lattice()
70 flexibleParams_.clear();
71 for (
int i = 0; i < np; i++) {
72 flexibleParams_.append(
true);
76 if (nFlexibleParams() == 0) isFlexible_ =
false;
78 flexibleParams_.clear();
79 for (
int i = 0; i < np; i++) {
80 flexibleParams_.append(
false);
95 interaction_.setNMonomer(system().mixture().nMonomer());
106 out <<
"Iterator times contributions:\n";
120 interaction_.update(system().interaction());
129 int AmIteratorBasis<D>::nElements()
131 const int nMonomer = system().mixture().nMonomer();
132 const int nBasis = system().domain().basis().nBasis();
134 int nEle = nMonomer*nBasis;
137 nEle += nFlexibleParams();
147 bool AmIteratorBasis<D>::hasInitialGuess()
148 {
return system().w().hasData(); }
158 const int nMonomer = system().mixture().nMonomer();
159 const int nBasis = system().domain().basis().nBasis();
162 for (
int i = 0; i < nMonomer; i++) {
163 for (
int k = 0; k < nBasis; k++) {
164 curr[i*nBasis+k] = currSys[i][k];
169 const int begin = nMonomer*nBasis;
170 const int nParam = system().domain().unitCell().nParameter();
172 = system().domain().unitCell().parameters();
174 for (
int i = 0; i < nParam; i++) {
175 if (flexibleParams_[i]) {
176 curr[begin + counter] = scaleStress_ * parameters[i];
189 void AmIteratorBasis<D>::evaluate()
193 system().compute(isFlexible_);
202 UTIL_CHECK(system().domain().basis().isInitialized());
203 const int nMonomer = system().mixture().nMonomer();
204 const int nBasis = system().domain().basis().nBasis();
205 const int n = nElements();
208 for (
int i = 0 ; i < n; ++i) {
213 for (
int i = 0; i < nMonomer; ++i) {
214 for (
int j = 0; j < nMonomer; ++j) {
215 double chi = interaction_.chi(i,j);
216 double p = interaction_.p(i,j);
219 for (
int k = 0; k < nBasis; ++k) {
220 int idx = i*nBasis + k;
221 resid[idx] += chi*c[k] - p*w[k];
227 if (system().mask().hasData()) {
229 double sumChiInv = interaction_.sumChiInverse();
230 for (
int i = 0; i < nMonomer; ++i) {
231 for (
int k = 0; k < nBasis; ++k) {
232 int idx = i*nBasis + k;
233 resid[idx] -= mask[k] / sumChiInv;
240 if (system().h().hasData()) {
241 for (
int i = 0; i < nMonomer; ++i) {
242 for (
int j = 0; j < nMonomer; ++j) {
243 double p = interaction_.p(i,j);
245 for (
int k = 0; k < nBasis; ++k) {
246 int idx = i*nBasis + k;
247 resid[idx] += p * h[k];
254 if (!system().mixture().isCanonical()) {
255 for (
int i = 0; i < nMonomer; ++i) {
256 resid[i*nBasis] -= 1.0 / interaction_.sumChiInverse();
260 for (
int i = 0; i < nMonomer; ++i) {
261 resid[i*nBasis] = 0.0;
267 const int nParam = system().domain().unitCell().nParameter();
280 for (
int i = 0; i < nParam ; i++) {
281 if (flexibleParams_[i]) {
282 double str = stress(i);
284 resid[nMonomer*nBasis + counter] = -1 * scaleStress_ * str;
299 UTIL_CHECK(system().domain().basis().isInitialized());
300 const int nMonomer = system().mixture().nMonomer();
301 const int nBasis = system().domain().basis().nBasis();
307 for (
int i = 0; i < nMonomer; i++) {
309 for (
int k = 0; k < nBasis; k++) {
310 wField[i][k] = newGuess[i*nBasis + k];
315 if (system().mixture().isCanonical()) {
317 for (
int i = 0; i < nMonomer; ++i) {
319 for (
int j = 0; j < nMonomer; ++j) {
320 chi = interaction_.chi(i,j);
321 wField[i][0] += chi * system().c().basis(j)[0];
325 if (system().h().hasData()) {
326 for (
int i = 0; i < nMonomer; ++i) {
327 wField[i][0] += system().h().basis(i)[0];
331 system().w().setBasis(wField);
334 const int nParam = system().domain().unitCell().nParameter();
335 const int begin = nMonomer*nBasis;
338 parameters = system().domain().unitCell().parameters();
340 const double coeff = 1.0 / scaleStress_;
342 for (
int i = 0; i < nParam; i++) {
343 if (flexibleParams_[i]) {
344 parameters[i] = coeff * newGuess[begin + counter];
350 system().setUnitCell(parameters);
358 void AmIteratorBasis<D>::outputToLog()
360 if (isFlexible_ && verbose() > 1) {
361 const int nParam = system().domain().unitCell().nParameter();
362 const int nMonomer = system().mixture().nMonomer();
363 const int nBasis = system().domain().basis().nBasis();
365 for (
int i = 0; i < nParam; i++) {
366 if (flexibleParams_[i]) {
367 double str = residual()[nMonomer*nBasis + counter] /
368 (-1.0 * scaleStress_);
370 <<
" Cell Param " << i <<
" = "
371 <<
Dbl(system().domain().unitCell().parameters()[i], 15)
402 for (
int i=0; i < n; ++i) {
404 if (std::isnan(a[i]) || std::isnan(b[i])) {
405 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
422 for (
int i = 0; i < n; i++) {
424 if (std::isnan(value)) {
425 throw NanException(
"AmIteratorBasis::dotProduct", __FILE__,
428 if (fabs(value) > max)
445 for (
int i = 0; i < n; i++) {
460 for (
int i = 0; i < n; i++) {
virtual void setup(bool isContinuation)
void readMixingParameters(std::istream &in, bool useLambdaRamp=true)
void outputTimers(std::ostream &out) const
void readErrorType(std::istream &in)
virtual void readParameters(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.
~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 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.
void setClassName(const char *className)
Set class name string.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
double max(Array< double > const &in)
Get maximum of array elements .
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.