53 read(in,
"mobility", mobility_);
56 int nMonomer = system().mixture().nMonomer();
57 IntVec<D> meshDimensions = system().domain().mesh().dimensions();
58 wp_.allocate(nMonomer);
59 wf_.allocate(nMonomer);
60 for (
int i=0; i < nMonomer; ++i) {
61 wp_[i].allocate(meshDimensions);
62 wf_[i].allocate(meshDimensions);
64 dci_.allocate(nMonomer-1);
65 eta_.allocate(nMonomer-1);
66 for (
int i=0; i < nMonomer - 1; ++i) {
67 dci_[i].allocate(meshDimensions);
68 eta_[i].allocate(meshDimensions);
70 dwc_.allocate(meshDimensions);
71 dwp_.allocate(meshDimensions);
78 int meshSize = system().domain().mesh().size();
79 int nMonomer = system().mixture().nMonomer();
82 for (
int i=0; i < nMonomer; ++i) {
88 for (
int i=0; i < nMonomer - 1; ++i) {
100 const int nMonomer = system().mixture().nMonomer();
101 const int meshSize = system().domain().mesh().size();
105 simulator().saveState();
108 for (i = 0; i < nMonomer; ++i) {
109 wp_[i] = system().w().rgrid(i);
114 dwp_ = simulator().wc(nMonomer-1);
117 const double vSystem = system().domain().unitCell().volume();
118 const double a = -1.0*mobility_;
119 const double b = sqrt(2.0*mobility_*
double(meshSize)/vSystem);
122 for (j = 0; j < nMonomer - 1; ++j) {
124 for (k = 0; k < meshSize; ++k) {
125 eta[k] = b*random().gaussian();
133 for (j = 0; j < nMonomer - 1; ++j) {
134 RField<D> const & dc = simulator().dc(j);
137 for (k = 0; k < meshSize; ++k) {
138 dwc_[k] = a*dc[k] + eta[k];
142 for (i = 0; i < nMonomer; ++i) {
144 evec = simulator().chiEvecs(j,i);
145 for (k = 0; k < meshSize; ++k) {
146 wp[k] += evec*dwc_[k];
152 system().setWRGrid(wp_);
155 bool isConverged =
false;
158 int compress = simulator().compressor().compress();
160 simulator().restoreState();
165 simulator().clearData();
166 simulator().computeWc();
167 simulator().computeCc();
168 simulator().computeDc();
172 RField<D> const & wp = simulator().wc(nMonomer-1);
173 for (k = 0; k < meshSize; ++k) {
174 dwp_[k] = wp[k] - dwp_[k];
178 for (i = 0; i < nMonomer; ++i) {
180 for (k = 0; k < meshSize; ++k) {
186 const double ha = 0.5*a;
187 for (j = 0; j < nMonomer - 1; ++j) {
188 RField<D> const & dcp = simulator().dc(j);
191 for (k = 0; k < meshSize; ++k) {
192 dwc_[k] = ha*( dci[k] + dcp[k]) + eta[k];
194 for (i = 0; i < nMonomer; ++i) {
196 evec = simulator().chiEvecs(j,i);
197 for (k = 0; k < meshSize; ++k) {
198 wf[k] += evec*dwc_[k];
204 system().setWRGrid(wf_);
207 int compress2 = simulator().compressor().compress();
209 simulator().restoreState();
215 simulator().clearState();
216 simulator().clearData();
217 simulator().computeWc();
218 simulator().computeCc();
219 simulator().computeDc();