PSCF v1.3
rpg/scft/iterator/AmIteratorBasis.tpp
1#ifndef RPG_AM_ITERATOR_BASIS_TPP
2#define RPG_AM_ITERATOR_BASIS_TPP
3
4/*
5* PSCF - Polymer Self-Consistent Field
6*
7* Copyright 2015 - 2025, The Regents of the University of Minnesota
8* Distributed under the terms of the GNU General Public License.
9*/
10
11#include "AmIteratorBasis.h"
12#include <rpg/system/System.h>
13
14#include <prdc/crystal/UnitCell.h>
15#include <prdc/crystal/Basis.h>
16
17#include <pscf/inter/Interaction.h>
18#include <pscf/iterator/NanException.h>
19
20#include <util/global.h>
21#include <cmath>
22
23namespace Pscf {
24namespace Rpg {
25
26 using namespace Util;
27 using namespace Prdc;
28
29 // Constructor
30 template <int D>
32 : Iterator<D>(system)
33 {
34 isSymmetric_ = true;
35 setClassName("AmIteratorBasis");
36 }
37
38 // Destructor
39 template <int D>
42
43 // Read parameter file block
44 template <int D>
46 {
47 // Call parent class readParameters
50
51 // Allocate local modified copy of Interaction class
52 interaction_.setNMonomer(system().mixture().nMonomer());
53
54 // Default parameter values
55 isFlexible_ = 1;
56 scaleStress_ = 10.0;
57
58 int np = system().domain().unitCell().nParameter();
59 UTIL_CHECK(np > 0);
60 UTIL_CHECK(np <= 6);
61 UTIL_CHECK(system().domain().unitCell().lattice() != UnitCell<D>::Null);
62
63 // Read in optional isFlexible value
64 readOptional(in, "isFlexible", isFlexible_);
65
66 // Populate flexibleParams_ based on isFlexible_ (all 0s or all 1s),
67 // then optionally overwrite with user input from param file
68 if (isFlexible_) {
69 flexibleParams_.clear();
70 for (int i = 0; i < np; i++) {
71 flexibleParams_.append(true); // Set all values to true
72 }
73 // Read optional flexibleParams_ array to overwrite current array
74 readOptionalFSArray(in, "flexibleParams", flexibleParams_, np);
75 if (nFlexibleParams() == 0) isFlexible_ = false;
76 } else { // isFlexible_ = false
77 flexibleParams_.clear();
78 for (int i = 0; i < np; i++) {
79 flexibleParams_.append(false); // Set all values to false
80 }
81 }
82
83 // Read optional scaleStress value
84 readOptional(in, "scaleStress", scaleStress_);
85 }
86
87 // Output timing results to log file.
88 template<int D>
89 void AmIteratorBasis<D>::outputTimers(std::ostream& out) const
90 {
91 // Output timing results, if requested.
92 out << "\n";
93 out << "Iterator times contributions:\n";
95 }
96
97 // -- Protected virtual function -- //
98
99 // Setup before entering iteration loop
100 template <int D>
101 void AmIteratorBasis<D>::setup(bool isContinuation)
102 {
103 // Call parent setup method
105
106 // Update chi matrix and related properties in member interaction_
107 interaction_.update(system().interaction());
108 }
109
110
111 // -- Private virtual functions used to implement AM algorithm -- //
112
113 // Set a vector equal to another (assign a = b)
114 template <int D>
115 void AmIteratorBasis<D>::setEqual(DArray<double>& a,
116 DArray<double> const & b)
117 { a = b; }
118
119 // Compute the inner product of two real vectors.
120 template <int D>
121 double AmIteratorBasis<D>::dotProduct(DArray<double> const & a,
122 DArray<double> const & b)
123 {
124 const int n = a.capacity();
125 UTIL_CHECK(b.capacity() == n);
126 double product = 0.0;
127 for (int i=0; i < n; ++i) {
128 // if either value is NaN, throw NanException
129 if (std::isnan(a[i]) || std::isnan(b[i])) {
130 throw NanException("AmIteratorBasis::dotProduct", __FILE__,
131 __LINE__, 0);
132 }
133 product += a[i] * b[i];
134 }
135 return product;
136 }
137
138 // Compute and return maximum element of residual vector.
139 template <int D>
140 double AmIteratorBasis<D>::maxAbs(DArray<double> const & a)
141 {
142 const int n = a.capacity();
143 double max = 0.0;
144 double value;
145 for (int i = 0; i < n; i++) {
146 value = a[i];
147 if (std::isnan(value)) { // if value is NaN, throw NanException
148 throw NanException("AmIteratorBasis::dotProduct", __FILE__,
149 __LINE__, 0);
150 }
151 if (fabs(value) > max)
152 max = fabs(value);
153 }
154 return max;
155 }
156
157 // Update the series of residual vectors.
158 template <int D>
159 void
160 AmIteratorBasis<D>::updateBasis(RingBuffer<DArray<double> > & basis,
161 RingBuffer<DArray<double> > const& hists)
162 {
163 // Make sure at least two histories are stored
164 UTIL_CHECK(hists.size() >= 2);
165
166 const int n = hists[0].capacity();
167 DArray<double> newbasis;
168 newbasis.allocate(n);
169
170 // New basis vector is difference between two most recent states
171 for (int i = 0; i < n; i++) {
172 newbasis[i] = hists[0][i] - hists[1][i];
173 }
174
175 basis.append(newbasis);
176 }
177
178 // Compute trial field so as to minimize L2 norm of residual.
179 template <int D>
180 void
181 AmIteratorBasis<D>::addHistories(DArray<double>& trial,
182 RingBuffer<DArray<double> > const & basis,
183 DArray<double> coeffs,
184 int nHist)
185 {
186 int n = trial.capacity();
187 for (int i = 0; i < nHist; i++) {
188 for (int j = 0; j < n; j++) {
189 // Not clear on the origin of the -1 factor
190 trial[j] += coeffs[i] * -1 * basis[i][j];
191 }
192 }
193 }
194
195 // Add predicted error to the trial field.
196 template <int D>
197 void AmIteratorBasis<D>::addPredictedError(DArray<double>& fieldTrial,
198 DArray<double> const & resTrial,
199 double lambda)
200 {
201 int n = fieldTrial.capacity();
202 for (int i = 0; i < n; i++) {
203 fieldTrial[i] += lambda * resTrial[i];
204 }
205 }
206
207 // Does the system have an initial field guess?
208 template <int D>
209 bool AmIteratorBasis<D>::hasInitialGuess()
210 {
211 return system().w().hasData();
212 }
213
214 // Compute the number of elements in the residual vector.
215 template <int D>
216 int AmIteratorBasis<D>::nElements()
217 {
218 const int nMonomer = system().mixture().nMonomer();
219 const int nBasis = system().domain().basis().nBasis();
220
221 int nEle = nMonomer*nBasis;
222
223 if (isFlexible_) {
224 nEle += nFlexibleParams();
225 }
226
227 return nEle;
228 }
229
230 // Get the current w fields and lattice parameters
231 template <int D>
232 void AmIteratorBasis<D>::getCurrent(DArray<double>& curr)
233 {
234 // Straighten out fields into linear arrays
235
236 const int nMonomer = system().mixture().nMonomer();
237 const int nBasis = system().domain().basis().nBasis();
238 const DArray< DArray<double> > & currSys = system().w().basis();
239
240 for (int i = 0; i < nMonomer; i++) {
241 for (int k = 0; k < nBasis; k++) {
242 curr[i*nBasis+k] = currSys[i][k];
243 }
244 }
245
246 if (isFlexible_) {
247 const int begin = nMonomer*nBasis;
248 const int nParam = system().domain().unitCell().nParameter();
249 FSArray<double,6> const & parameters
250 = system().domain().unitCell().parameters();
251 int counter = 0;
252 for (int i = 0; i < nParam; i++) {
253 if (flexibleParams_[i]) {
254 curr[begin + counter] = scaleStress_ * parameters[i];
255 counter++;
256 }
257 }
258 UTIL_CHECK(counter == nFlexibleParams());
259 }
260
261 }
262
263 // Perform the main system computation (solve the MDE)
264 template <int D>
265 void AmIteratorBasis<D>::evaluate()
266 {
267 // Solve MDEs for current omega field
268 // (computes stress if isFlexible_ == true)
269 system().compute(isFlexible_);
270 }
271
272 // Compute the residual for the current system state
273 template <int D>
274 void AmIteratorBasis<D>::getResidual(DArray<double>& resid)
275 {
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();
280
281 // Initialize residuals
282 for (int i = 0 ; i < n; ++i) {
283 resid[i] = 0.0;
284 }
285
286 // Compute SCF residual vector elements
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);
291 DArray<double> const & c = system().c().basis(j);
292 DArray<double> const & w = system().w().basis(j);
293 for (int k = 0; k < nBasis; ++k) {
294 int idx = i*nBasis + k;
295 resid[idx] += chi*c[k] - p*w[k];
296 }
297 }
298 }
299
300 // If iterator has mask, account for it in residual values
301 if (system().mask().hasData()) {
302 DArray<double> const & mask = system().mask().basis();
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;
308 }
309 }
310 }
311
312 // If iterator has external fields, account for them in the values
313 // of the residuals
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);
318 DArray<double> const & h = system().h().basis(j);
319 for (int k = 0; k < nBasis; ++k) {
320 int idx = i*nBasis + k;
321 resid[idx] += p * h[k];
322 }
323 }
324 }
325 }
326
327 // If not canonical, account for incompressibility
328 if (!system().mixture().isCanonical()) {
329 for (int i = 0; i < nMonomer; ++i) {
330 resid[i*nBasis] -= 1.0 / interaction_.sumChiInverse();
331 }
332 } else {
333 // Explicitly set homogeneous residual components
334 for (int i = 0; i < nMonomer; ++i) {
335 resid[i*nBasis] = 0.0;
336 }
337 }
338
339 // If variable unit cell, compute stress residuals
340 if (isFlexible_) {
341 const int nParam = system().domain().unitCell().nParameter();
342
343 // Note:
344 // Combined -1 factor and stress scaling here. This is okay:
345 // - Residuals only show up as dot products (U, v, norm)
346 // or with their absolute value taken (max), so the
347 // sign on a given residual vector element is not relevant
348 // as long as it is consistent across all vectors
349 // - The scaling is applied here and to the unit cell param
350 // storage, so that updating is done on the same scale,
351 // and then undone right before passing to the unit cell.
352
353 int counter = 0;
354 for (int i = 0; i < nParam ; i++) {
355 if (flexibleParams_[i]) {
356 double str = stress(i);
357
358 resid[nMonomer*nBasis + counter] = -1 * scaleStress_ * str;
359 counter++;
360 }
361 }
362 UTIL_CHECK(counter == nFlexibleParams());
363 }
364
365 }
366
367 // Update the current system field coordinates
368 template <int D>
369 void AmIteratorBasis<D>::update(DArray<double>& newGuess)
370 {
371 UTIL_CHECK(system().domain().basis().isInitialized());
372 const int nMonomer = system().mixture().nMonomer();
373 const int nBasis = system().domain().basis().nBasis();
374
375 DArray< DArray<double> > wField;
376 wField.allocate(nMonomer);
377
378 // Restructure in format of monomers, basis functions
379 for (int i = 0; i < nMonomer; i++) {
380 wField[i].allocate(nBasis);
381 for (int k = 0; k < nBasis; k++) {
382 wField[i][k] = newGuess[i*nBasis + k];
383 }
384 }
385
386 // If canonical, explicitly set homogeneous field components
387 if (system().mixture().isCanonical()) {
388 double chi;
389 for (int i = 0; i < nMonomer; ++i) {
390 wField[i][0] = 0.0; // initialize to 0
391 for (int j = 0; j < nMonomer; ++j) {
392 chi = interaction_.chi(i,j);
393 wField[i][0] += chi * system().c().basis(j)[0];
394 }
395 }
396 // If iterator has external fields, include them in homogeneous field
397 if (system().h().hasData()) {
398 for (int i = 0; i < nMonomer; ++i) {
399 wField[i][0] += system().h().basis(i)[0];
400 }
401 }
402 }
403 system().w().setBasis(wField);
404
405 if (isFlexible_) {
406 const int nParam = system().domain().unitCell().nParameter();
407 const int begin = nMonomer*nBasis;
408
409 FSArray<double,6> parameters;
410 parameters = system().domain().unitCell().parameters();
411
412 const double coeff = 1.0 / scaleStress_;
413 int counter = 0;
414 for (int i = 0; i < nParam; i++) {
415 if (flexibleParams_[i]) {
416 parameters[i] = coeff * newGuess[begin + counter];
417 counter++;
418 }
419 }
420 UTIL_CHECK(counter == nFlexibleParams());
421
422 system().setUnitCell(parameters);
423 }
424 }
425
426 // Output relevant system details to the iteration log file.
427 template<int D>
428 void AmIteratorBasis<D>::outputToLog()
429 {
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();
434 int counter = 0;
435 for (int i = 0; i < nParam; i++) {
436 if (flexibleParams_[i]) {
437 double str = residual()[nMonomer*nBasis + counter] /
438 (-1.0 * scaleStress_);
439 Log::file()
440 << " Cell Param " << i << " = "
441 << Dbl(system().domain().unitCell().parameters()[i], 15)
442 << " , stress = "
443 << Dbl(str, 15)
444 << "\n";
445 counter++;
446 }
447 }
448 }
449 }
450
451}
452}
453#endif
Exception thrown when not-a-number (NaN) is encountered.
Base template for UnitCell<D> classes, D=1, 2 or 3.
Definition UnitCell.h:56
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.
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.
Definition Array.h:159
Dynamically allocatable contiguous array template.
Definition DArray.h:32
void allocate(int capacity)
Allocate the underlying C array.
Definition DArray.h:199
Wrapper for a double precision number, for formatted ostream output.
Definition Dbl.h:40
A fixed capacity (static) contiguous array with a variable logical size.
Definition FSArray.h:38
static std::ostream & file()
Get log ostream by reference.
Definition Log.cpp:59
Class for storing history of previous values in an array.
Definition RingBuffer.h:27
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
Periodic fields and crystallography.
Definition CField.cpp:11
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.
Definition param_pc.dox:1
float product(float a, float b)
Product for float Data.
Definition product.h:22