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 <rpc/solvers/Mixture.h> 
   14#include <rpc/field/Domain.h> 
   15#include <prdc/crystal/Basis.h> 
   16#include <prdc/crystal/UnitCell.h> 
   17#include <pscf/inter/Interaction.h> 
   18#include <pscf/iterator/NanException.h> 
   19#include <util/containers/DArray.h> 
   59      UnitCell<D> const & unitCell = system().domain().unitCell();
 
   72         flexibleParams_.clear();
 
   73         for (
int i = 0; i < np; i++) {
 
   74            flexibleParams_.append(
true); 
 
   78         if (nFlexibleParams() == 0) isFlexible_ = 
false;
 
   80         flexibleParams_.clear();
 
   81         for (
int i = 0; i < np; i++) {
 
   82            flexibleParams_.append(
false); 
 
   94      const int nMonomer = system().mixture().nMonomer();
 
   95      interaction_.setNMonomer(nMonomer);
 
 
  106      out << 
"Iterator times contributions:\n";
 
 
  117      interaction_.update(system().interaction());
 
 
  126   int AmIteratorBasis<D>::nElements()
 
  128      const int nMonomer = system().mixture().nMonomer();
 
  129      const int nBasis = system().domain().basis().nBasis();
 
  131      int nEle = nMonomer*nBasis;
 
  133         nEle += nFlexibleParams();
 
  142   bool AmIteratorBasis<D>::hasInitialGuess()
 
  143   {  
return system().w().hasData(); }
 
  151      const int nMonomer = system().mixture().nMonomer();
 
  152      const int nBasis = system().domain().basis().nBasis();
 
  153      const int nEle = nElements();
 
  158      for (
int i = 0; i < nMonomer; i++) {
 
  161         for (
int k = 0; k < nBasis; k++) {
 
  162            state[begin + k] = field[k];
 
  169         UnitCell<D> const & unitCell = system().domain().unitCell();
 
  171         const int nParam = unitCell.nParameter();
 
  172         begin = nMonomer*nBasis;
 
  174         for (
int i = 0; i < nParam; i++) {
 
  175            if (flexibleParams_[i]) {
 
  176               state[begin + counter] = scaleStress_*parameters[i];
 
  189   void AmIteratorBasis<D>::evaluate()
 
  190   {  system().compute(isFlexible_); }
 
  198      const int nMonomer = system().mixture().nMonomer();
 
  199      const int nBasis = system().domain().basis().nBasis();
 
  200      const int n = nElements();
 
  208      for (i = 0 ; i < n; ++i) {
 
  213      for (i = 0; i < nMonomer; ++i) {
 
  215         for (j = 0; j < nMonomer; ++j) {
 
  216            chi = interaction_.chi(i,j);
 
  217            p = interaction_.p(i,j);
 
  220            if (system().h().hasData()) {
 
  222               for (k = 0; k < nBasis; ++k) {
 
  223                  resid[begin + k] += chi*c[k] + p*(h[k] - w[k]);
 
  226               for (k = 0; k < nBasis; ++k) {
 
  227                  resid[begin + k] += chi*c[k] - p*w[k];
 
  234      double shift = -1.0 / interaction_.sumChiInverse();
 
  235      if (system().mask().hasData()) {
 
  237         for (i = 0; i < nMonomer; ++i) {
 
  239            for (k = 0; k < nBasis; ++k) {
 
  240               resid[begin + k] += shift*mask[k];
 
  244         for (i = 0; i < nMonomer; ++i) {
 
  245            resid[i*nBasis] += shift;
 
  250      if (system().mixture().isCanonical()) {
 
  251         for (i = 0; i < nMonomer; ++i) {
 
  252            resid[i*nBasis] = 0.0;
 
  268         double coeff = -1.0 * scaleStress_;
 
  269         const int nParam = system().domain().unitCell().nParameter();
 
  270         begin = nMonomer*nBasis;
 
  272         for (i = 0; i < nParam ; i++) {
 
  273            if (flexibleParams_[i]) {
 
  274               resid[begin + counter] = coeff * stress(i);
 
  290      const int nMonomer = system().mixture().nMonomer();
 
  291      const int nBasis = system().domain().basis().nBasis();
 
  296      for (
int i = 0; i < nMonomer; i++) {
 
  302      for (
int i = 0; i < nMonomer; i++) {
 
  304         for (
int k = 0; k < nBasis; k++) {
 
  305            wFields[i][k] = newState[begin + k];
 
  310      if (system().mixture().isCanonical()) {
 
  313         for (
int i = 0; i < nMonomer; ++i) {
 
  318         double chi, wAve, cAve;
 
  319         for (
int i = 0; i < nMonomer; ++i) {
 
  321            for (
int j = 0; j < nMonomer; ++j) {
 
  322               chi = interaction_.chi(i,j);
 
  323               cAve = system().c().basis(j)[0];
 
  326            wFields[i][0] = wAve;
 
  330         if (system().h().hasData()) {
 
  331            for (
int i = 0; i < nMonomer; ++i) {
 
  332               wFields[i][0] += system().h().basis(i)[0];
 
  338      system().w().setBasis(wFields);
 
  345         parameters = system().domain().unitCell().parameters();
 
  348         const double coeff = 1.0 / scaleStress_;
 
  349         const int nParam = system().domain().unitCell().nParameter();
 
  350         const int begin = nMonomer*nBasis;
 
  352         for (
int i = 0; i < nParam; i++) {
 
  353            if (flexibleParams_[i]) {
 
  354               parameters[i] = coeff * newState[begin + counter];
 
  361         system().setUnitCell(parameters);
 
  369   void AmIteratorBasis<D>::outputToLog()
 
  371      if (isFlexible() && verbose() > 1) {
 
  373         UnitCell<D> const & unitCell = system().domain().unitCell();
 
  375         const int nMonomer = system().mixture().nMonomer();
 
  376         const int nBasis = system().domain().basis().nBasis();
 
  377         const int begin = nMonomer*nBasis;
 
  379         for (
int i = 0; i < nParam; i++) {
 
  380            if (flexibleParams_[i]) {
 
  381               str = -1.0 * residual()[begin + counter] / scaleStress_;
 
  383                  << 
" Cell Param  " << i << 
" = " 
  384                  << 
Dbl(unitCell.parameters()[i], 15)
 
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)
int nParameter() const
Get the number of parameters in the unit cell.
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 readParameters(std::istream &in) override
Read all parameters and initialize.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
void outputTimers(std::ostream &out) const override
Output timing results to log file.
AmIteratorBasis(System< D > &system)
Constructor.
~AmIteratorBasis()
Destructor.
void setup(bool isContinuation) override
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.
Periodic fields and crystallography.
Real periodic fields, SCFT and PS-FTS (CPU).
PSCF package top-level namespace.