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