File: slr.h

package info (click to toggle)
ergo 3.8-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 17,396 kB
  • sloc: cpp: 94,740; ansic: 17,015; sh: 7,559; makefile: 1,402; yacc: 127; lex: 110; awk: 23
file content (404 lines) | stat: -rw-r--r-- 14,786 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
/* Ergo, version 3.8, a program for linear scaling electronic structure
 * calculations.
 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
 * and Anastasia Kruchinina.
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Primary academic reference:
 * Ergo: An open-source program for linear-scaling electronic structure
 * calculations,
 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
 * Kruchinina,
 * SoftwareX 7, 107 (2018),
 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
 * 
 * For further information about Ergo, see <http://www.ergoscf.org>.
 */

/** @file slr.h

    @brief Contains a Simple Linear Response implementation based on
    the MO orbital picture.

    @author: Pawel Salek <em>responsible</em>
*/

#if !defined(SLR_HEADER)
#define SLR_HEADER
/* Copyright(c) Pawel Salek 2006. */

#include <unistd.h>

#include "realtype.h"

namespace LR {
class VarVector;
/** template based proxy object that uses bool-valued policies to
    perform the assignments.  */
template<bool MultByS,bool SwapXY>
  class VarVectorProxyOp {
 public:
  const VarVector& vec;
  ergo_real scalar;
  explicit VarVectorProxyOp(const VarVector& a, ergo_real s=1.0) : vec(a), scalar(s) {}
};


/** Vector of variables parametrising the solution to the linear
    response equations.  It provides some handy operations for
    compact notation. */
class VarVector {
  ergo_real *v;  /**< vector coefficients */
 public:
  char *fName;    /**< Present in secondary storage (i.e. disk) under
                   * given file name. */
  int  nvar;      /**< nvar := nocc*nvirt. Vector length is 2*nvar */
  unsigned onDisk:1;   /**< valid representation on disk */
  unsigned inMemory:1; /**< valid representation in memory */

  VarVector(const VarVector& a) : fName(NULL), nvar(a.nvar),
    onDisk(0), inMemory(1) {
    if(nvar) {
      v = new ergo_real[nvar*2];
      memcpy(v, a.v, 2*nvar*sizeof(ergo_real));
    } else v = NULL;
  }

  VarVector() : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {}

  /** Creates a vector from a full matrix. */
  VarVector(int nbast, int nocc, const ergo_real *full_mat)
    : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1)  {
    setFromFull(nbast, nocc, full_mat);
  }

  ~VarVector() {
    if(v) delete v; 
    if(fName) {
      unlink(fName);
      delete fName;
    }
}
  ergo_real* x() const { return v; } /**< return the X part of the vector. */
  ergo_real* y() const { return v + nvar; } /**< returns the Y part. */
  void symorth(void);
  void setSize(int n) { /**< sets the size. Reallocates if needed. */
    if(nvar != n) {
      if(v)
	delete v;
      v = new ergo_real[2*n];
      nvar = n;
      onDisk = 0;
    }
  }
  int size() const { return nvar; }
  void print(const char *comment = NULL) {
    if(comment) puts(comment);
    for(int i=0; i<nvar*2; i++) printf("%15.10f\n", (double)v[i]);
  }
  void setFromFull(int nbast, int nocc, const ergo_real *full_mat);
  void setFull(int nbast, int nocc, ergo_real *full_mat) const;
  const ergo_real& operator[](int i) const { return v[i]; }
  ergo_real& operator[](int i)             { return v[i]; }
  VarVector& operator=(ergo_real scalar) {
    for(int i=0; i<2*nvar; i++) v[i] = scalar;
    onDisk = 0;
    return *this;
  };
  VarVector& operator=(const VarVector& b) {
    if(this == &b) return *this;
    if(nvar != b.nvar) { 
      nvar = b.nvar; 
      if(v) delete v;
      v = new ergo_real[2*nvar];
    }
    memcpy(v, b.v, 2*nvar*sizeof(v[0]));
    onDisk = 0;
    return *this;
  }

    VarVector&
    operator=(const VarVectorProxyOp<false, false >& proxy) {
      if (&proxy.vec == this)
        throw "VarVector self-assignment not-implemented";
      if(nvar != proxy.vec.nvar) {
        if(v) delete v; 
        nvar = proxy.vec.nvar;
        v = new ergo_real[2*nvar];
      }

      for(int i=0; i<2*nvar; i++)
        v[i]      = proxy.scalar*proxy.vec[i];
      onDisk = 0;
      return *this;
    }
    VarVector&
    operator=(const VarVectorProxyOp<false, true >& proxy) {
      if (&proxy.vec == this)
        throw "VarVector self-assignment not-implemented";
      if(nvar != proxy.vec.nvar) {
        if(v) delete v; 
        nvar = proxy.vec.nvar;
        v = new ergo_real[2*nvar];
      }
      for(int i=0; i<nvar; i++) {
        v[i]      = proxy.scalar*proxy.vec[i+nvar];
        v[i+nvar] = proxy.scalar*proxy.vec[i];
      }
      onDisk = 0;
      return *this;
    }
    
    /** Load the object to memory. */
    void load(const char* tmpdir);
    /** Save the object. Probably does not need be done explicitly. */
    void save(const char* tmpdir);
    /** Releases the memory, saving if necessary. */
    void release(const char* tmpdir);

    friend ergo_real operator*(const VarVector& a, const VarVector& b);
    friend ergo_real operator*(const VarVector &a,
                               const VarVectorProxyOp<false,false>& b);
    friend ergo_real operator*(const VarVector &a,
                               const VarVectorProxyOp<true,false>& b);
};

/** E2Evaluator interface provides a way to perform a linear
    transformation of supplied transition density matrix @param
    dmat. The result is returned in @param fmat. */
class E2Evaluator {
 public:
  virtual bool transform(const ergo_real *dmat, ergo_real *fmat) = 0;
  virtual ~E2Evaluator() {}
};

/* ################################################################### */
/** a collection of vectors, usually handled at once. */
class VarVectorCollection {
  VarVector *vecs;
  unsigned *ages;
  unsigned currentAge;
  int nVecs;
  int nAllocated;
  bool diskMode;
 public:
  explicit VarVectorCollection(int nSize=0) : vecs(NULL), ages(NULL), currentAge(0),
    nVecs(0), nAllocated(0), diskMode(false) {
    if(nSize) setSize(nSize);
  }
  ~VarVectorCollection();
  void setSize(int sz);
  VarVector& operator[](int i);
  int size() const { return nVecs; }
  bool getDiskMode() const { return diskMode; }
  void setDiskMode(bool x) { diskMode = x; }
  /** Make sure there is space for one vector. */
  void release(); 
  /** Release all vectors from the memory, saving if necessary. */
  void releaseAll();
  static const char *tmpdir;
};
 
/* ################################################################### */
/** Abstract interface to a one electron operator. */
class OneElOperator {
 public:
  virtual void getOper(ergo_real *result) = 0;
  virtual ~OneElOperator() {}
};

/** a class implementing dynamic resized two dimensional arrays. */
class SmallMatrix {
  ergo_real *mat;
  int nsize;
 protected:
  struct RowProxy {
    ergo_real *toprow;
    explicit RowProxy(ergo_real *r) : toprow(r) {}
    ergo_real& operator[](int j) const {
      //printf("  returning row %i -> %p\n", j, toprow + j);
      return toprow[j]; }
  };
 public:
  explicit SmallMatrix(int sz) : mat(new ergo_real[sz*sz]), nsize(sz) {}
  ~SmallMatrix() { delete [] mat; }
  const RowProxy operator[](int i) {
    //printf("Returning column %i -> %p\n", i, mat + i*nsize);
    return RowProxy(mat + i*nsize); }
  void expand(int newSize);
};


/* ################################################################### */
/** Linear Response iterative solver using a variant of the Davidson
    method. */
class LRSolver {
 public:

  LRSolver(int nbast, int nocc,
           const ergo_real *fock_matrix,
           const ergo_real *s);
  virtual ~LRSolver() {/* FIXME: delete esub etc */
    if(cmo)   delete cmo; 
    if(fdiag) delete fdiag;
    delete [] xSub;
  }

  /** Computes the residual vector. The residual vector is created by
      solving the problem in the subspace and then using the solution
      coefficients to form the approximate solution vector. This trial
      vector is then substituted to the equation and the residual is
      defined as the difference between the transformed trial vector
      and the expected solution.  */
  virtual bool getResidual(VarVectorCollection& residualv) = 0;

  /** Computes the initial vector the subspace is to be seeded
      with. Allocates @param vecs and returns the number of
      vectors. */
  virtual int getInitialGuess(VarVectorCollection& vecs) = 0;

  /** returns the preconditioning shift. Proper preconditioning is
      vital for the quick convergence. */
  virtual ergo_real getPreconditionerShift(int i) const = 0;

  /** expands above the default limit */
  virtual void increaseSubspaceLimit(int newSize);

  /** Solves the problem defined by the subclass. */
  bool solve(E2Evaluator& e, bool diskMode = false);
  void computeExactE2Diag(E2Evaluator& e2);
  ergo_real convThreshold;   /**< iterative method convergence threshold */
  int       maxSubspaceSize; /**< current subspace size limit. */
 protected:
  static const int MVEC = 200; /**< default limit for subspace size */
  VarVector e2diag; /**< approximation to the diagonal of E2 operator */
  int subspaceSize; /**< current subspace size */
  SmallMatrix eSub; /**< E[2] matrix projected onto subspace */
  SmallMatrix sSub; /**< S[2] matrix projected onto subspace */
  ergo_real *xSub;  /**< solution vector projected onto subspace */

  /** Computes a vector built of base vectors with specified vectors.
      r := Av - f*Sv */
  void getAvMinusFreqSv(ergo_real f, ergo_real *weights,
                        VarVector& r);

  /** Projects vector @param full on the reduced subspace, returns the
      result in @param w which is a preallocated vector of projection
      coefficients (weights) with size equal at least to the subspace
      size. */
  void projectOnSubspace(const VarVector& full, ergo_real *w)/* const*/;
  /** Build full fector from the reduced form */
  void buildVector(const ergo_real *w, VarVector& full) /* const */;

  /** Transform square operator to the vector form */
  void operToVec(OneElOperator& oper, VarVector& res) const;

  /** setE2diag is called by the constructor to fill in the
   approximation of the E[2] operator diagonal.  It returns
   E_LUMO-E_HOMO which is useful for other things. */
  ergo_real setE2diag(int nbast, int nocc,
                      const ergo_real *fock_matrix,
                      const ergo_real *s);
  
  int nbast; /**< number of basis functions */
  int nocc;  /**< number of occupied orbitals */
  VarVectorCollection vects;  /**< base vectors */
  virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
  void mo2ao(int nbast, const ergo_real *mo, ergo_real *ao) const;
  void ao2mo(int nbast, const ergo_real *ao, ergo_real *mo) const;
 private:
  /** vects and Avects members store the trial vectors and their
      transformed versions. Only every second vector is stored, the
      paired vectors are recovered with help of swapXY() function. */
  VarVectorCollection Avects; /**< transformed base vectors */
  ergo_real *fdiag;       /**< the eigenvalues of the Fock
                             matrix. Used by load_F_MO. */
  ergo_real *cmo;         /**< the MO coefficients. */

  void load_F_MO(ergo_real *fmat) const;
  bool lintrans(E2Evaluator& e2, const VarVector& v, VarVector& Av) const;
};

/* ################################################################### */
/** Iterative Set Of Linear Equations solver, extending the generic
    LRSolver. */
class SetOfEqSolver : public LRSolver {
  ergo_real frequency; /**< frequency for which the SOE is to be solved. */
  VarVector rhs; /**< RHS of the SOE */
 public:
  /** Creates the set-of-equations solver. The KS and overlap matrix
      may be deleted immediately after the object creation. */
  SetOfEqSolver(int nbast, int nocc,
                const ergo_real *fock_matrix,
                const ergo_real *s,
		ergo_real freq)
    : LRSolver(nbast, nocc, fock_matrix, s), frequency(freq),
    rhsSub(new ergo_real[MVEC]) {};
  void setRHS(OneElOperator& op);    /**< initializes the rhs field */
  virtual ~SetOfEqSolver() { delete [] rhsSub; }
  virtual ergo_real getPreconditionerShift(int) const { return frequency; }
  virtual int getInitialGuess(VarVectorCollection& vecs);
  virtual bool getResidual(VarVectorCollection& residualv);
  /** expands above the default limit */
  virtual void increaseSubspaceLimit(int newSize);
  ergo_real getPolarisability(OneElOperator& oper) /* const */;
 protected:
  ergo_real *rhsSub; /**< RHS vector projected onto subspace */
  virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
  ergo_real multiplyXtimesVec(const VarVector& rhs) /* const */;
  ergo_real xTimesRHS;
};


/* ################################################################### */
/** Iterative Eigenvalue solver, extending the generic LRSolver. */
class EigenSolver : public LRSolver {
  ergo_real *ritzVals; /**< recent ritz values in the subspace. */
  ergo_real *transMoms2; /**< most recent SQUARED transition moments. */
  int nStates;    /**< number of excited states to compute */
  int nConverged; /**< number of already converged eigenvalues */
  ergo_real *last_ev; /**< most recent eigenvectors in the reduced space */
 public:
  EigenSolver(int nbast, int nocc,
              const ergo_real *fock_matrix,
              const ergo_real *overlap, int n)
    : LRSolver(nbast, nocc, NULL, NULL), 
    ritzVals(new ergo_real[MVEC]), transMoms2(new ergo_real[MVEC]),
    nStates(n), nConverged(0), last_ev(NULL) {
    ritzVals[0] = setE2diag(nbast, nocc, fock_matrix, overlap);
    for(int i=1; i<n; i++) ritzVals[i] = ritzVals[0];
  }
  virtual ~EigenSolver() {
    if(last_ev) delete last_ev;
    delete [] ritzVals;
    delete [] transMoms2;
  }
  virtual ergo_real getPreconditionerShift(int i) const { 
    return ritzVals[nConverged+i];
  }
  virtual int getInitialGuess(VarVectorCollection& vecs);
  virtual bool getResidual(VarVectorCollection& residualv);
  /** expands above the default limit */
  virtual void increaseSubspaceLimit(int newSize);

  ergo_real getFreq(int i) const { return ritzVals[i]; }
  void computeMoments(OneElOperator& dipx,
                      OneElOperator& dipy,
                      OneElOperator& dipz);
  ergo_real getTransitionMoment2(int i) const { return transMoms2[i]; }
};

} /* End of namespace LR */
#endif /* !defined(SLR_HEADER) */