JDFTx  1.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Pulay.h
Go to the documentation of this file.
1 /*-------------------------------------------------------------------
2 Copyright 2015 Ravishankar Sundararaman
3 
4 This file is part of JDFTx.
5 
6 JDFTx is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10 
11 JDFTx is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 
16 You should have received a copy of the GNU General Public License
17 along with JDFTx. If not, see <http://www.gnu.org/licenses/>.
18 -------------------------------------------------------------------*/
19 
20 #ifndef JDFTX_CORE_PULAY_H
21 #define JDFTX_CORE_PULAY_H
22 
23 #include <core/PulayParams.h>
24 #include <core/Util.h>
25 #include <electronic/matrix.h>
26 #include <cfloat>
27 
33 template<typename Variable> class Pulay
36 {
37 public:
38  Pulay(const PulayParams& pp);
39 
45  double minimize(double Eprev=+DBL_MAX, std::vector<string> extraNames=std::vector<string>(), std::vector<double> extraThresh=std::vector<double>());
46 
47  void loadState(const char* filename);
48  void saveState(const char* filename) const;
49  void clearState();
50 
52  virtual double sync(double x) const { return x; }
53 
54 protected:
55  //----- Interface specification -----
56 
65  virtual double cycle(double dEprev, std::vector<double>& extraValues)=0;
66 
67  virtual void report(int iter) {}
68  virtual void axpy(double alpha, const Variable& X, Variable& Y) const=0;
69  virtual double dot(const Variable& X, const Variable& Y) const=0;
70  virtual size_t variableSize() const=0;
71  virtual void readVariable(Variable&, FILE*) const=0;
72  virtual void writeVariable(const Variable&, FILE*) const=0;
73  virtual Variable getVariable() const=0;
74  virtual void setVariable(const Variable&)=0;
75  virtual Variable precondition(const Variable&) const=0;
76  virtual Variable applyMetric(const Variable&) const=0;
77 
78 private:
79  const PulayParams& pp;
80  std::vector<Variable> pastVariables;
81  std::vector<Variable> pastResiduals;
82  matrix overlap;
83 };
84 
85 
86 //---------------------- Implementation ----------------------------
88 
89 #include <core/Minimize_linmin.h>
90 #include <memory>
91 
92 //Norm convergence check (eigenvalue-difference or residual)
93 //Make sure value is within tolerance for nCheck consecutive cycles
94 class NormCheck
95 { unsigned nCheck; double threshold;
96  std::deque<bool> history;
97 public:
98  NormCheck(unsigned nCheck, double threshold) : nCheck(nCheck), threshold(fabs(threshold)) {}
99  bool checkConvergence(double norm)
100  { history.push_back(fabs(norm)<threshold);
101  if(history.size()==nCheck+1) history.pop_front(); //discard old unneeded elements
102  if(history.size()==nCheck)
103  { for(bool converged: history)
104  if(!converged)
105  return false;
106  return true;
107  }
108  else return false;
109  }
110 };
111 
112 template<typename Variable> Pulay<Variable>::Pulay(const PulayParams& pp)
113 : pp(pp), overlap(pp.history, pp.history)
114 {
115 }
116 
117 template<typename Variable> double Pulay<Variable>::minimize(double Eprev, std::vector<string> extraNames, std::vector<double> extraThresh)
118 {
119  double E = sync(Eprev); Eprev = 0.;
120  double dE = E-Eprev;
121  assert(extraNames.size()==extraThresh.size());
122 
123  //Initialize convergence checkers:
124  EdiffCheck ediffCheck(2, pp.energyDiffThreshold); ediffCheck.checkConvergence(E); //store the initial energy in the check's history
125  NormCheck resCheck(2, pp.residualThreshold);
126  std::vector<std::shared_ptr<NormCheck> > extraCheck(extraNames.size());
127  for(size_t iExtra=0; iExtra<extraNames.size(); iExtra++)
128  extraCheck[iExtra] = std::make_shared<NormCheck>(2, extraThresh[iExtra]);
129 
130  for(int iter=0; iter<pp.nIterations; iter++)
131  {
132  //If history is full, remove oldest member
133  assert(pastResiduals.size() == pastVariables.size());
134  if((int)pastResiduals.size() >= pp.history)
135  { size_t ndim = pastResiduals.size();
136  if(ndim>1) overlap.set(0,ndim-1, 0,ndim-1, overlap(1,ndim, 1,ndim));
137  pastVariables.erase(pastVariables.begin());
138  pastResiduals.erase(pastResiduals.begin());
139  }
140 
141  //Cache the old energy and variables
142  Eprev = E;
143  pastVariables.push_back(getVariable());
144 
145  //Perform cycle:
146  std::vector<double> extraValues(extraThresh.size());
147  E = sync(cycle(dE, extraValues));
148  dE = E - Eprev;
149  for(auto& v: extraValues) v = sync(v);
150 
151  //Calculate and cache residual:
152  double residualNorm = 0.;
153  { Variable residual = getVariable(); axpy(-1., pastVariables.back(), residual);
154  pastResiduals.push_back(residual);
155  residualNorm = sync(sqrt(dot(residual,residual)));
156  }
157 
158  //Print energy and convergence parameters:
159  fprintf(pp.fpLog, "%sCycle: %2i %s: ", pp.linePrefix, iter, pp.energyLabel);
160  fprintf(pp.fpLog, pp.energyFormat, E);
161  fprintf(pp.fpLog, " d%s: %+.3e", pp.energyLabel, dE);
162  fprintf(pp.fpLog, " |Residual|: %.3e", residualNorm);
163  for(size_t iExtra=0; iExtra<extraNames.size(); iExtra++)
164  fprintf(pp.fpLog, " |%s|: %.3e", extraNames[iExtra].c_str(), extraValues[iExtra]);
165  fprintf(pp.fpLog, " t[s]: %9.2lf", clock_sec());
166  fprintf(pp.fpLog, "\n"); fflush(pp.fpLog);
167 
168  //Optional reporting:
169  report(iter);
170 
171  //Check for convergence and update variable:
172  if(std::isnan(E))
173  { fprintf(pp.fpLog, "%sE=%le. Stopping ...\n\n", pp.linePrefix, E);
174  return E;
175  }
176  bool converged = false;
177  if(!converged && ediffCheck.checkConvergence(E))
178  { fprintf(pp.fpLog, "%sConverged (|Delta E|<%le for 2 iters).\n\n", pp.linePrefix, pp.energyDiffThreshold);
179  converged = true;
180  }
181  if(!converged && resCheck.checkConvergence(residualNorm))
182  { fprintf(pp.fpLog, "%sConverged (|Residual|<%le for 2 iters).\n\n", pp.linePrefix, pp.residualThreshold);
183  converged = true;
184  }
185  if(!converged && extraNames.size())
186  for(size_t iExtra=0; iExtra<extraNames.size(); iExtra++)
187  if(extraCheck[iExtra]->checkConvergence(extraValues[iExtra]))
188  { fprintf(pp.fpLog, "%sConverged (|%s|<%le for 2 iters).\n\n", pp.linePrefix, extraNames[iExtra].c_str(), extraThresh[iExtra]);
189  converged = true;
190  break;
191  }
192  fflush(pp.fpLog);
193  if(converged || killFlag) break; //converged or manually interrupted
194 
195  //---- DIIS/Pulay mixing -----
196 
197  //Update the overlap matrix
198  size_t ndim = pastResiduals.size();
199  Variable MlastResidual = applyMetric(pastResiduals.back());
200  for(size_t j=0; j<ndim; j++)
201  { double thisOverlap = dot(pastResiduals[j], MlastResidual);
202  overlap.set(j, ndim-1, thisOverlap);
203  overlap.set(ndim-1, j, thisOverlap);
204  }
205 
206  //Invert the residual overlap matrix to get the minimum of residual
207  matrix cOverlap(ndim+1, ndim+1); //Add row and column to enforce normalization constraint
208  cOverlap.set(0, ndim, 0, ndim, overlap(0, ndim, 0, ndim));
209  for(size_t j=0; j<ndim; j++)
210  { cOverlap.set(j, ndim, 1);
211  cOverlap.set(ndim, j, 1);
212  }
213  cOverlap.set(ndim, ndim, 0);
214  matrix cOverlap_inv = inv(cOverlap);
215 
216  //Update variable:
217  complex* coefs = cOverlap_inv.data();
218  Variable v;
219  for(size_t j=0; j<ndim; j++)
220  { double alpha = coefs[cOverlap_inv.index(j, ndim)].real();
221  axpy(alpha, pastVariables[j], v);
222  axpy(alpha, precondition(pastResiduals[j]), v);
223  }
224  setVariable(v);
225  }
226  return E;
227 }
228 
229 template<typename Variable> void Pulay<Variable>::loadState(const char* filename)
230 {
231  size_t nBytesCycle = 2 * variableSize(); //number of bytes per history entry
232  size_t nBytesFile = fileSize(filename);
233  size_t ndim = nBytesFile / nBytesCycle;
234  size_t dimOffset = 0;
235  if(int(ndim) > pp.history)
236  { dimOffset = ndim - pp.history;
237  ndim = pp.history;
238  }
239  if(nBytesFile % nBytesCycle != 0)
240  die("Pulay history file '%s' does not contain an integral multiple of the mixed variables and residuals.\n", filename);
241  fprintf(pp.fpLog, "%sReading %lu past variables and residuals from '%s' ... ", pp.linePrefix, ndim, filename); logFlush();
242  pastVariables.resize(ndim);
243  pastResiduals.resize(ndim);
244  FILE* fp = fopen(filename, "r");
245  if(dimOffset) fseek(fp, dimOffset*nBytesCycle, SEEK_SET);
246  for(size_t idim=0; idim<ndim; idim++)
247  { readVariable(pastVariables[idim], fp);
248  readVariable(pastResiduals[idim], fp);
249  }
250  fclose(fp);
251  fprintf(pp.fpLog, "done.\n"); fflush(pp.fpLog);
252  //Compute overlaps of loaded history:
253  for(size_t i=0; i<ndim; i++)
254  { Variable Mresidual_i = applyMetric(pastResiduals[i]);
255  for(size_t j=0; j<=i; j++)
256  { double thisOverlap = dot(pastResiduals[j], Mresidual_i);
257  overlap.set(i,j, thisOverlap);
258  overlap.set(j,i, thisOverlap);
259  }
260  }
261 }
262 
263 template<typename Variable> void Pulay<Variable>::saveState(const char* filename) const
264 {
265  if(mpiUtil->isHead())
266  { FILE* fp = fopen(filename, "w");
267  for(size_t idim=0; idim<pastVariables.size(); idim++)
268  { writeVariable(pastVariables[idim], fp);
269  writeVariable(pastResiduals[idim], fp);
270  }
271  fclose(fp);
272  }
273 }
274 
275 template<typename Variable> void Pulay<Variable>::clearState()
276 { pastVariables.clear();
277  pastResiduals.clear();
278 }
279 
281 
282 #endif //JDFTX_CORE_PULAY_H
Pulay mixing to optimize self-consistent field optimization
Definition: Pulay.h:35
Definition: Minimize_linmin.h:27
ScalarField sqrt(const ScalarField &)
Elementwise square root (preserve input)
virtual Variable applyMetric(const Variable &) const =0
Apply metric to variable/residual.
complex dot(const Tptr &X, const Tptr &Y)
Definition: Operators.h:195
General complex matrix.
Definition: matrix.h:57
#define logFlush()
fflush() for log files
Definition: Util.h:111
virtual double dot(const Variable &X, const Variable &Y) const =0
Euclidean dot product. Metric applied separately for efficiency.
bool isHead() const
whether this is the root process (makes code more readable)
Definition: MPIUtil.h:39
int index(int i, int j) const
Index into data()
Definition: matrix.h:65
off_t fileSize(const char *filename)
Get the size of a file.
virtual Variable precondition(const Variable &) const =0
Apply preconditioner to variable/residual.
virtual void setVariable(const Variable &)=0
Set the state of system to specified variable.
virtual void readVariable(Variable &, FILE *) const =0
Read variable from stream.
Parameters to control Pulay mixing.
Definition: PulayParams.h:26
virtual void report(int iter)
Override to perform optional reporting.
Definition: Pulay.h:67
double clock_sec()
Elapsed time in microseconds (from start of program)
void clearState()
remove past variables and residuals
bool killFlag
Flag set by signal handlers - all compute loops should quit cleanly when this is set.
void loadState(const char *filename)
Load the state from a single binary file.
#define die(...)
Quit with an error message (formatted using printf()). Must be called from all processes.
Definition: Util.h:114
Miscellaneous utilities.
ScalarField inv(const ScalarField &)
Elementwise reciprocal (preserve input)
void saveState(const char *filename) const
Save the state to a single binary file.
virtual Variable getVariable() const =0
Write variable to stream.
virtual double cycle(double dEprev, std::vector< double > &extraValues)=0
Single cycle of the self-consistency loop. In each subsequent cycle, Pulay will try to zero the diffe...
complex * data()
Return a pointer to the actual data Return a CPU pointer to the actual data, will move data from GPU ...
virtual double sync(double x) const
Override to synchronize scalars over MPI processes (if the same minimization is happening in sync ove...
Definition: Pulay.h:52
Complex number (need to define our own because we need operators for gpu code as well) ...
Definition: scalar.h:49
double minimize(double Eprev=+DBL_MAX, std::vector< string > extraNames=std::vector< string >(), std::vector< double > extraThresh=std::vector< double >())
Minimize energy using a self-consistent iteration.
void axpy(double alpha, const Tptr &X, Tptr &Y)
Generic axpy for complex data types (Note: null pointers are treated as zero)
Definition: Operators.h:157
virtual void axpy(double alpha, const Variable &X, Variable &Y) const =0
Scaled accumulate on variable.
virtual size_t variableSize() const =0
Number of bytes per variable.
#define assert(expr)
A custom assertion with stack trace (NOTE: enabled in release modes as well)
Definition: Util.h:100