PSCF v1.3
rpc/scft/sweep/Sweep.tpp
1#ifndef RPC_SWEEP_TPP
2/*
3* PSCF - Polymer Self-Consistent Field
4*
5* Copyright 2015 - 2025, The Regents of the University of Minnesota
6* Distributed under the terms of the GNU General Public License.
7*/
8
9#include "Sweep.h"
10#include <rpc/system/System.h>
11#include <rpc/scft/iterator/Iterator.h>
12#include <rpc/scft/ScftThermo.h>
13#include <prdc/environment/Environment.h>
14#include <pscf/inter/Interaction.h>
15#include <pscf/sweep/SweepTmpl.tpp>
16#include <util/misc/FileMaster.h>
17#include <util/misc/ioUtil.h>
18
19namespace Pscf {
20namespace Rpc {
21
22 using namespace Util;
23
24 // Maximum number of previous states = order of continuation + 1
25 # define RPC_HISTORY_CAPACITY 3
26
27 /*
28 * Default constructor (for unit testing).
29 */
30 template <int D>
32 : SweepTmpl< BasisFieldState<D> >(RPC_HISTORY_CAPACITY),
33 writeCRGrid_(false),
34 writeCBasis_(false),
35 writeWRGrid_(false),
36 systemPtr_(nullptr)
37 {}
38
39 /*
40 * Constructor, creates association with parent system.
41 */
42 template <int D>
44 : SweepTmpl< BasisFieldState<D> >(RPC_HISTORY_CAPACITY),
45 writeCRGrid_(false),
46 writeCBasis_(false),
47 writeWRGrid_(false),
48 systemPtr_(&sys)
49 {
50 // Get specialized sweep parameters from Environment
51 if (system().hasEnvironment()) {
52 addParameterTypes(system().environment().getParameterTypes());
53 }
54 }
55
56 /*
57 * Destructor.
58 */
59 template <int D>
62
63 /*
64 * Set association with a parent system (for unit testing).
65 */
66 template <int D>
68 { systemPtr_ = &system; }
69
70 /*
71 * Read parameters
72 */
73 template <int D>
74 void Sweep<D>::readParameters(std::istream& in)
75 {
76 // Call the base class's readParameters function.
78
79 // Read optional flags indicating which field types to output
80 readOptional(in, "writeCRGrid", writeCRGrid_);
81 readOptional(in, "writeCBasis", writeCBasis_);
82 readOptional(in, "writeWRGrid", writeWRGrid_);
83 }
84
85 /*
86 * Check allocation of one state object, allocate if necessary.
87 */
88 template <int D>
90 {
92 state.setSystem(system());
93 state.allocate();
94 state.unitCell() = system().domain().unitCell();
95 }
96
97 /*
98 * Setup operations at the beginning of a sweep.
99 */
100 template <int D>
102 {
103 initialize();
104 checkAllocation(trial_);
105
106 // Open log summary file
107 std::string fileName = baseFileName_;
108 fileName += "sweep.log";
109 system().fileMaster().openOutputFile(fileName, logFile_);
110 logFile_ << " step ds free_energy pressure"
111 << std::endl;
112 };
113
114 /*
115 * Set non-adjustable system parameters to new values.
116 *
117 * \param s path length coordinate, in range [0,1]
118 */
119 template <int D>
121 {
122 // Empty default implementation allows Sweep<D> to be compiled.
123 UTIL_THROW("Calling unimplemented function Sweep::setParameters");
124 };
125
126 /*
127 * Create guess for adjustable variables by polynomial extrapolation.
128 */
129 template <int D>
130 void Sweep<D>::extrapolate(double sNew)
131 {
132 UTIL_CHECK(historySize() > 0);
133
134 // If historySize() == 1, do nothing: Use previous system state
135 // as trial for the new state.
136
137 if (historySize() > 1) {
139
140 // Does the iterator allow a flexible unit cell ?
141 bool isFlexible = system().iterator().isFlexible();
142
143 // Compute coefficients of polynomial extrapolation to sNew
144 setCoefficients(sNew);
145
146 // Set extrapolated trial w fields
147 double coeff;
148 int nMonomer = system().mixture().nMonomer();
149 int nBasis = system().domain().basis().nBasis();
150 DArray<double>* newFieldPtr;
151 DArray<double>* oldFieldPtr;
152 int i, j, k;
153 for (i=0; i < nMonomer; ++i) {
154 newFieldPtr = &(trial_.field(i));
155
156 // Previous state k = 0 (most recent)
157 oldFieldPtr = &state(0).field(i);
158 coeff = c(0);
159 for (j=0; j < nBasis; ++j) {
160 (*newFieldPtr)[j] = coeff*(*oldFieldPtr)[j];
161 }
162
163 // Previous states k >= 1 (older)
164 for (k = 1; k < historySize(); ++k) {
165 oldFieldPtr = &state(k).field(i);
166 coeff = c(k);
167 for (j=0; j < nBasis; ++j) {
168 (*newFieldPtr)[j] += coeff*(*oldFieldPtr)[j];
169 }
170 }
171 }
172
173 // Make sure unitCellParameters_ is up to date with system
174 // (if we are sweeping in a lattice parameter, then the system
175 // parameters will be up-to-date but unitCellParameters_ wont be)
176 FSArray<double, 6> oldParameters = unitCellParameters_;
177 unitCellParameters_ = system().domain().unitCell().parameters();
178
179 // If isFlexible, extrapolate the flexible unit cell parameters
180 if (isFlexible) {
181
182 double coeff;
183 double parameter;
184
185 const FSArray<bool,6> flexParams =
186 system().iterator().flexibleParams();
187 const int nParameter
188 = system().domain().unitCell().nParameter();
189 UTIL_CHECK(flexParams.size() == nParameter);
190
191 // Add contributions from all previous states
192 for (k = 0; k < historySize(); ++k) {
193 coeff = c(k);
194 for (int i = 0; i < nParameter; ++i) {
195 if (flexParams[i]) {
196 if (k == 0) {
197 unitCellParameters_[i] = 0;
198 }
199 parameter = state(k).unitCell().parameter(i);
200 unitCellParameters_[i] += coeff*parameter;
201 }
202 }
203 }
204 }
205
206 // Reset trial_.unitCell() object
207 trial_.unitCell().setParameters(unitCellParameters_);
208
209 // Check if any unit cell parameters have changed
210 bool newCellParams(true);
211 for (int i = 0; i < oldParameters.size(); i++) {
212 if (fabs(oldParameters[i] - unitCellParameters_[i]) < 1e-10) {
213 newCellParams = false;
214 break;
215 }
216 }
217
218 // Transfer data from trial_ state to parent system
219 trial_.setSystemState(newCellParams);
220 }
221
222
223 };
224
225 /*
226 * Call current iterator to solve SCFT problem.
227 *
228 * Return 0 for sucessful solution, 1 on failure to converge.
229 */
230 template <int D>
231 int Sweep<D>::solve(bool isContinuation)
232 { return system().iterate(isContinuation); };
233
234 /*
235 * Reset system to previous solution after iterature failure.
236 *
237 * The implementation of this function should reset the system state
238 * to correspond to that stored in state(0).
239 */
240 template <int D>
242 {
243 bool isFlexible = system().iterator().isFlexible();
244 state(0).setSystemState(isFlexible);
245 }
246
247 /*
248 * Update state(0) and output data after successful convergence
249 *
250 * The implementation of this function should copy the current
251 * system state into state(0) and output any desired information
252 * about the current converged solution.
253 */
254 template <int D>
256 {
257 state(0).setSystem(system());
258 state(0).getSystemState();
259
260 // Output converged solution to several files
261 outputSolution();
262
263 // Output summary to log file
264 outputSummary(logFile_);
265
266 };
267
268 template <int D>
269 void Sweep<D>::outputSolution()
270 {
271 std::ofstream out;
272 std::string outFileName;
273 std::string indexString = toString(nAccept() - 1);
274
275 // Open parameter file, with thermodynamic properties at end
276 outFileName = baseFileName_;
277 outFileName += indexString;
278 outFileName += ".stt";
279 system().fileMaster().openOutputFile(outFileName, out);
280
281 // Write data file, with thermodynamic properties at end
282 system().writeParamNoSweep(out);
283 out << std::endl;
284 system().scft().write(out);
285 out.close();
286
287 FieldIo<D> const & fieldIo = system().domain().fieldIo();
288 UnitCell<D> const & unitCell = system().domain().unitCell();
289
290 // Write w fields
291 outFileName = baseFileName_;
292 outFileName += indexString;
293 outFileName += "_w";
294 outFileName += ".bf";
295 UTIL_CHECK(system().w().hasData());
296 UTIL_CHECK(system().w().isSymmetric());
297 fieldIo.writeFieldsBasis(outFileName,
298 system().w().basis(), unitCell);
299
300 // Optionally write c rgrid files
301 if (writeCRGrid_) {
302 outFileName = baseFileName_;
303 outFileName += indexString;
304 outFileName += "_c";
305 outFileName += ".rf";
306 fieldIo.writeFieldsRGrid(outFileName,
307 system().c().rgrid(), unitCell);
308 }
309
310 // Optionally write c basis files
311 if (writeCBasis_) {
312 outFileName = baseFileName_;
313 outFileName += indexString;
314 outFileName += "_c";
315 outFileName += ".bf";
316 UTIL_CHECK(system().c().hasData());
317 fieldIo.writeFieldsBasis(outFileName,
318 system().c().basis(), unitCell);
319
320 }
321
322 // Optionally write w rgrid files
323 if (writeWRGrid_) {
324 outFileName = baseFileName_;
325 outFileName += indexString;
326 outFileName += "_w";
327 outFileName += ".rf";
328 fieldIo.writeFieldsRGrid(outFileName,
329 system().w().rgrid(), unitCell);
330 }
331
332 }
333
334 template <int D>
335 void Sweep<D>::outputSummary(std::ostream& out)
336 {
337 int i = nAccept() - 1;
338 double sNew = s(0);
339 if (!system().scft().hasData()) system().scft().compute();
340 out << Int(i,5) << Dbl(sNew)
341 << Dbl(system().scft().fHelmholtz(),16)
342 << Dbl(system().scft().pressure(),16);
343 out << std::endl;
344 }
345
346 template <int D>
348 { logFile_.close(); }
349
350} // namespace Rpc
351} // namespace Pscf
352#endif
Base template for UnitCell<D> classes, D=1, 2 or 3.
Definition UnitCell.h:56
FieldState for fields in symmetry-adapted basis format.
File input/output operations and format conversions for fields.
void writeFieldsRGrid(std::ostream &out, DArray< RField< D > > const &fields, UnitCell< D > const &unitCell, bool writeHeader=true, bool isSymmetric=true, bool writeMeshSize=true) const override
Write array of RField objects (fields on r-space grid) to a stream.
void writeFieldsBasis(std::ostream &out, DArray< DArray< double > > const &fields, UnitCell< D > const &unitCell) const
Write an array of fields in basis format to an output stream.
virtual int solve(bool isContinuation)
Call current iterator to solve SCFT problem.
virtual void extrapolate(double sNew)
Create a guess for adjustable variables by continuation.
void setSystem(System< D > &system)
Set association with parent System.
Sweep()
Default Constructor.
virtual void getSolution()
Update state(0) and output data after successful convergence.
virtual void checkAllocation(BasisFieldState< D > &state)
Check allocation of fields in one state, allocate if necessary.
System< D > & system()
Return the parent system by reference.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
virtual void reset()
Reset system to previous solution after iterature failure.
bool writeCRGrid_
Should concentration fields be written to file in r-grid format?
virtual void cleanup()
Cleanup operation at the beginning of a sweep.
bool writeCBasis_
Should concentration fields be written to file in basis format?
virtual void readParameters(std::istream &in)
Read parameters from param file.
virtual void setParameters(double sNew)=0
Set system parameters to new values.
bool writeWRGrid_
Should converged w fields be written to file in r-grid format?
bool hasSystem()
Does an association with the parent System exist?
virtual void setup()
Setup operation at the beginning of a sweep.
Main class, representing one complete system.
void addParameterTypes(GArray< ParameterType > paramTypes)
BasisFieldState< D > & state(int i)
Definition SweepTmpl.h:127
Dynamically allocatable contiguous array template.
Definition DArray.h:32
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
int size() const
Return logical size of this array (i.e., number of elements).
Definition FSArray.h:207
Wrapper for an int, for formatted ostream output.
Definition Int.h:37
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Definition global.h:49
std::string toString(int n)
Return string representation of an integer.
Definition ioUtil.cpp:52
Real periodic fields, SCFT and PS-FTS (CPU).
Definition param_pc.dox:2
PSCF package top-level namespace.
Definition param_pc.dox:1