PSCF v1.3.2
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#include <prdc/crystal/UnitCell.h>
14#include <prdc/crystal/Basis.h>
15#include <pscf/inter/Interaction.h>
16#include <pscf/iterator/NanException.h>
17#include <util/global.h>
18#include <cmath>
19
20namespace Pscf {
21namespace Rpg {
22
23 using namespace Util;
24 using namespace Prdc;
25
26 // Public member functions
27
28 /*
29 * Constructor.
30 */
31 template <int D>
33 : Iterator<D>(system)
34 {
35 isSymmetric_ = true;
36 ParamComposite::setClassName("AmIteratorBasis");
37 }
38
39 /*
40 * Destructor.
41 */
42 template <int D>
45
46 /*
47 * Read parameter file block.
48 */
49 template <int D>
51 {
52 // Call parent class readParameters
55
56 // Read optional isFlexible value
57 isFlexible_ = 1; // default value
58 readOptional(in, "isFlexible", isFlexible_);
59
60 // Get and check the number of unit cell parameters
61 int np = system().domain().unitCell().nParameter();
62 UTIL_CHECK(np > 0);
63 UTIL_CHECK(np <= 6);
64 UTIL_CHECK(system().domain().unitCell().lattice()
66
67 // Populate flexibleParams_ based on isFlexible_ (all 0s or all 1s),
68 // then optionally overwrite with user input from param file
69 if (isFlexible_) {
70 flexibleParams_.clear();
71 for (int i = 0; i < np; i++) {
72 flexibleParams_.append(true); // Set all values to true
73 }
74 // Read optional flexibleParams_ array to overwrite current array
75 readOptionalFSArray(in, "flexibleParams", flexibleParams_, np);
76 if (nFlexibleParams() == 0) isFlexible_ = false;
77 } else { // isFlexible_ = false
78 flexibleParams_.clear();
79 for (int i = 0; i < np; i++) {
80 flexibleParams_.append(false); // Set all values to false
81 }
82 }
83
84 // Read optional scaleStress value
85 scaleStress_ = 10.0; // default value
86 readOptional(in, "scaleStress", scaleStress_);
87
88 // Optionally read mixing parameters (lambda, useLambdaRamp, r)
90
91 // Optionally read mixing parameters (lambda, useLambdaRamp, r)
93
94 // Allocate local modified copy of Interaction class
95 interaction_.setNMonomer(system().mixture().nMonomer());
96 }
97
98 /*
99 * Output timing results to log file.
100 */
101 template<int D>
102 void AmIteratorBasis<D>::outputTimers(std::ostream& out) const
103 {
104 // Output timing results, if requested.
105 out << "\n";
106 out << "Iterator times contributions:\n";
108 }
109
110 // Protected virtual function
111
112 // Setup before entering iteration loop
113 template <int D>
114 void AmIteratorBasis<D>::setup(bool isContinuation)
115 {
116 // Call parent setup method
117 AmTmpl::setup(isContinuation);
118
119 // Update chi matrix and related properties in member interaction_
120 interaction_.update(system().interaction());
121 }
122
123 // Private virtual functions that interact with parent system
124
125 /*
126 * Compute the number of elements in the residual vector.
127 */
128 template <int D>
129 int AmIteratorBasis<D>::nElements()
130 {
131 const int nMonomer = system().mixture().nMonomer();
132 const int nBasis = system().domain().basis().nBasis();
133
134 int nEle = nMonomer*nBasis;
135
136 if (isFlexible_) {
137 nEle += nFlexibleParams();
138 }
139
140 return nEle;
141 }
142
143 /*
144 * Does the system have an initial field guess?
145 */
146 template <int D>
147 bool AmIteratorBasis<D>::hasInitialGuess()
148 { return system().w().hasData(); }
149
150 /*
151 * Get the current w fields and lattice parameters.
152 */
153 template <int D>
154 void AmIteratorBasis<D>::getCurrent(DArray<double>& curr)
155 {
156 // Straighten out fields into linear arrays
157
158 const int nMonomer = system().mixture().nMonomer();
159 const int nBasis = system().domain().basis().nBasis();
160 const DArray< DArray<double> > & currSys = system().w().basis();
161
162 for (int i = 0; i < nMonomer; i++) {
163 for (int k = 0; k < nBasis; k++) {
164 curr[i*nBasis+k] = currSys[i][k];
165 }
166 }
167
168 if (isFlexible_) {
169 const int begin = nMonomer*nBasis;
170 const int nParam = system().domain().unitCell().nParameter();
171 FSArray<double,6> const & parameters
172 = system().domain().unitCell().parameters();
173 int counter = 0;
174 for (int i = 0; i < nParam; i++) {
175 if (flexibleParams_[i]) {
176 curr[begin + counter] = scaleStress_ * parameters[i];
177 counter++;
178 }
179 }
180 UTIL_CHECK(counter == nFlexibleParams());
181 }
182
183 }
184
185 /*
186 * Perform the main system computation (solve the MDE).
187 */
188 template <int D>
189 void AmIteratorBasis<D>::evaluate()
190 {
191 // Solve MDEs for current omega field
192 // (computes stress if isFlexible_ == true)
193 system().compute(isFlexible_);
194 }
195
196 /*
197 * Compute the residual for the current system state.
198 */
199 template <int D>
200 void AmIteratorBasis<D>::getResidual(DArray<double>& resid)
201 {
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();
206
207 // Initialize residuals
208 for (int i = 0 ; i < n; ++i) {
209 resid[i] = 0.0;
210 }
211
212 // Compute SCF residual vector elements
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);
217 DArray<double> const & c = system().c().basis(j);
218 DArray<double> const & w = system().w().basis(j);
219 for (int k = 0; k < nBasis; ++k) {
220 int idx = i*nBasis + k;
221 resid[idx] += chi*c[k] - p*w[k];
222 }
223 }
224 }
225
226 // If iterator has mask, account for it in residual values
227 if (system().mask().hasData()) {
228 DArray<double> const & mask = system().mask().basis();
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;
234 }
235 }
236 }
237
238 // If iterator has external fields, account for them in the values
239 // of the residuals
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);
244 DArray<double> const & h = system().h().basis(j);
245 for (int k = 0; k < nBasis; ++k) {
246 int idx = i*nBasis + k;
247 resid[idx] += p * h[k];
248 }
249 }
250 }
251 }
252
253 // If not canonical, account for incompressibility
254 if (!system().mixture().isCanonical()) {
255 for (int i = 0; i < nMonomer; ++i) {
256 resid[i*nBasis] -= 1.0 / interaction_.sumChiInverse();
257 }
258 } else {
259 // Explicitly set homogeneous residual components
260 for (int i = 0; i < nMonomer; ++i) {
261 resid[i*nBasis] = 0.0;
262 }
263 }
264
265 // If variable unit cell, compute stress residuals
266 if (isFlexible_) {
267 const int nParam = system().domain().unitCell().nParameter();
268
269 // Note:
270 // Combined -1 factor and stress scaling here. This is okay:
271 // - Residuals only show up as dot products (U, v, norm)
272 // or with their absolute value taken (max), so the
273 // sign on a given residual vector element is not relevant
274 // as long as it is consistent across all vectors
275 // - The scaling is applied here and to the unit cell param
276 // storage, so that updating is done on the same scale,
277 // and then undone right before passing to the unit cell.
278
279 int counter = 0;
280 for (int i = 0; i < nParam ; i++) {
281 if (flexibleParams_[i]) {
282 double str = stress(i);
283
284 resid[nMonomer*nBasis + counter] = -1 * scaleStress_ * str;
285 counter++;
286 }
287 }
288 UTIL_CHECK(counter == nFlexibleParams());
289 }
290
291 }
292
293 /*
294 * Update the current system field coordinates.
295 */
296 template <int D>
297 void AmIteratorBasis<D>::update(DArray<double>& newGuess)
298 {
299 UTIL_CHECK(system().domain().basis().isInitialized());
300 const int nMonomer = system().mixture().nMonomer();
301 const int nBasis = system().domain().basis().nBasis();
302
303 DArray< DArray<double> > wField;
304 wField.allocate(nMonomer);
305
306 // Restructure in format of monomers, basis functions
307 for (int i = 0; i < nMonomer; i++) {
308 wField[i].allocate(nBasis);
309 for (int k = 0; k < nBasis; k++) {
310 wField[i][k] = newGuess[i*nBasis + k];
311 }
312 }
313
314 // If canonical, explicitly set homogeneous field components
315 if (system().mixture().isCanonical()) {
316 double chi;
317 for (int i = 0; i < nMonomer; ++i) {
318 wField[i][0] = 0.0; // initialize to 0
319 for (int j = 0; j < nMonomer; ++j) {
320 chi = interaction_.chi(i,j);
321 wField[i][0] += chi * system().c().basis(j)[0];
322 }
323 }
324 // If system has external fields, include them in homogeneous part
325 if (system().h().hasData()) {
326 for (int i = 0; i < nMonomer; ++i) {
327 wField[i][0] += system().h().basis(i)[0];
328 }
329 }
330 }
331 system().w().setBasis(wField);
332
333 if (isFlexible_) {
334 const int nParam = system().domain().unitCell().nParameter();
335 const int begin = nMonomer*nBasis;
336
337 FSArray<double,6> parameters;
338 parameters = system().domain().unitCell().parameters();
339
340 const double coeff = 1.0 / scaleStress_;
341 int counter = 0;
342 for (int i = 0; i < nParam; i++) {
343 if (flexibleParams_[i]) {
344 parameters[i] = coeff * newGuess[begin + counter];
345 counter++;
346 }
347 }
348 UTIL_CHECK(counter == nFlexibleParams());
349
350 system().setUnitCell(parameters);
351 }
352 }
353
354 /*
355 * Output relevant system details to the iteration log file.
356 */
357 template<int D>
358 void AmIteratorBasis<D>::outputToLog()
359 {
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();
364 int counter = 0;
365 for (int i = 0; i < nParam; i++) {
366 if (flexibleParams_[i]) {
367 double str = residual()[nMonomer*nBasis + counter] /
368 (-1.0 * scaleStress_);
369 Log::file()
370 << " Cell Param " << i << " = "
371 << Dbl(system().domain().unitCell().parameters()[i], 15)
372 << " , stress = "
373 << Dbl(str, 15)
374 << "\n";
375 counter++;
376 }
377 }
378 }
379 }
380
381 #if 0
382 // Private virtual functions for vector math
383
384 /*
385 * Set a vector equal to another (assign a = b).
386 */
387 template <int D>
388 void AmIteratorBasis<D>::setEqual(DArray<double>& a,
389 DArray<double> const & b)
390 { a = b; }
391
392 /*
393 * Compute the inner product of two real vectors.
394 */
395 template <int D>
396 double AmIteratorBasis<D>::dotProduct(DArray<double> const & a,
397 DArray<double> const & b)
398 {
399 const int n = a.capacity();
400 UTIL_CHECK(b.capacity() == n);
401 double product = 0.0;
402 for (int i=0; i < n; ++i) {
403 // if either value is NaN, throw NanException
404 if (std::isnan(a[i]) || std::isnan(b[i])) {
405 throw NanException("AmIteratorBasis::dotProduct", __FILE__,
406 __LINE__, 0);
407 }
408 product += a[i] * b[i];
409 }
410 return product;
411 }
412
413 /*
414 * Compute and return maximum element of residual vector.
415 */
416 template <int D>
417 double AmIteratorBasis<D>::maxAbs(DArray<double> const & a)
418 {
419 const int n = a.capacity();
420 double max = 0.0;
421 double value;
422 for (int i = 0; i < n; i++) {
423 value = a[i];
424 if (std::isnan(value)) { // if value is NaN, throw NanException
425 throw NanException("AmIteratorBasis::dotProduct", __FILE__,
426 __LINE__, 0);
427 }
428 if (fabs(value) > max)
429 max = fabs(value);
430 }
431 return max;
432 }
433
434 /*
435 * Compute the vector difference a = b - c
436 */
437 template <int D>
438 void AmIteratorBasis<D>::subVV(DArray<double>& a,
439 DArray<double> const & b,
440 DArray<double> const & c)
441 {
442 const int n = a.capacity();
443 UTIL_CHECK(n == b.capacity());
444 UTIL_CHECK(n == c.capacity());
445 for (int i = 0; i < n; i++) {
446 a[i] = b[i] - c[i];
447 }
448 }
449
450 /*
451 * Composite a += b*c for vectors a and b, scalar c
452 */
453 template <int D>
454 void AmIteratorBasis<D>::addEqVc(DArray<double>& a,
455 DArray<double> const & b,
456 double c)
457 {
458 const int n = a.capacity();
459 UTIL_CHECK(n == b.capacity());
460 for (int i = 0; i < n; i++) {
461 a[i] += c*b[i];
462 }
463 }
464 #endif
465
466}
467}
468#endif
void readMixingParameters(std::istream &in, bool useLambdaRamp=true)
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.
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.
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
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.
Definition global.h:68
double max(Array< double > const &in)
Get maximum of array elements .
Definition Reduce.cpp:20
Periodic fields and crystallography.
Definition CField.cpp:11
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.
float product(float a, float b)
Product for float Data.
Definition product.h:22