1 /*
2    CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3    Copyright (C) 2013-2018 Sebastian Wouters
4 
5    This program 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; either version 2 of the License, or
8    (at your option) any later version.
9 
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU General Public License for more details.
14 
15    You should have received a copy of the GNU General Public License along
16    with this program; if not, write to the Free Software Foundation, Inc.,
17    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19 
20 #ifndef TENSOROPERATOR_CHEMPS2_H
21 #define TENSOROPERATOR_CHEMPS2_H
22 
23 #include "Tensor.h"
24 #include "TensorT.h"
25 #include "SyBookkeeper.h"
26 
27 namespace CheMPS2{
28 /** TensorOperator class.
29     \author Sebastian Wouters <sebastianwouters@gmail.com>
30     \date October 15, 2015
31 
32     The TensorOperator class is a storage and update class for tensor operators with a given:\n
33     - spin (two_j)
34     - particle number (n_elec)
35     - point group irrep (n_irrep).
36 
37     It replaces the previous classes TensorDiag, TensorSwap, TensorS0Abase, TensorS1Bbase, TensorF0Cbase, TensorF1Dbase, TensorA, TensorB, TensorC, and TensorD. Their storage and update functions have a common origin. The boolean prime_last denotes whether in which convention the tensor operator is stored:\n
38     - \f$ \braket{ j m J M | j' m' } \braket{ j ( N, I ) || J ( n\_elec, n\_irrep ) || j' ( N + n\_elec, I \times n\_irrep ) } \f$ (prime_last == true)
39     - \f$ \braket{ j' m' J M | j m } \braket{ j ( N, I ) || J ( n\_elec, n\_irrep ) || j' ( N + n\_elec, I \times n\_irrep ) } \f$ (prime_last == false).
40 
41     This determines the specific reduced update formulae when contracting with the Clebsch-Gordan coefficients of the reduced MPS tensors. */
42    class TensorOperator : public Tensor{
43 
44       public:
45 
46          //! Constructor
47          /** \param boundary_index The boundary index
48              \param two_j Twice the spin of the tensor operator
49              \param n_elec How many electrons there are more in the symmetry sector of the lower leg compared to the upper leg
50              \param n_irrep The (real-valued abelian) point group irrep difference between the symmetry sectors of the lower and upper legs (see Irreps.h)
51              \param moving_right If true: sweep from left to right. If false: sweep from right to left.
52              \param prime_last Convention in which the tensor operator is stored (see class information)
53              \param jw_phase Whether or not to include a Jordan-Wigner phase due to the fermion anti-commutation relations
54              \param bk_up   Symmetry bookkeeper of the upper MPS
55              \param bk_down Symmetry bookkeeper of the lower MPS */
56          TensorOperator( const int boundary_index, const int two_j, const int n_elec, const int n_irrep, const bool moving_right, const bool prime_last, const bool jw_phase, const SyBookkeeper * bk_up, const SyBookkeeper * bk_down );
57 
58          //! Destructor
59          virtual ~TensorOperator();
60 
61          //! Get the number of symmetry blocks
62          /** \return The number of symmetry blocks */
63          int gNKappa() const;
64 
65          //! Get the pointer to the storage
66          /** return pointer to the storage */
67          double * gStorage();
68 
69          //! Get the index corresponding to a certain tensor block
70          /** \param N1 The up particle number sector
71              \param TwoS1 The up spin symmetry sector
72              \param I1 The up irrep sector
73              \param N2 The down particle number sector
74              \param TwoS2 The down spin symmetry sector
75              \param I2 The down irrep sector
76              \return The kappa corresponding to the input parameters; -1 means no such block */
77          int gKappa( const int N1, const int TwoS1, const int I1, const int N2, const int TwoS2, const int I2 ) const;
78 
79          //! Get the storage jump corresponding to a certain tensor block
80          /** \param kappa The symmetry block
81              \return kappa2index[kappa], the memory jumper to a certain block */
82          int gKappa2index( const int kappa ) const;
83 
84          //! Get the pointer to the storage of a certain tensor block
85          /** \param N1 The up particle number sector
86              \param TwoS1 The up spin symmetry sector
87              \param I1 The up irrep sector
88              \param N2 The down particle number sector
89              \param TwoS2 The down spin symmetry sector
90              \param I2 The down irrep sector
91              \return Pointer to the storage of the specified tensor block; NULL means no such block */
92          double * gStorage( const int N1, const int TwoS1, const int I1, const int N2, const int TwoS2, const int I2 );
93 
94          //! Get the boundary index
95          /** \return the index */
96          int gIndex() const;
97 
98          //! Get twice the spin of the tensor operator
99          /** \return Twice the spin of the tensor operatorg */
100          int get_2j() const;
101 
102          //! Get how many electrons there are more in the symmetry sector of the lower leg compared to the upper leg
103          /** \return How many electrons there are more in the symmetry sector of the lower leg compared to the upper leg */
104          int get_nelec() const;
105 
106          //! Get the (real-valued abelian) point group irrep difference between the symmetry sectors of the lower and upper legs (see Irreps.h)
107          /** \return The (real-valued abelian) point group irrep difference between the symmetry sectors of the lower and upper legs (see Irreps.h) */
108          int get_irrep() const;
109 
110          //! Clear and update
111          /** \param previous The previous TensorOperator needed for the update
112              \param mps_tensor_up   The upper MPS tensor needed for the update
113              \param mps_tensor_down The lower MPS tensor needed for the update
114              \param workmem Work memory */
115          void update( TensorOperator * previous, TensorT * mps_tensor_up, TensorT * mps_tensor_down, double * workmem );
116 
117          //! daxpy for TensorOperator
118          /** \param alpha The prefactor
119              \param to_add The TensorOperator x which should be added: this <-- this + alpha * to_add */
120          void daxpy( double alpha, TensorOperator * to_add );
121 
122          //! daxpy_transpose for C- and D-tensors (with special spin-dependent factors)
123          /** \param alpha The prefactor
124              \param to_add The TensorOperator x which should be added: this <-- this + alpha * special_spin_dependent_factor * to_add^T */
125          void daxpy_transpose_tensorCD( const double alpha, TensorOperator * to_add );
126 
127          //! Set all storage variables to 0.0
128          void clear();
129 
130          //! Make the in-product of two TensorOperator
131          /** \param buddy The second tensor
132              \param trans If trans == 'N' a regular ddot is taken. If trans == 'T' and n_elec==0, the in-product with buddy's transpose is made.
133              \return The in-product */
134          double inproduct( TensorOperator * buddy, const char trans ) const;
135 
136       protected:
137 
138          //! The bookkeeper of the upper MPS
139          const SyBookkeeper * bk_up;
140 
141          //! The bookkeeper of the lower MPS
142          const SyBookkeeper * bk_down;
143 
144          //! Twice the spin of the tensor operator
145          int two_j;
146 
147          //! How many electrons there are more in the symmetry sector of the lower leg compared to the upper leg
148          int n_elec;
149 
150          //! The (real-valued abelian) point group irrep difference between the symmetry sectors of the lower and upper legs (see Irreps.h)
151          int n_irrep;
152 
153          //! Whether or not moving right
154          bool moving_right;
155 
156          //! The up particle number sector
157          int * sector_nelec_up;
158 
159          //! The up spin symmetry sector
160          int * sector_irrep_up;
161 
162          //! The up spin symmetry sector
163          int * sector_spin_up;
164 
165          //! The down spin symmetry sector (pointer points to sectorTwoS1 if two_j == 0)
166          int * sector_spin_down;
167 
168          //! Update moving right
169          /** \param ikappa The tensor block which should be updated
170              \param previous The previous TensorOperator needed for the update
171              \param mps_tensor_up   The upper MPS tensor needed for the update
172              \param mps_tensor_down The lower MPS tensor needed for the update
173              \param workmem Work memory */
174          void update_moving_right( const int ikappa, TensorOperator * previous, TensorT * mps_tensor_up, TensorT * mps_tensor_down, double * workmem );
175 
176          //! Update moving left
177          /** \param ikappa The tensor block which should be updated
178              \param previous The previous TensorOperator needed for the update
179              \param mps_tensor_up   The upper MPS tensor needed for the update
180              \param mps_tensor_down The lower MPS tensor needed for the update
181              \param workmem Work memory */
182          void update_moving_left( const int ikappa, TensorOperator * previous, TensorT * mps_tensor_up, TensorT * mps_tensor_down, double * workmem );
183 
184          //! Convention in which the tensor operator is stored (see class information)
185          bool prime_last;
186 
187          //! Whether or not to include a Jordan-Wigner phase due to the fermion anti-commutation relations
188          bool jw_phase;
189 
190       private:
191 
192    };
193 }
194 
195 #endif
196