PSCF v1.3
rpc/fts/simulator/Simulator.tpp
1#ifndef RPC_SIMULATOR_TPP
2#define RPC_SIMULATOR_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 "Simulator.h"
12
13#include <rpc/system/System.h>
14#include <rpc/fts/compressor/Compressor.h>
15#include <rpc/fts/compressor/CompressorFactory.h>
16#include <rpc/fts/perturbation/Perturbation.h>
17#include <rpc/fts/perturbation/PerturbationFactory.h>
18#include <rpc/fts/ramp/Ramp.h>
19#include <rpc/fts/ramp/RampFactory.h>
20
21#include <pscf/inter/Interaction.h>
22
23#include <util/misc/Timer.h>
24#include <util/random/Random.h>
25#include <util/global.h>
26
27// Gnu scientifie library
28#include <gsl/gsl_eigen.h>
29
30namespace Pscf {
31namespace Rpc {
32
33 using namespace Util;
34
35 /*
36 * Constructor.
37 */
38 template <int D>
39 Simulator<D>::Simulator(System<D>& system)
40 : random_(),
41 hamiltonian_(0.0),
45 iStep_(0),
46 iTotalStep_(0),
47 seed_(0),
48 hasHamiltonian_(false),
49 hasWc_(false),
50 hasCc_(false),
51 hasDc_(false),
52 systemPtr_(&system),
53 compressorFactoryPtr_(nullptr),
54 compressorPtr_(nullptr),
55 perturbationFactoryPtr_(nullptr),
56 perturbationPtr_(nullptr),
57 rampFactoryPtr_(nullptr),
58 rampPtr_(nullptr),
59 isAllocated_(false)
60 {
61 setClassName("Simulator");
62 compressorFactoryPtr_ = new CompressorFactory<D>(system);
63 perturbationFactoryPtr_ = new PerturbationFactory<D>(*this);
64 rampFactoryPtr_ = new RampFactory<D>(*this);
65 }
66
67 /*
68 * Destructor.
69 */
70 template <int D>
72 {
73 if (compressorFactoryPtr_) {
74 delete compressorFactoryPtr_;
75 }
76 if (compressorPtr_ ) {
77 delete compressorPtr_;
78 }
79 if (perturbationFactoryPtr_) {
80 delete perturbationFactoryPtr_;
81 }
82 if (perturbationPtr_) {
83 delete perturbationPtr_;
84 }
85 if (rampFactoryPtr_) {
86 delete rampFactoryPtr_;
87 }
88 if (rampPtr_) {
89 delete rampPtr_;
90 }
91 }
92
93 /*
94 * Allocate required memory.
95 */
96 template <int D>
98 {
99 UTIL_CHECK(!isAllocated_);
100
101 const int nMonomer = system().mixture().nMonomer();
102
103 // Allocate projected chi matrix chiP_ and associated arrays
104 chiP_.allocate(nMonomer, nMonomer);
105 chiEvals_.allocate(nMonomer);
106 chiEvecs_.allocate(nMonomer, nMonomer);
107 sc_.allocate(nMonomer);
108
109 // Allocate memory for eignevector components of w and c fields
110 wc_.allocate(nMonomer);
111 cc_.allocate(nMonomer);
112 const IntVec<D> dimensions = system().domain().mesh().dimensions();
113 for (int i = 0; i < nMonomer; ++i) {
114 wc_[i].allocate(dimensions);
115 cc_[i].allocate(dimensions);
116 }
117
118 // Allocate memory for components of d (functional derivative)
119 dc_.allocate(nMonomer-1);
120 for (int i = 0; i < nMonomer - 1; ++i) {
121 dc_[i].allocate(dimensions);
122 }
123
124 // Allocate state_, if necessary.
125 if (!state_.isAllocated) {
126 const IntVec<D> dimensions = system().domain().mesh().dimensions();
127 state_.allocate(nMonomer, dimensions);
128 }
129
130 isAllocated_ = true;
131 }
132
133 /*
134 * Default implementation - designed to be used by subclasses to read the
135 * initial common part of the parameter file block.
136 */
137 template <int D>
138 void Simulator<D>::readParameters(std::istream &in)
139 {
140 // Optionally read a random number generator seed
141 readRandomSeed(in);
142
143 bool isEnd = false;
144
145 // Read required Compressor block
146 readCompressor(in, isEnd);
147
148 // Optionally read a Perturbation
149 readPerturbation(in, isEnd);
150
151 // Optionally read a Ramp
152 readRamp(in, isEnd);
153 }
154
155 /*
156 * Perform a field theoretic simulation (unimplemented by base class).
157 */
158 template <int D>
160 { UTIL_THROW("Error: Unimplemented function Simulator<D>::simulate"); }
161
162 /*
163 * Open, read and analyze a trajectory file (unimplemented by base class).
164 */
165 template <int D>
166 void Simulator<D>::analyze(int min, int max,
167 std::string classname,
168 std::string filename)
169 { UTIL_THROW("Error: Unimplemented function Simulator<D>::analyze"); }
170
171 /*
172 * Clear all local state data (eigen-components of w and Hamiltonian)
173 */
174 template <int D>
176 {
177 hasHamiltonian_ = false;
178 hasWc_ = false;
179 hasCc_ = false;
180 hasDc_ = false;
181 }
182
183 /*
184 * Compute field theoretic Hamiltonian H[W].
185 */
186 template <int D>
188 {
189 UTIL_CHECK(isAllocated_);
190 UTIL_CHECK(system().w().hasData());
191 UTIL_CHECK(system().c().hasData());
193 hasHamiltonian_ = false;
194
195 Mixture<D> const & mixture = system().mixture();
196 Domain<D> const & domain = system().domain();
197
198 const int nMonomer = mixture.nMonomer();
199 const int meshSize = domain.mesh().size();
200
201 const int np = mixture.nPolymer();
202 const int ns = mixture.nSolvent();
203 double phi, mu;
204
205 // Compute polymer ideal gas contributions to lnQ
206 double lnQ = 0.0;
207 if (np > 0) {
208 Polymer<D> const * polymerPtr;
209 double length;
210 for (int i = 0; i < np; ++i) {
211 polymerPtr = &mixture.polymer(i);
212 phi = polymerPtr->phi();
213 mu = polymerPtr->mu();
215 length = polymerPtr->length();
216 } else {
217 length = (double)polymerPtr->nBead();
218 }
219 // Recall: mu = ln(phi/q)
220 if (phi > 1.0E-08) {
221 lnQ += phi*( -mu + 1.0 )/length;
222 }
223 }
224 }
225
226 // Compute solvent ideal gas contributions to lnQ
227 if (ns > 0) {
228 Solvent<D> const * solventPtr;
229 double size;
230 for (int i = 0; i < ns; ++i) {
231 solventPtr = &mixture.solvent(i);
232 phi = solventPtr->phi();
233 mu = solventPtr->mu();
234 size = solventPtr->size();
235 // Recall: mu = ln(phi/q)
236 if (phi > 1.0E-8) {
237 lnQ += phi*( -mu + 1.0 )/size;
238 }
239 }
240 }
241
242 // Subtract average of pressure field wc_[nMonomer-1]
243 RField<D> const & Wc = wc_[nMonomer-1];
244 for (int i = 0; i < meshSize; ++i) {
245 lnQ += Wc[i]/double(meshSize);
246 }
247 // lnQ now contains a value per monomer
248
249 // Initialize field contribution HW
250
251 // Compute quadratic field contribution to HW
252 double HW = 0.0;
253 double prefactor, w, s;
254 int i, j;
255 for (j = 0; j < nMonomer - 1; ++j) {
256 RField<D> const & Wc = wc_[j];
257 prefactor = -0.5*double(nMonomer)/chiEvals_[j];
258 s = sc_[j];
259 for (i = 0; i < meshSize; ++i) {
260 w = Wc[i] - s;
261 HW += prefactor*w*w;
262 }
263 }
264
265 // Normalize HW to equal a value per monomer
266 HW /= double(meshSize);
267
268 // Add constant term K/2 per monomer (K=s=e^{T}chi e/M^2)
269 HW += 0.5*sc_[nMonomer - 1];
270
271 // Compute number of monomers in the system (nMonomerSystem)
272 const double vSystem = domain.unitCell().volume();
273 const double vMonomer = mixture.vMonomer();
274 const double nMonomerSystem = vSystem / vMonomer;
275
276 // Compute final Hamiltonian components
277 fieldHamiltonian_ = nMonomerSystem * HW;
278 idealHamiltonian_ = -1.0 * nMonomerSystem * lnQ;
280
281 if (hasPerturbation()) {
284 } else {
286 }
287
288 hasHamiltonian_ = true;
289 }
290
291 template <int D>
293 {
294 UTIL_CHECK(isAllocated_);
295
296 const int nMonomer = system().mixture().nMonomer();
297 DMatrix<double> const & chi = system().interaction().chi();
298 double d = 1.0/double(nMonomer);
299 int i, j, k;
300
301 // Compute orthogonal projection matrix P
303 P.allocate(nMonomer, nMonomer);
304 for (i = 0; i < nMonomer; ++i) {
305 for (j = 0; j < nMonomer; ++j) {
306 P(i,j) = -d;
307 }
308 P(i,i) += 1.0;
309 }
310
311 // Compute T = chi*P (temporary matrix)
313 T.allocate(nMonomer, nMonomer);
314 for (i = 0; i < nMonomer; ++i) {
315 for (j = 0; j < nMonomer; ++j) {
316 T(i, j) = 0.0;
317 for (k = 0; k < nMonomer; ++k) {
318 T(i,j) += chi(i,k)*P(k,j);
319 }
320 }
321 }
322
323 // Compute chiP = = P*chi*P = P*T
324 DMatrix<double> chiP;
325 chiP.allocate(nMonomer, nMonomer);
326 for (i = 0; i < nMonomer; ++i) {
327 for (j = 0; j < nMonomer; ++j) {
328 chiP(i, j) = 0.0;
329 for (k = 0; k < nMonomer; ++k) {
330 chiP(i,j) += P(i,k)*T(k,j);
331 }
332 }
333 }
334
335 // Eigenvalue calculations use data structures and
336 // functions from the Gnu Scientific Library (GSL)
337
338 // Allocate GSL matrix A that will hold a copy of chiP
339 gsl_matrix* A = gsl_matrix_alloc(nMonomer, nMonomer);
340
341 // Copy DMatrix<double> chiP to gsl_matrix A
342 for (i = 0; i < nMonomer; ++i) {
343 for (j = 0; j < nMonomer; ++j) {
344 gsl_matrix_set(A, i, j, chiP(i, j));
345 }
346 }
347
348 // Compute eigenvalues and eigenvectors of chiP (or A)
349 gsl_eigen_symmv_workspace* work = gsl_eigen_symmv_alloc(nMonomer);
350 gsl_vector* Avals = gsl_vector_alloc(nMonomer);
351 gsl_matrix* Avecs = gsl_matrix_alloc(nMonomer, nMonomer);
352 int error;
353 error = gsl_eigen_symmv(A, Avals, Avecs, work);
354 UTIL_CHECK(error == 0);
355
356 // Requirements:
357 // - A has exactly one zero eigenvalue, with eigenvector (1,...,1)
358 // - All other eigenvalues must be negative.
359
360 // Copy eigenpairs with non-null eigenvalues
361 int iNull = -1; // index for null eigenvalue
362 int nNull = 0; // number of null eigenvalue
363 k = 0; // re-ordered index for non-null eigenvalue
364 double val;
365 for (i = 0; i < nMonomer; ++i) {
366 val = gsl_vector_get(Avals, i);
367 if (std::abs(val) < 1.0E-8) {
368 ++nNull;
369 iNull = i;
370 UTIL_CHECK(nNull <= 1);
371 } else {
372 chiEvals_[k] = val;
373 UTIL_CHECK(val < 0.0);
374 for (j = 0; j < nMonomer; ++j) {
375 chiEvecs_(k, j) = gsl_matrix_get(Avecs, j, i);
376 }
377 if (chiEvecs_(k, 0) < 0.0) {
378 for (j = 0; j < nMonomer; ++j) {
379 chiEvecs_(k, j) = -chiEvecs_(k, j);
380 }
381 }
382 ++k;
383 }
384 }
385 UTIL_CHECK(nNull == 1);
386 UTIL_CHECK(iNull >= 0);
387
388 // Set eigenpair with zero eigenvalue
389 i = nMonomer - 1;
390 chiEvals_[i] = 0.0;
391 for (j = 0; j < nMonomer; ++j) {
392 chiEvecs_(i, j) = gsl_matrix_get(Avecs, j, iNull);
393 }
394 if (chiEvecs_(i, 0) < 0) {
395 for (j = 0; j < nMonomer; ++j) {
396 chiEvecs_(i, j) = -chiEvecs_(i, j);
397 }
398 }
399
400 // Normalize all eigenvectors so that the sum of squares = nMonomer
401 double vec, norm, prefactor;
402 for (i = 0; i < nMonomer; ++i) {
403 norm = 0.0;
404 for (j = 0; j < nMonomer; ++j) {
405 vec = chiEvecs_(i, j);
406 norm += vec*vec;
407 }
408 prefactor = sqrt( double(nMonomer)/norm );
409 for (j = 0; j < nMonomer; ++j) {
410 chiEvecs_(i, j) *= prefactor;
411 }
412 }
413
414 // Check final eigenvector is (1, ..., 1)
415 for (j = 0; j < nMonomer; ++j) {
416 UTIL_CHECK(std::abs(chiEvecs_(nMonomer-1, j) - 1.0) < 1.0E-8);
417 }
418
419 // Compute vector s in monomer basis
421 s.allocate(nMonomer);
422 for (i = 0; i < nMonomer; ++i) {
423 s[i] = 0.0;
424 for (j = 0; j < nMonomer; ++j) {
425 s[i] += chi(i,j);
426 }
427 s[i] = s[i]/double(nMonomer);
428 }
429
430 // Compute components of s in eigenvector basis -> sc_
431 for (i = 0; i < nMonomer; ++i) {
432 sc_[i] = 0.0;
433 for (j = 0; j < nMonomer; ++j) {
434 sc_[i] += chiEvecs_(i,j)*s[j];
435 }
436 sc_[i] = sc_[i]/double(nMonomer);
437 }
438
439 #if 0
440 // Debugging output
441 for (i = 0; i < nMonomer; ++i) {
442 Log::file() << "Eigenpair " << i << "\n";
443 Log::file() << "value = " << chiEvals_[i] << "\n";
444 Log::file() << "vector = [ ";
445 for (j = 0; j < nMonomer; ++j) {
446 Log::file() << chiEvecs_(i, j) << " ";
447 }
448 Log::file() << "]\n";
449 Log::file() << " sc[i] = " << sc_[i] << std::endl;
450 }
451 #endif
452
453 }
454
455 /*
456 * Compute the eigenvector components of the w fields, using the
457 * eigenvectors chiEvecs_ of the projected chi matrix as a basis.
458 */
459 template <int D>
461 {
462 UTIL_CHECK(isAllocated_);
463
464 const int nMonomer = system().mixture().nMonomer();
465 const int meshSize = system().domain().mesh().size();
466 int i, j, k;
467
468 // Loop over eigenvectors (j is an eigenvector index)
469 for (j = 0; j < nMonomer; ++j) {
470
471 // Loop over grid points to zero out field wc_[j]
472 RField<D>& Wc = wc_[j];
473 for (i = 0; i < meshSize; ++i) {
474 Wc[i] = 0.0;
475 }
476
477 // Loop over monomer types (k is a monomer index)
478 for (k = 0; k < nMonomer; ++k) {
479 double vec = chiEvecs_(j, k)/double(nMonomer);
480
481 // Loop over grid points
482 RField<D> const & Wr = system().w().rgrid(k);
483 for (i = 0; i < meshSize; ++i) {
484 Wc[i] += vec*Wr[i];
485 }
486
487 }
488 }
489
490 hasWc_ = true;
491 }
492
493 /*
494 * Compute the eigenvector components of the c-fields, using the
495 * eigenvectors chiEvecs_ of the projected chi matrix as a basis.
496 */
497 template <int D>
499 {
500 // Preconditions
501 UTIL_CHECK(isAllocated_);
502 UTIL_CHECK(system().w().hasData());
503 UTIL_CHECK(system().c().hasData());
504
505 const int nMonomer = system().mixture().nMonomer();
506 const int meshSize = system().domain().mesh().size();
507 int i, j, k;
508
509 // Loop over eigenvectors (i is an eigenvector index)
510 for (i = 0; i < nMonomer; ++i) {
511
512 // Set cc_[i] to zero
513 RField<D>& Cc = cc_[i];
514 for (k = 0; k < meshSize; ++k) {
515 Cc[k] = 0.0;
516 }
517
518 // Loop over monomer types
519 for (j = 0; j < nMonomer; ++j) {
520 RField<D> const & Cr = system().c().rgrid(j);
521 double vec = chiEvecs_(i, j);
522
523 // Loop over grid points
524 for (k = 0; k < meshSize; ++k) {
525 Cc[k] += vec*Cr[k];
526 }
527
528 }
529 }
530
531 hasCc_ = true;
532 }
533
534 /*
535 * Compute d fields, i.e., functional derivatives of H[W].
536 */
537 template <int D>
539 {
540 // Preconditions
541 UTIL_CHECK(isAllocated_);
542 if (!hasWc_) computeWc();
543 if (!hasCc_) computeCc();
544
545 // Local constants and variables
546 const int meshSize = system().domain().mesh().size();
547 const int nMonomer = system().mixture().nMonomer();
548 const double vMonomer = system().mixture().vMonomer();
549 const double a = 1.0/vMonomer;
550 double b, s;
551 int i, k;
552
553 // Compute derivatives for standard Hamiltonian
554 // Loop over composition eigenvectors (exclude the last)
555 for (i = 0; i < nMonomer - 1; ++i) {
556 RField<D>& Dc = dc_[i];
557 RField<D> const & Wc = wc_[i];
558 RField<D> const & Cc = cc_[i];
559 b = -1.0*double(nMonomer)/chiEvals_[i];
560 s = sc_[i];
561 // Loop over grid points
562 for (k = 0; k < meshSize; ++k) {
563 Dc[k] = a*( b*(Wc[k] - s) + Cc[k] );
564 }
565 }
566
567 // Add derivatives arising from a perturbation (if any).
568 if (hasPerturbation()) {
569 perturbation().incrementDc(dc_);
570 }
571
572 hasDc_ = true;
573 }
574
575 /*
576 * Save the current state prior to a next move.
577 *
578 * Invoked before each move.
579 */
580 template <int D>
582 {
583 UTIL_CHECK(system().w().hasData());
584 UTIL_CHECK(hasWc());
585 UTIL_CHECK(state_.isAllocated);
586 UTIL_CHECK(!state_.hasData);
587
588 // Set fields
589 int nMonomer = system().mixture().nMonomer();
590
591 for (int i = 0; i < nMonomer; ++i) {
592 state_.w[i] = system().w().rgrid(i);
593 state_.wc[i] = wc(i);
594 }
595
596 // Save cc based on ccSavePolicy
597 if (state_.needsCc) {
598 UTIL_CHECK(hasCc());
599 UTIL_CHECK(state_.cc.isAllocated());
600 for (int i = 0; i < nMonomer; ++i) {
601 state_.cc[i] = cc(i);
602 }
603 }
604
605 // Save dc based on dcSavePolicy
606 if (state_.needsDc) {
607 UTIL_CHECK(hasDc());
608 UTIL_CHECK(state_.dc.isAllocated());
609 for (int i = 0; i < nMonomer - 1; ++i) {
610 state_.dc[i] = dc(i);
611 }
612 }
613
614 // Save Hamiltonian based on hamiltonianSavePolicy
615 if (state_.needsHamiltonian){
617 state_.hamiltonian = hamiltonian();
618 state_.idealHamiltonian = idealHamiltonian();
619 state_.fieldHamiltonian = fieldHamiltonian();
620 state_.perturbationHamiltonian = perturbationHamiltonian();
621 }
622
623 if (hasPerturbation()) {
624 perturbation().saveState();
625 }
626
627 state_.hasData = true;
628 }
629
630 /*
631 * Restore a saved fts state.
632 *
633 * Invoked after an attempted Monte-Carlo move is rejected
634 * or an fts move fails to converge
635 */
636 template <int D>
638 {
639 UTIL_CHECK(state_.isAllocated);
640 UTIL_CHECK(state_.hasData);
641 const int nMonomer = system().mixture().nMonomer();
642
643 // Restore fields
644 system().w().setRGrid(state_.w);
645
646 // Restore Hamiltonian and components
647 if (state_.needsHamiltonian){
648 hamiltonian_ = state_.hamiltonian;
649 idealHamiltonian_ = state_.idealHamiltonian;
650 fieldHamiltonian_ = state_.fieldHamiltonian;
651 perturbationHamiltonian_ = state_.perturbationHamiltonian;
652 hasHamiltonian_ = true;
653 }
654
655 for (int i = 0; i < nMonomer; ++i) {
656 wc_[i] = state_.wc[i];
657 }
658 hasWc_ = true;
659
660 if (state_.needsCc) {
661 for (int i = 0; i < nMonomer; ++i) {
662 cc_[i] = state_.cc[i];
663 }
664 hasCc_ = true;
665 }
666
667 if (state_.needsDc) {
668 for (int i = 0; i < nMonomer - 1; ++i) {
669 dc_[i] = state_.dc[i];
670 }
671 hasDc_ = true;
672 }
673
674 if (hasPerturbation()) {
675 perturbation().restoreState();
676 }
677
678 state_.hasData = false;
679 }
680
681 /*
682 * Clear the saved Monte-Carlo state.
683 *
684 * Invoked when an attempted Monte-Carlo move is accepted.
685 */
686 template <int D>
688 { state_.hasData = false; }
689
690
691 /*
692 * Output all timer results.
693 */
694 template<int D>
695 void Simulator<D>::outputTimers(std::ostream& out) const
696 {
697 UTIL_CHECK(compressorPtr_);
698 outputMdeCounter(out);
699 compressorPtr_->outputTimers(out);
700 }
701
702 /*
703 * Output modified diffusion equation (MDE) counter.
704 */
705 template<int D>
706 void Simulator<D>::outputMdeCounter(std::ostream& out) const
707 {
708 UTIL_CHECK(compressorPtr_);
709 out << "MDE counter "
710 << compressorPtr_->mdeCounter() << std::endl;
711 out << std::endl;
712 }
713
714 /*
715 * Clear all timers.
716 */
717 template<int D>
719 {
720 UTIL_CHECK(compressorPtr_);
721 compressorPtr_->clearTimers();
722 }
723
724 // Protected Functions
725
726 /*
727 * Optionally read a random number generator seed.
728 */
729 template<int D>
730 void Simulator<D>::readRandomSeed(std::istream& in)
731 {
732 // Optionally read a random number generator seed
733 seed_ = 0;
734 readOptional(in, "seed", seed_);
735
736 // Set random number generator seed
737 // Default value seed_ = 0 uses the clock time.
738 random().setSeed(seed_);
739 }
740
741 /*
742 * Optionally read the Compressor parameter file block.
743 */
744 template<int D>
745 void Simulator<D>::readCompressor(std::istream& in, bool& isEnd)
746 {
747 if (isEnd) return;
748 UTIL_CHECK(compressorFactoryPtr_);
750 std::string className;
751 compressorPtr_ =
752 compressorFactory().readObjectOptional(in, *this,
753 className, isEnd);
754 if (!compressorPtr_ && ParamComponent::echo()) {
755 Log::file() << indent() << " Compressor{ [absent] }\n";
756 }
757 }
758
759 // Functions related to an associated Perturbation
760
761 /*
762 * Optionally read a Perturbation parameter file block.
763 */
764 template<int D>
765 void Simulator<D>::readPerturbation(std::istream& in, bool& isEnd)
766 {
767 if (isEnd) return;
768 UTIL_CHECK(perturbationFactoryPtr_);
770 std::string className;
771 perturbationPtr_ =
772 perturbationFactory().readObjectOptional(in, *this,
773 className, isEnd);
774 if (!perturbationPtr_ && ParamComponent::echo()) {
775 Log::file() << indent() << " Perturbation{ [absent] }\n";
776 }
777 }
778
779 /*
780 * Set the associated Perturbation<D> object.
781 */
782 template<int D>
784 {
785 UTIL_CHECK(ptr);
786 perturbationPtr_ = ptr;
787 }
788
789 // Functions associated with associated Ramp
790
791 /*
792 * Optionally read a parameter file block for an associated Ramp.
793 */
794 template<int D>
795 void Simulator<D>::readRamp(std::istream& in, bool& isEnd)
796 {
797 if (isEnd) return;
798 UTIL_CHECK(!rampPtr_);
799
800 std::string className;
801 rampPtr_ =
802 rampFactory().readObjectOptional(in, *this, className, isEnd);
803 if (!rampPtr_ && ParamComponent::echo()) {
804 Log::file() << indent() << " Ramp{ [absent] }\n";
805 }
806 }
807
808 /*
809 * Set the associated Ramp<D> object.
810 */
811 template<int D>
813 {
814 UTIL_CHECK(ptr);
815 rampPtr_ = ptr;
816 }
817
818}
819}
820#endif
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
Field of real double precision values on an FFT mesh.
Definition cpu/RField.h:29
SolventT & solvent(int id)
Get a solvent solver object.
PolymerT & polymer(int id)
Get a polymer solver object by non-const reference.
Factory for subclasses of Compressor.
Spatial domain for a periodic structure with real fields, on a CPU.
Mesh< D > & mesh()
Get the Mesh by non-const reference.
UnitCell< D > & unitCell()
Get the UnitCell by non-const reference.
Solver and descriptor for a mixture of polymers and solvents.
int nPolymer() const
Get number of polymer species.
int nMonomer() const
Get number of monomer types.
int nSolvent() const
Get number of solvent (point particle) species.
double vMonomer() const
Get monomer reference volume (set to 1.0 by default).
Factory for subclasses of Perturbation.
Base class for additive perturbations of standard FTS Hamiltonian.
Descriptor and solver for one polymer species.
double phi() const
Get the overall volume fraction for this species.
Definition Species.h:149
int nBead() const
Total number of beads in the polymer (bead model).
double mu() const
Get the chemical potential for this species (units kT=1).
Definition Species.h:155
double length() const
Sum of the lengths of all blocks in the polymer (thread model).
Factory for subclasses of Ramp.
Class that varies parameters during a simulation (abstract).
bool hasCc_
Have eigen-components of the current c fields been computed ?
void readPerturbation(std::istream &in, bool &isEnd)
Optionally read an associated perturbation.
Perturbation< D > & perturbation()
Get the perturbation factory by non-const reference.
System< D > & system()
Get parent system by reference.
bool hasWc_
Have eigen-components of the current w fields been computed ?
void computeDc()
Compute functional derivatives of the Hamiltonian.
SimState< D > state_
Previous state saved during at the beginning of a step.
bool hasHamiltonian_
Has the Hamiltonian been computed for the current w and c fields?
virtual void analyze(int min, int max, std::string classname, std::string filename)
Read and analyze a trajectory file.
void setClassName(const char *className)
Set class name string.
double perturbationHamiltonian() const
Get the perturbation to the standard Hamiltonian (if any).
void clearState()
Clear the saved copy of the fts state.
PerturbationFactory< D > & perturbationFactory()
Get the perturbation factory by reference.
void computeWc()
Compute eigenvector components of the current w fields.
DArray< RField< D > > const & cc() const
Get all eigenvector components of the current c fields.
void setPerturbation(Perturbation< D > *ptr)
Set the associated perturbation.
void clearData()
Clear field eigen-components and hamiltonian components.
bool hasDc() const
Are the current d fields valid ?
void analyzeChi()
Perform eigenvalue analysis of projected chi matrix.
double fieldHamiltonian_
Field contribution (H_W) to Hamiltonian.
virtual void outputMdeCounter(std::ostream &out) const
Output MDE counter.
long seed_
Random number generator seed.
virtual void readParameters(std::istream &in)
Read parameters for a simulation.
void allocate()
Allocate required memory.
double idealHamiltonian_
Ideal gas contribution (-lnQ) to Hamiltonian H[W].
virtual void simulate(int nStep)
Perform a field theoretic Monte-Carlo simulation.
void setRamp(Ramp< D > *ptr)
Set the associated ramp.
DArray< RField< D > > wc_
Eigenvector components of w fields on a real space grid.
DArray< RField< D > > dc_
Components of d fields on a real space grid.
bool hasCc() const
Are eigen-components of current c fields valid ?
RampFactory< D > & rampFactory()
Get the ramp factory by reference.
double hamiltonian_
Total field theoretic Hamiltonian H[W] (extensive value).
double perturbationHamiltonian_
Perturbation to the standard Hamiltonian (if any).
virtual void clearTimers()
Clear timers.
bool hasWc() const
Are eigen-components of current w fields valid ?
CompressorFactory< D > & compressorFactory()
Get the compressor factory by reference.
bool hasDc_
Have functional derivatives of H[W] been computed ?
void restoreState()
Restore the saved copy of the fts move state.
void saveState()
Save a copy of the fts move state.
long iTotalStep_
Step counter - total number of attempted BD or MC steps.
double hamiltonian() const
Get the Hamiltonian used in PS-FTS.
bool hasCompressor() const
Does this Simulator have a Compressor?
void readRandomSeed(std::istream &in)
Optionally read a random number generator seed.
bool hasPerturbation() const
Does this Simulator have a Perturbation?
DArray< RField< D > > cc_
Eigenvector components of c fields on a real space grid.
virtual void outputTimers(std::ostream &out) const
Output timing results.
void readCompressor(std::istream &in, bool &isEnd)
Read the compressor block of the parameter file.
double fieldHamiltonian() const
Get the quadratic field contribution to the Hamiltonian.
long iStep_
Step counter - attempted steps for which compressor converges.
void readRamp(std::istream &in, bool &isEnd)
Optionally read an associated ramp.
bool hasHamiltonian() const
Has the Hamiltonian been computed for current w and c fields?
DArray< RField< D > > const & dc() const
Get all of the current d fields.
Random random_
Random number generator.
DArray< RField< D > > const & wc() const
Get all eigenvector components of the current w fields.
void computeHamiltonian()
Compute the Hamiltonian used in PS-FTS.
double idealHamiltonian() const
Get ideal gas contribution to the Hamiltonian.
Random & random()
Get random number generator by reference.
void computeCc()
Compute eigenvector components of the current c fields.
Solver and descriptor for a solvent species.
double phi() const
Get the overall volume fraction for this species.
Definition Species.h:149
double mu() const
Get the chemical potential for this species (units kT=1).
Definition Species.h:155
double size() const
Get the size (number of monomers) in this solvent.
Main class, representing one complete system.
Dynamically allocatable contiguous array template.
Definition DArray.h:32
void allocate(int capacity)
Allocate the underlying C array.
Definition DArray.h:199
Dynamically allocated Matrix.
Definition DMatrix.h:25
void allocate(int capacity1, int capacity2)
Allocate memory for a matrix.
Definition DMatrix.h:170
static std::ostream & file()
Get log ostream by reference.
Definition Log.cpp:59
std::string indent() const
Return indent string for this object (string of spaces).
static bool echo()
Get echo parameter.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
std::string className() const
Get 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
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Definition global.h:49
bool isThread()
Is the thread model in use ?
Real periodic fields, SCFT and PS-FTS (CPU).
Definition param_pc.dox:2
PSCF package top-level namespace.
Definition param_pc.dox:1