1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program 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  * This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
30 /** @file dft_common.h
31 
32     @brief Common DFT routines. Mostly functional mixing.
33 
34     @author: Pawel Salek <em>responsible</em>
35 */
36 
37 #ifndef _DFT_COMMON_H_
38 #define _DFT_COMMON_H_
39 
40 #include <stdlib.h>
41 #include <vector>
42 
43 #ifdef __cplusplus
44 #define EXTERN_C extern "C"
45 #else
46 #define EXTERN_C
47 #endif
48 
49 #include "realtype.h"
50 #include "basisinfo.h"
51 #include "matrix_typedefs.h"
52 #include "functionals.h"
53 #include "grid_atomic.h"
54 
55 /** A vector of first order derivatives with respect to two
56  * parameters: density rho and SQUARE of the gradient of density grho.
57  * zeta_i = |nabla rho_i|^2
58  */
59 
60 typedef struct {
61     real fR;  /* d/drho F     */
62     real fZ;  /* d/zeta F     */
63 } FirstDrv;
64 
65 /* SecondDrv:  matrix  of  second  order functional  derivatives  with
66  * respect  to two  parameters: density  rho and  SQUARE  of the
67  * density gradient zeta.  The derivatives are computed for alpha-alpha
68  * or beta-beta spin-orbital block (i.e. include triplet flag).
69  */
70 typedef struct {
71     real fR; /* d/drho  F */
72     real fZ; /* d/dzeta F */
73     real fRR; /* d/drho^2 F */
74     real fRZ; /* d/(drho dzeta) F */
75     real fZZ; /* d/dzeta^2 F */
76     /* additional derivatives required by  */
77     /* general linear response scheme     */
78     real fRG; /* d/(drho dgamma) F */
79     real fZG; /* d/(dzeta dgamma) F */
80     real fGG; /* d/dgamma^2 F */
81     real fG;  /* d/dgamma F */
82 } SecondDrv;
83 
84 
85 EXTERN_C void dftpot0_(FirstDrv *ds, const real* weight, const FunDensProp* dp);
86 EXTERN_C void dftpot1_(SecondDrv *ds, const real* w, const FunDensProp* dp,
87 		       const int* triplet);
88 
89 EXTERN_C void dft_init(void);
90 EXTERN_C int dft_setfunc(const char *line);
91 
92 class ShellTree;
93 
94 /** Ergo specific implementation of molecule-grid interface. */
95 class ErgoMolInfo : public GridGenMolInfo {
96   const BasisInfoStruct& bis;
97   const Molecule&        molecule;
98  public:
99   ErgoMolInfo(const BasisInfoStruct& bis_,  const Molecule& mol);
100   virtual ~ErgoMolInfo();
101 
102   virtual void getAtom(int icent, int *cnt, real (*coor)[3],
103                        int *charge, int *mult) const;
104   virtual void setShellRadii(real *shellRadii) const;
105   virtual void getBlocks(const real *center, real cellsz,
106                          const real *rshell,
107                          int *nblcnt, int (*iblcks)[2]) const;
108   void getBlocks1(const real *center, real cellsz,
109                   const real *rshell,
110                   int *nblcnt, int (*iblcks)[2]) const;
111   virtual void getExps(int *maxl, int **nucbas, real (**aa)[2]) const;
112   ShellTree *shellTree;
113 };
114 
115 EXTERN_C void ergoShellsToOrbs(const int *nshlbl, const int (*shlblock)[2],
116                                int *norbbl, int (*orbblock)[2],
117                                const BasisInfoStruct& bis);
118 
119 EXTERN_C int dft_get_num_threads();
120 EXTERN_C void dft_set_num_threads(int nThreads);
121 
122 
123 EXTERN_C void dft_init(void);
124 
125 #define dal_new(sz,tp) (tp*)dal_malloc_((sz)*sizeof(tp),__FUNCTION__, __LINE__)
126 void* dal_malloc_(size_t sz, const char *func, unsigned line);
127 /* dal_malloc: usage discouraged */
128 #define dal_malloc(sz) dal_malloc_((sz),__FUNCTION__, __LINE__)
129 
130 /* useful  constants for BLAS interfacing */
131 extern int  ZEROI, ONEI, THREEI, FOURI;
132 extern real ZEROR, ONER, TWOR, FOURR;
133 
134 /** Class Box provides an ability to determine box containing all
135     Objects.  The class Object must provide field center[] and method
136     radius(). */
137 class Box {
138 public:
139   real getDistanceTo(const real* v) const;
140   int getMaxDim() const;
size(int dim)141   real size(int dim) const { return hi[dim]-lo[dim]; }
142 
overlapsWith(const real * center,real radius)143   bool overlapsWith(const real *center, real radius) const {
144     real d = getDistanceTo(center);
145     return d < radius;
146   }
147 
148   /** Determines whether given point is inside the box. In order to
149       avoid double counting, the points that are overlap with the
150       lower limits are included but those that overlap with the higher
151       limit are excluded. */
contains(const real * p)152   bool contains(const real *p) const {
153 #if 0
154     printf("B:(%8.2f %8.2f %8.2f)-(%8.2f %8.2f %8.2f): %8.2f %8.2f %8.2f ",
155            lo[0], lo[1], lo[2], hi[0], hi[1], hi[2],
156            p[0], p[1], p[2]);
157 #endif
158     for(int i=0; i<3; i++)
159       if(p[i]<lo[i] || p[i] >= hi[i]) {
160         //puts("F");
161         return false;
162       }
163     //puts("T");
164     return true;
165   }
166 
167   real lo[3];
168   real hi[3];
169 };
170 
171 template<typename Iterator>
getBoundingBox(Box & box,Iterator start,Iterator end)172   void getBoundingBox(Box& box, Iterator start, Iterator end)
173 {
174   static const ergo_real OFF = 0.1;
175   if(start == end)
176     throw "BoundingBox called for empty set";
177 
178   real r = start->radius() + OFF;
179   for(int i=0; i<3; i++) {
180     box.lo[i] = start->center[i]-r;
181     box.hi[i] = start->center[i]+r;
182   }
183 
184   for(++start; start != end; ++start) {
185     real r = start->radius() + OFF;
186     for(int i=0; i<3; i++) {
187       real l = start->center[i]-r; if (l<box.lo[i]) box.lo[i] = l;
188       real h = start->center[i]+r; if (h>box.hi[i]) box.hi[i] = h;
189     }
190   }
191 }
192 
193 
194 int sync_threads(bool release, int nThreads);
195 
196 #endif /* _DFT_COMMON_H_ */
197