1 /* Copyright (c) 2015  Gerald Knizia
2  *
3  * This file is part of the IboView program (see: http://www.iboview.org)
4  *
5  * IboView is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, version 3.
8  *
9  * IboView is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with bfint (LICENSE). If not, see http://www.gnu.org/licenses/
16  *
17  * Please see IboView documentation in README.txt for:
18  * -- A list of included external software and their licenses. The included
19  *    external software's copyright is not touched by this agreement.
20  * -- Notes on re-distribution and contributions to/further development of
21  *    the IboView software
22  */
23 
24 #ifndef CT8K_BASISSET_H
25 #define CT8K_BASISSET_H
26 
27 #include <vector>
28 #include "Ir.h"
29 #include "CxPodArray.h"
30 #include "CtBasisShell.h"
31 
32 namespace ct {
33 
34 enum FBasisContext {
35    BASIS_Orbital,
36    BASIS_JFit,
37    BASIS_JkFit,
38    BASIS_Mp2Fit,
39    BASIS_CcsdFit,
40    BASIS_F12RI,
41    BASIS_Guess,
42    BASIS_IaoFit,
43    BASIS_Multipoles
44 };
45 
46 std::string BasisContextName( FBasisContext Context );
47 
48 struct FMatrixView;
49 struct FAtomSet;
50 
51 // IR driver-level interface to a basis set. Semi-lightweight.
52 // supposed to be generated from the program's native basis set
53 // format (FBasisSet here) on demand for interfacing between IR
54 // and the program.
55 struct FRawBasis : public FIntrusivePtrDest
56 {
57    typedef TArray<ir::FRawShell>
58       FShellArray;
59    FShellArray
60       Shells;
61    TArray<size_t>
62       // offset of Shell[i] with respect to underlying basis in terms of
63       // individual functions. Contains one additional element giving the
64       // total number of functions.
65       ShellOffsets,
66       // shells of atom #i start at CenterOffsets[i].
67       CenterOffsets,
68       // [i]: center index of shell [i].
69       ShellCenters;
70    size_t
71       // largest number of functions in a shell occuring in the basis.
72       // May be useful for allocating intermediates.
73       nMaxFnPerShell;
nFnFRawBasis74    size_t nFn() const { return ShellOffsets.back(); }
nShFRawBasis75    size_t nSh() const { return Shells.size(); }
nFnFRawBasis76    size_t nFn(size_t iSh) const { return ShellOffsets[iSh+1]-ShellOffsets[iSh]; }
iFnFRawBasis77    size_t iFn(size_t iSh) const { return ShellOffsets[iSh]; }
78    // index of first shell on center iCen
iCenShFRawBasis79    size_t iCenSh(size_t iCen) const { return CenterOffsets[iCen]; }
nCenShFRawBasis80    size_t nCenSh(size_t iCen) const { return CenterOffsets[iCen+1] - CenterOffsets[iCen]; }
81    // index of first function center iCen
iCenFnFRawBasis82    size_t iCenFn(size_t iCen) const { return iFn(iCenSh(iCen)); }
nCenFnFRawBasis83    size_t nCenFn(size_t iCen) const { return iCenFn(iCen+1) - iCenFn(iCen); }
84    // total number of centers
nCenFRawBasis85    size_t nCen() const { return CenterOffsets.size() - 1; }
86 
87    // center index of a given shell.
iShCenFRawBasis88    size_t iShCen(size_t iSh) const { return ShellCenters[iSh]; }
89 
90    // construct derived data (ShellOffsets, CenterOffsets, nMaxFnPerShell) from
91    // this->Shells and the shell center indices given as argument. pShCen[iSh]
92    // should give the center id for Shells[iSh], and those must be consecutive.
93    void Finalize(size_t const *pShCen);
94 };
95 typedef boost::intrusive_ptr<FRawBasis>
96    FRawBasisPtr;
97 
98 // program native format of a basis set.
99 struct FBasisSet : FIntrusivePtrDest
100 {
101    typedef std::vector<FBasisShell>
102       FBasisShellArray;
103    FBasisShellArray
104       Shells;  // Gaussian shells centered on different atoms.
105    FBasisContext
106       Context; // for informative purposes only.
107    std::string
108       Name; // for informative purposes only.
109 
110    void Print(std::ostream &out) const;
111    size_t nFn() const;
112    size_t nPrimitiveGtos() const; // for decoration purposes only.
113    size_t nFnOfLargestShell() const;
114    uint nMaxL() const; // find largest angular momentum in the basis.
115 
116    // rotate and/or translate the basis (i.e., the positions of the atoms) in real space. R is a
117    // 3x3 matrix and 'd' is a 3-vector such that r' = R (r - d) gives the transformed atomic positions.
118    void Transform_3x4(double const *R, double const *d, FMemoryStack &Mem);
119 
120    // transform a set of MO coefficients (in C) in such a way that if the molecule and its
121    // basis are rotated in real space by the 3x3 matrix R (via FBasisSet Transform_3x4), then
122    // this function computes the MO coefficients corresponding to the new alignment.
123    void TransformMoCoeffs_3x4(FMatrixView C, double const *R, FMemoryStack &Mem) const;
124 
125    // construct from atom set and context
126    FBasisSet(FAtomSet const &AtomSet, FBasisContext Context_);
127    // construct from explicit list of shells. Name and Context for informative purposes only.
128    FBasisSet(FBasisShell const *pShells_, size_t nShells_, FBasisContext Context_, std::string const &Name);
129 public:
130    FRawBasisPtr
131       // Only valid if RebuildInterfacingData has been called after the basis
132       // was last modified. If constructed regularly, this is done
133       // automatically.
134       pRawBasis;
135 
136    // rebuilds the data for low-level interfaces (i.e., the data in this section)
137    // from this->Shells
138    void RebuildInterfacingData();
139 
140    FRawBasisPtr MakeRawBasis();
141 private:
142    // makes *this a basis set corresponding to the basis descriptions in
143    // AtomSet. Attempts to load the neccessary external data. Throws
144    // std::runtime_error if failed.
145    void LoadFromAtomSet( FAtomSet const &AtomSet, FBasisContext Context );
146 
147    void MakeAtomOffsets( size_t *&pAtomShellOffsets, size_t *&pAtomBfnOffsets, size_t nAtoms_, FMemoryStack &Mem ) const;
148    void MakeShellOffsets( size_t *&pShellOffsets, FMemoryStack &Mem ) const;
149    void Finalize();
150 };
151 
152 void Trafo3x4(FVector3 &InOut, double const *R, double const *d);
153 
154 
155 inline std::ostream &operator << ( std::ostream &out, FBasisSet const &BasisSet ) {
156    BasisSet.Print(out);
157    return out;
158 }
159 
160 typedef boost::intrusive_ptr<FBasisSet>
161    FBasisSetPtr;
162 
163 // generalized 2-index integral kernel: effectively a functor which
164 // can calculate shell doublets for two shell objects. Used for making
165 // one-electron matrices.
166 struct FKrn2i : public FIntrusivePtrDest
167 {
168    virtual void EvalInt2e2c(double *pOut, size_t StrideA, size_t StrideC, ir::FRawShell const *pA, ir::FRawShell const *pC, double Prefactor, bool Add, ct::FMemoryStack &Mem) const = 0;
169    // Accumulate derivative contraction
170    //     Grad[Nuc,xyz] += Prefactor * \sum_{ac} rdm[a,c] d/d[Nuc,xyz] (a|krn|c).
171    // Parameters:
172    //   - pGrad[3*iCenter + ixyz]: Location in which the output gradient for center iCenter and ixyz=0...2 is stored.
173    //   - pRdmAC[StrideA * ia + StrideC * ic]: gradient coefficients for *pA, *pC,
174    //     where ia/ic go over the contracted functions of *pA, *pC.
175    //   - iCenA, iCenC: center indices in pGrad itself which belong
176    // default implementation just crashes and says "can't do this".
177    virtual void AccCoreDeriv1d(double *pGrad, double const *pRdmAC, size_t StrideA, size_t StrideC, ir::FRawShell const *pA, size_t iCenterA, ir::FRawShell const *pC, size_t iCenterC, double Prefactor, ct::FMemoryStack &Mem) const;
178 };
179 typedef boost::intrusive_ptr<FKrn2i>
180    FKrn2iPtr;
181 
182 void MakeIntMatrix( FMatrixView &Out, FBasisSet const &RowBasis,
183    FBasisSet const &ColBasis, FKrn2i const &Krn2i, FMemoryStack &Mem, double Prefactor=1., bool Add=false);
184 void MakeIntMatrix( FMatrixView &Out, FRawBasis const &RowBasis,
185    FRawBasis const &ColBasis, FKrn2i const &Krn2i, FMemoryStack &Mem, double Prefactor=1., bool Add=false);
186 // increment gradient in pGrad[ixyz + 3 *iCenter] by \sum{a,b} Rdm[a,b] d/d[Eu] Krn2i[a,b].
187 // If &RowBasis == &ColBasis it is assumed that Rdm is symmetric. If it is not symmetric, pls
188 // symmetrize first (not checked!).
189 void AccGradient2ix( double *pGrad, FMatrixView const &Rdm, FRawBasis const &RowBasis,
190    FRawBasis const &ColBasis, FKrn2i const &Krn2i, FMemoryStack &Mem, double Prefactor=1.);
191 
192 // evaluate <a|krn Laplace^LaplaceOrder|b> by directly calling
193 // a IR integral kernel.
194 struct FKrn2i_Direct : public FKrn2i
195 {
196    FKrn2i_Direct(ir::FIntegralKernel const *pIrKernel_, uint LaplaceOrder=0, double Prefactor=1.);
197    void EvalInt2e2c(double *pOut, size_t StrideA, size_t StrideC, ir::FRawShell const *pA, ir::FRawShell const *pC, double Prefactor, bool Add, ct::FMemoryStack &Mem) const;
198 //    void EvalInt2e2c1d(double *pOutDa, double *pOutDb, size_t StrideA, size_t StrideC, size_t StrideDeriv, ir::FRawShell const *pA, ir::FRawShell const *pC, double Prefactor, bool Add, ct::FMemoryStack &Mem) const;
199    void AccCoreDeriv1d(double *pGrad, double const *pRdmAC, size_t StrideA, size_t StrideC, ir::FRawShell const *pA, size_t iCenterA, ir::FRawShell const *pC, size_t iCenterC, double Prefactor, ct::FMemoryStack &Mem) const;
200 protected:
201    ir::FIntegralKernel const
202       *m_pIrKernel;
203    uint
204       m_LaplaceOrder;
205    double
206       m_Prefactor;
207 };
208 
209 // evaluate <a| v(r) |b> where v(r) is sum of the fields emitted via pIrKernel
210 // for the sum of the supplied spherical point multipoles.
211 // Can be used, for example, to make nuclear attraction integrals.
212 struct FKrn2i_PointMultipoles : public FKrn2i
213 {
214    // each point should provide 2*l+1 coefficients; one for each of its spherical
215    // components. Slm order is equal to the order used for basis functions.
216    // (in particular, dipole components are x,y,z, not x,z,y.).
217    struct FPoint {
218       FVector3 vPos;
219       uint l;
220       ptrdiff_t iCenter;
221       // ^- for gradients: identify the atomic center index on which this multipole sits. Set to -1 if None.
222       //    note: even if having external lattices, one could just append their gradient array to the QM system gradient array
223       //    and then do both in one go.
224    };
225    FKrn2i_PointMultipoles(ir::FIntegralKernel const *pIrKernel_, FPoint const *pPoints, double const *pCoeffs, size_t nPoints);
226    void EvalInt2e2c(double *pOut, size_t StrideA, size_t StrideC, ir::FRawShell const *pA, ir::FRawShell const *pC, double Prefactor, bool Add, ct::FMemoryStack &Mem) const;
227 protected:
228    ir::FIntegralKernel const
229       *m_pIrKernel;
230    TArray<ir::FRawShell>
231       m_PointShells;
232    TArray<double>
233       m_Data;
234    TArray<ptrdiff_t>
235       m_PointCenters;
236    double
237       m_Exp;
238    uint
239       m_MaxL;
240 };
241 
242 
243 
244 } // namespace ct
245 
246 #endif // CT8K_BASISSET_H
247