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