PSCF v1.3.2
rpc/scft/iterator/AmIteratorBasis.tpp
1#ifndef RPC_AM_ITERATOR_BASIS_TPP
2#define RPC_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 <rpc/system/System.h>
13#include <rpc/solvers/Mixture.h>
14#include <rpc/field/Domain.h>
15#include <prdc/crystal/UnitCell.h>
16#include <pscf/inter/Interaction.h>
17#include <pscf/iterator/NanException.h>
18#include <util/containers/DArray.h>
19#include <util/global.h>
20#include <cmath>
21
22namespace Pscf {
23namespace Rpc {
24
25 using namespace Util;
26 using namespace Prdc;
27
28 // Public member functions
29
30 /*
31 * Constructor.
32 */
33 template <int D>
35 : Iterator<D>(system)
36 {
37 isSymmetric_ = true;
38 ParamComposite::setClassName("AmIteratorBasis");
39 }
40
41 /*
42 * Destructor.
43 */
44 template <int D>
47
48 /*
49 * Read parameters from file.
50 */
51 template <int D>
53 {
54 // Use base class methods to read parameters
57
58 UnitCell<D> const & unitCell = system().domain().unitCell();
59 UTIL_CHECK(unitCell.lattice() != UnitCell<D>::Null);
60 int np = unitCell.nParameter();
61 UTIL_CHECK(np > 0);
62 UTIL_CHECK(np <= 6);
63
64 // Read optional isFlexible boolean (true by default)
65 isFlexible_ = 1; // Default
66 readOptional(in, "isFlexible", isFlexible_);
67
68 // Populate flexibleParams_ based on isFlexible_ (all 0s or all 1s),
69 // then optionally overwrite with user input from param file
70 if (isFlexible_) {
71 flexibleParams_.clear();
72 for (int i = 0; i < np; i++) {
73 flexibleParams_.append(true); // Set all values to true
74 }
75 // Read optional flexibleParams_ array to overwrite current array
76 readOptionalFSArray(in, "flexibleParams", flexibleParams_, np);
77 if (nFlexibleParams() == 0) isFlexible_ = false;
78 } else { // isFlexible_ = false
79 flexibleParams_.clear();
80 for (int i = 0; i < np; i++) {
81 flexibleParams_.append(false); // Set all values to false
82 }
83 }
84
85 // Read optional scaleStress value
86 scaleStress_ = 10.0; // Default value
87 readOptional(in, "scaleStress", scaleStress_);
88
89 // Read option mixing parameters (lambda, useLambdaRamp, and r)
91
92 // Allocate local modified copy of Interaction class
93 const int nMonomer = system().mixture().nMonomer();
94 interaction_.setNMonomer(nMonomer);
95
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 out << "\n";
105 out << "Iterator times contributions:\n";
107 }
108
109 // Protected virtual function
110
111 // Setup before entering iteration loop
112 template <int D>
113 void AmIteratorBasis<D>::setup(bool isContinuation)
114 {
115 AmTmpl::setup(isContinuation);
116 interaction_.update(system().interaction());
117 }
118
119 // Private virtual functions to exchange data with parent system
120
121 /*
122 * Compute and return number of elements in a residual vector.
123 */
124 template <int D>
125 int AmIteratorBasis<D>::nElements()
126 {
127 const int nMonomer = system().mixture().nMonomer();
128 const int nBasis = system().domain().basis().nBasis();
129
130 int nEle = nMonomer*nBasis;
131 if (isFlexible()) {
132 nEle += nFlexibleParams();
133 }
134 return nEle;
135 }
136
137 /*
138 * Does the system have an initial field guess?
139 */
140 template <int D>
141 bool AmIteratorBasis<D>::hasInitialGuess()
142 { return system().w().hasData(); }
143
144 /*
145 * Get the current field vector (w fields and lattice parameters).
146 */
147 template <int D>
148 void AmIteratorBasis<D>::getCurrent(DArray<double>& state)
149 {
150 const int nMonomer = system().mixture().nMonomer();
151 const int nBasis = system().domain().basis().nBasis();
152 const int nEle = nElements();
153 UTIL_CHECK(state.capacity() == nEle);
154
155 // Copy w-fieldd into linear array
156 int begin;
157 for (int i = 0; i < nMonomer; i++) {
158 DArray<double> const & field = system().w().basis(i);
159 begin = i*nBasis;
160 for (int k = 0; k < nBasis; k++) {
161 state[begin + k] = field[k];
162 }
163 }
164
165 // Add elements associated with unit cell parameters (if any)
166 if (isFlexible()) {
167 UTIL_CHECK(nFlexibleParams() > 0);
168 UnitCell<D> const & unitCell = system().domain().unitCell();
169 FSArray<double,6> const & parameters = unitCell.parameters();
170 const int nParam = unitCell.nParameter();
171 begin = nMonomer*nBasis;
172 int counter = 0;
173 for (int i = 0; i < nParam; i++) {
174 if (flexibleParams_[i]) {
175 state[begin + counter] = scaleStress_*parameters[i];
176 counter++;
177 }
178 }
179 UTIL_CHECK(counter == nFlexibleParams());
180 }
181
182 }
183
184 /*
185 * Perform the main system computation (solve the MDE).
186 */
187 template <int D>
188 void AmIteratorBasis<D>::evaluate()
189 { system().compute(isFlexible_); }
190
191 /*
192 * Compute the residual for the current system state.
193 */
194 template <int D>
195 void AmIteratorBasis<D>::getResidual(DArray<double>& resid)
196 {
197 const int nMonomer = system().mixture().nMonomer();
198 const int nBasis = system().domain().basis().nBasis();
199 const int n = nElements();
200 UTIL_CHECK(resid.capacity() == n);
201
202 // Local variables
203 double chi, p;
204 int i, j, k, begin;
205
206 // Initialize residual vector to zero
207 for (i = 0 ; i < n; ++i) {
208 resid[i] = 0.0;
209 }
210
211 // Compute SCF residual vector elements
212 for (i = 0; i < nMonomer; ++i) {
213 begin = i*nBasis;
214 for (j = 0; j < nMonomer; ++j) {
215 chi = interaction_.chi(i,j);
216 p = interaction_.p(i,j);
217 DArray<double> const & c = system().c().basis(j);
218 DArray<double> const & w = system().w().basis(j);
219 if (system().h().hasData()) {
220 DArray<double> const & h = system().h().basis(j);
221 for (k = 0; k < nBasis; ++k) {
222 resid[begin + k] += chi*c[k] + p*(h[k] - w[k]);
223 }
224 } else {
225 for (k = 0; k < nBasis; ++k) {
226 resid[begin + k] += chi*c[k] - p*w[k];
227 }
228 }
229 }
230 }
231
232 // Add term proportional to sum of all concentrations
233 double shift = -1.0 / interaction_.sumChiInverse();
234 if (system().mask().hasData()) {
235 DArray<double> const & mask = system().mask().basis();
236 for (i = 0; i < nMonomer; ++i) {
237 begin = i*nBasis;
238 for (k = 0; k < nBasis; ++k) {
239 resid[begin + k] += shift*mask[k];
240 }
241 }
242 } else {
243 for (i = 0; i < nMonomer; ++i) {
244 resid[i*nBasis] += shift;
245 }
246 }
247
248 // If canonical ensemble, zero homogeneous residual components
249 if (system().mixture().isCanonical()) {
250 for (i = 0; i < nMonomer; ++i) {
251 resid[i*nBasis] = 0.0;
252 }
253 }
254
255 // If flexible unit cell, then compute stress residuals
256 if (isFlexible()) {
257
258 // Combined -1 factor and stress scaling here. This is okay:
259 // - residuals only show up as dot products (U, v, norm)
260 // or with their absolute value taken (max), so the
261 // sign on a given residual vector element is not relevant
262 // as long as it is consistent across all vectors
263 // - The scaling is applied here and to the unit cell param
264 // storage, so that updating is done on the same scale,
265 // and then undone right before passing to the unit cell.
266
267 double coeff = -1.0 * scaleStress_;
268 const int nParam = system().domain().unitCell().nParameter();
269 begin = nMonomer*nBasis;
270 int counter = 0;
271 for (i = 0; i < nParam ; i++) {
272 if (flexibleParams_[i]) {
273 resid[begin + counter] = coeff * stress(i);
274 counter++;
275 }
276 }
277 UTIL_CHECK(counter == nFlexibleParams());
278 }
279
280 }
281
282 /*
283 * Update the current system field vector.
284 */
285 template <int D>
286 void AmIteratorBasis<D>::update(DArray<double>& newState)
287 {
288 // Constants
289 const int nMonomer = system().mixture().nMonomer();
290 const int nBasis = system().domain().basis().nBasis();
291
292 // Allocate wFields container
293 DArray< DArray<double> > wFields;
294 wFields.allocate(nMonomer);
295 for (int i = 0; i < nMonomer; i++) {
296 wFields[i].allocate(nBasis);
297 }
298
299 // Copy w fields from newState to wFields container
300 int begin;
301 for (int i = 0; i < nMonomer; i++) {
302 begin = i*nBasis;
303 for (int k = 0; k < nBasis; k++) {
304 wFields[i][k] = newState[begin + k];
305 }
306 }
307
308 // If canonical, explicitly set homogeneous field components
309 if (system().mixture().isCanonical()) {
310
311 // Set homogeneous components of all w fields to zero
312 for (int i = 0; i < nMonomer; ++i) {
313 wFields[i][0] = 0.0;
314 }
315
316 // Add average values arising from interactions
317 double chi, wAve, cAve;
318 for (int i = 0; i < nMonomer; ++i) {
319 wAve = 0.0;
320 for (int j = 0; j < nMonomer; ++j) {
321 chi = interaction_.chi(i,j);
322 cAve = system().c().basis(j)[0];
323 wAve += chi * cAve;
324 }
325 wFields[i][0] = wAve;
326 }
327
328 // If external fields exist, add their spatial averages
329 if (system().h().hasData()) {
330 for (int i = 0; i < nMonomer; ++i) {
331 wFields[i][0] += system().h().basis(i)[0];
332 }
333 }
334 }
335
336 // Set fields in system w container
337 system().w().setBasis(wFields);
338
339 // If flexible, update unit cell parameters
340 if (isFlexible()) {
341
342 // Initialize parameters array with current values
343 FSArray<double, 6> parameters;
344 parameters = system().domain().unitCell().parameters();
345
346 // Reset any parameters that are flexible
347 const double coeff = 1.0 / scaleStress_;
348 const int nParam = system().domain().unitCell().nParameter();
349 const int begin = nMonomer*nBasis;
350 int counter = 0;
351 for (int i = 0; i < nParam; i++) {
352 if (flexibleParams_[i]) {
353 parameters[i] = coeff * newState[begin + counter];
354 counter++;
355 }
356 }
357 UTIL_CHECK(counter == nFlexibleParams());
358
359 // Set system unit cell parameters
360 system().setUnitCell(parameters);
361 }
362 }
363
364 /*
365 * Output relevant system details to the iteration log.
366 */
367 template<int D>
368 void AmIteratorBasis<D>::outputToLog()
369 {
370 if (isFlexible() && verbose() > 1) {
371 double str;
372 UnitCell<D> const & unitCell = system().domain().unitCell();
373 const int nParam = unitCell.nParameter();
374 const int nMonomer = system().mixture().nMonomer();
375 const int nBasis = system().domain().basis().nBasis();
376 const int begin = nMonomer*nBasis;
377 int counter = 0;
378 for (int i = 0; i < nParam; i++) {
379 if (flexibleParams_[i]) {
380 str = -1.0 * residual()[begin + counter] / scaleStress_;
381 Log::file()
382 << " Cell Param " << i << " = "
383 << Dbl(unitCell.parameters()[i], 15)
384 << " , stress = "
385 << Dbl(str, 15)
386 << "\n";
387 counter++;
388 }
389 }
390 }
391 }
392
393}
394}
395#endif
void readMixingParameters(std::istream &in, bool useLambdaRamp=true)
int nParameter() const
Get the number of parameters in the unit cell.
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 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.
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.
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
Periodic fields and crystallography.
Definition CField.cpp:11
Real periodic fields, SCFT and PS-FTS (CPU).
Definition param_pc.dox:2
PSCF package top-level namespace.