IM-SRG++  0
ModelSpace.hh
1 #ifndef ModelSpace_h
2 #define ModelSpace_h 1
3 
4 //#define ARMA_NO_DEBUG
5 #include <vector>
6 #include <unordered_map>
7 #include <map>
8 #include <armadillo>
9 #include "IMSRGProfiler.hh"
10 #ifndef SQRT2
11  #define SQRT2 1.4142135623730950488
12 #endif
13 #define OCC_CUT 1e-6
14 
15 
16 using namespace std;
17 
18 typedef unsigned long long int index_t;
19 
20 class ModelSpace; //forward declaration so Ket can use ModelSpace
21 
22 //struct Orbit
23 class Orbit
24 {
25  public:
26  // Fields
27  int n;
28  int l;
29  int j2;
30  int tz2;
31  double occ; // particle=0, hole=1
32 // int ph; // particle=0, hole=1
33 // int io; // inside=0, outside=1
34  int cvq; // core=0, valence=1, qspace=2
35  int index;
36 
37  //Constructors
38  ~Orbit();
39  Orbit();
40  Orbit(int n ,int l, int j, int t, double occ, int cvq, int index);
41 // Orbit(int n ,int l, int j, int t, int ph, int cvq, int index);
42  Orbit(const Orbit&);
43 // void swap(Orbit&) throw();
44 // Orbit& operator=( const Orbit& );
45  // Methods
46 };
47 
48 
49 
50 
51 class Ket // | pq >
52 {
53 
54  public:
55  // Constructors
56  ~Ket();
57  Ket();
58  Ket(Orbit& op, Orbit& oq);
59  // Methods
60  int Phase(int J);
61  int delta_pq(){return dpq;};
62 
63  // Fields
64  Orbit* op;
65  Orbit* oq;
66  int p;
67  int q;
68  private:
69  // Fields
70  int dpq;
71  int phase_prefactor;
72 
73 };
74 
75 
76 
77 
79 {
80  public:
81  //Fields
82  int J;
83  int parity;
84  int Tz;
85 
86  // Constructors
87  ~TwoBodyChannel();
89  TwoBodyChannel(int j, int p, int t, ModelSpace* ms);
90  TwoBodyChannel(int N, ModelSpace* ms);
91  void Initialize(int N, ModelSpace* ms);
92 
93 
94  //Methods
95  int GetNumberKets() const {return NumberKets;};
96  int GetLocalIndex(int ketindex) const { return KetMap[ketindex];}; // modelspace ket index => local ket index
97  int GetLocalIndex(int p, int q) const ;
98  int GetKetIndex(int i) const { return KetList[i];}; // local ket index => modelspace ket index
99  Ket& GetKet(int i) ; // get pointer to ket using local index
100 
101  arma::uvec KetIndex_pp ;
102  arma::uvec KetIndex_hh ;
103  arma::uvec KetIndex_ph ;
104  arma::uvec KetIndex_cc ;
105  arma::uvec KetIndex_vc ;
106  arma::uvec KetIndex_qc ;
107  arma::uvec KetIndex_vv ;
108  arma::uvec KetIndex_qv ;
109  arma::uvec KetIndex_qq ;
110  arma::vec Ket_occ_hh;
111  arma::vec Ket_unocc_hh;
112  arma::vec Ket_occ_ph;
113  arma::vec Ket_unocc_ph;
114 
115 
116 
117  arma::uvec GetKetIndexFromList(vector<index_t>& vec_in);
118  arma::uvec& GetKetIndex_pp();
119  arma::uvec& GetKetIndex_hh();
120  arma::uvec& GetKetIndex_ph();
121  arma::uvec& GetKetIndex_cc();
122  arma::uvec& GetKetIndex_vc();
123  arma::uvec& GetKetIndex_qc();
124  arma::uvec& GetKetIndex_vv();
125  arma::uvec& GetKetIndex_qv();
126  arma::uvec& GetKetIndex_qq();
127 
128 
129 // private:
130  //Fields
131  ModelSpace * modelspace;
132  int NumberKets; // Number of pq configs that participate in this channel
133  vector<int> KetList; // eg [2, 4, 7, ...] Used for looping over all the kets in the channel
134  vector<int> KetMap; // eg [ -1, -1, 0, -1, 1, -1, -1, 2 ...] Used for asking what is the local index of this ket. -1 means the ket doesn't participate in this channel
135  //Methods
136  virtual bool CheckChannel_ket(Orbit* op, Orbit* oq) const; // check if |pq> participates in this channel
137  bool CheckChannel_ket(Ket &ket) const {return CheckChannel_ket(ket.op,ket.oq);}; // check if |pq> participates in this channel
138 
139 };
140 
141 
142 
144 {
145  public:
148  TwoBodyChannel_CC(int j, int p, int t, ModelSpace* ms);
149  TwoBodyChannel_CC(int N, ModelSpace* ms);
150  bool CheckChannel_ket(Orbit* op, Orbit* oq) const; // check if |pq> participates in this channel
151 // bool CheckChannel_ket(int p, int q) const; // check if |pq> participates in this channel
152 };
153 
154 
155 
156 
158 {
159 
160  public:
161 
162 
163  // Constructors
164  ~ModelSpace();
165  ModelSpace();
166  ModelSpace(const ModelSpace&); // copy constructor
167  ModelSpace( ModelSpace&&); // move constructor
168  ModelSpace(int emax, vector<string> hole_list, vector<string> valence_list);
169  ModelSpace(int emax, vector<string> hole_list, vector<string> core_list, vector<string> valence_list);
170  ModelSpace(int emax, string reference, string valence);
171  ModelSpace(int emax, string reference);
172 
173  // Overloaded operators
174  ModelSpace operator=(const ModelSpace&);
175  ModelSpace operator=(ModelSpace&&);
176 
177  // Methods
178 
179  void Init(int emax, string reference, string valence);
180  void Init(int emax, map<index_t,double> hole_list, string valence);
181  void Init(int emax, map<index_t,double> hole_list, vector<index_t> core_list, vector<index_t> valence_list);
182  void Init(int emax, vector<string> hole_list, vector<string> core_list, vector<string> valence_list);
183  void Init_occ_from_file(int emax, string valence, string occ_file);
184 
185 // vector<index_t> GetOrbitsAZ(int A, int Z);
186  map<index_t,double> GetOrbitsAZ(int A, int Z);
187  void GetAZfromString(string str, int& A, int& Z);
188  vector<index_t> String2Index( vector<string> vs );
189  string Index2String(index_t ind);
190  void Get0hwSpace(int Aref, int Zref, vector<index_t>& core_list, vector<index_t>& valence_list);
191  void ParseCommaSeparatedValenceSpace(string valence, vector<index_t>& core_list, vector<index_t>& valence_list);
192 
193  void SetupKets();
194  void AddOrbit(Orbit orb);
195  void AddOrbit(int n, int l, int j2, int tz2, double occ, int io);
196  // Setter/Getters
197  Orbit& GetOrbit(int i) {return (Orbit&) Orbits[i];};
198 // Orbit& GetOrbit(int i) const {return (Orbit&) Orbits[i];};
199  Ket& GetKet(int i) const {return (Ket&) Kets[i];};
200  Ket& GetKet(int p, int q) const {return (Ket&) Kets[Index2(p,q)];};
201  int GetOrbitIndex(int n, int l, int j2, int tz2) const {return Index1(n,l,j2,tz2);};
202  int GetKetIndex(int p, int q) const {return Index2(p,q);}; // convention is p<=q
203  int GetKetIndex(Ket * ket) const {return Index2(ket->p,ket->q);}; // convention is p<=q
204  int GetNumberOrbits() const {return norbits;};
205  int GetNumberKets() const {return Kets.size();};
206  void SetHbarOmega(double hw) {hbar_omega = hw;};
207  void SetTargetMass(int A) {target_mass = A;};
208  void SetTargetZ(int Z) {target_Z = Z;};
209  double GetHbarOmega() const {return hbar_omega;};
210  int GetTargetMass() const {return target_mass;};
211  int GetTargetZ() const {return target_Z;};
212  int GetAref() const {return Aref;};
213  int GetZref() const {return Zref;};
214  int GetNumberTwoBodyChannels() const {return TwoBodyChannels.size();};
215  TwoBodyChannel& GetTwoBodyChannel(int ch) const {return (TwoBodyChannel&) TwoBodyChannels[ch];};
216  TwoBodyChannel_CC& GetTwoBodyChannel_CC(int ch) const {return (TwoBodyChannel_CC&) TwoBodyChannels_CC[ch];};
217  inline int GetTwoBodyJmax() const {return TwoBodyJmax;};
218  inline int GetThreeBodyJmax() const {return ThreeBodyJmax;};
219  void SetReference(vector<index_t>);
220  void SetReference(map<index_t,double>);
221  void SetReference(string);
222 
223  int GetEmax(){return Emax;};
224  int GetE2max(){return E2max;};
225  int GetE3max(){return E3max;};
226  int GetLmax2(){return Lmax2;};
227  int GetLmax3(){return Lmax3;};
228  void SetEmax(int e){Emax=e;};
229  void SetE2max(int e){E2max=e;};
230  void SetE3max(int e){E3max=e;};
231  void SetLmax2(int l){Lmax2=l;};
232  void SetLmax3(int l){Lmax3=l;};
233 
234  double GetSixJ(double j1, double j2, double j3, double J1, double J2, double J3);
235  double GetNineJ(double j1, double j2, double j3, double j4, double j5, double j6, double j7, double j8, double j9);
236  double GetMoshinsky( int N, int Lam, int n, int lam, int n1, int l1, int n2, int l2, int L); // Inconsistent notation. Not ideal.
237  bool SixJ_is_empty(){ return SixJList.empty(); };
238 
239  int GetOrbitIndex(string);
240  int GetTwoBodyChannelIndex(int j, int p, int t);
241  inline int phase(int x) {return (x%2)==0 ? 1 : -1;};
242  inline int phase(double x) {return phase(int(x));};
243 
244  inline int Index1(int n, int l, int j2, int tz2) const {return(2*n+l)*(2*n+l+3) + 1-j2 + (tz2+1)/2 ;};
245 // inline int Index2(int p, int q) const {return q*(q+1)/2 + p;};
246  inline int Index2(int p, int q) const {return p*(2*norbits-1-p)/2 + q;};
247 
248  void PreCalculateMoshinsky();
249  void ClearVectors();
250  void CalculatePandyaLookup(int rank_J, int rank_T, int parity); // construct a lookup table for more efficient pandya transformation
251 // map<array<int,2>,vector<array<int,2>>>& GetPandyaLookup(int rank_J, int rank_T, int parity);
252  map<array<int,2>,array<vector<int>,2>>& GetPandyaLookup(int rank_J, int rank_T, int parity);
253 
254 
255  // Data members
256  vector<index_t> holes; // in the reference Slater determinant
257 // map<index_t,double> holes; // in the reference Slater determinant
258  vector<index_t> particles; // above the reference Slater determinant
259  vector<index_t> core; // core for decoupling
260  vector<index_t> valence; // valence space for decoupling
261  vector<index_t> qspace; // above the valence space for decoupling
262  vector<index_t> proton_orbits;
263  vector<index_t> neutron_orbits;
264 
265  vector<index_t> KetIndex_pp;
266  vector<index_t> KetIndex_ph;
267  vector<index_t> KetIndex_hh;
268  vector<index_t> KetIndex_cc;
269  vector<index_t> KetIndex_vc;
270  vector<index_t> KetIndex_qc;
271  vector<index_t> KetIndex_vv;
272  vector<index_t> KetIndex_qv;
273  vector<index_t> KetIndex_qq;
274  vector<double> Ket_occ_hh;
275  vector<double> Ket_unocc_hh;
276  vector<double> Ket_occ_ph;
277  vector<double> Ket_unocc_ph;
278 
279 // array< array< vector<index_t>, 2>,3> MonopoleKets; //List of kets of a given Tz,parity
280  array< array< unordered_map<index_t,index_t>, 2>,3> MonopoleKets; //List of kets of a given Tz,parity
281 
282  int Emax;
283  int E2max;
284  int E3max;
285  int Lmax2;
286  int Lmax3;
287  int OneBodyJmax;
288  int TwoBodyJmax;
289  int ThreeBodyJmax;
290  map<array<int,3>,vector<index_t> > OneBodyChannels;
291 
292  vector<unsigned int> SortedTwoBodyChannels;
293  vector<unsigned int> SortedTwoBodyChannels_CC;
294 
295  static map<string,vector<string>> ValenceSpaces;
296 // map< array<int,3>, map< array<int,2>,vector<array<int,2>> > > PandyaLookup;
297 
298 
299 // private:
300  // Fields
301  int norbits;
302  double hbar_omega;
303  int target_mass;
304  int target_Z;
305  int Aref;
306  int Zref;
307  int nTwoBodyChannels;
308  vector<Orbit> Orbits;
309  vector<Ket> Kets;
310  vector<TwoBodyChannel> TwoBodyChannels;
311  vector<TwoBodyChannel_CC> TwoBodyChannels_CC;
312  map< array<int,3>, map< array<int,2>,array<vector<int>,2> > > PandyaLookup;
313  bool moshinsky_has_been_precalculated;
314  IMSRGProfiler profiler;
315 // map<long int,double> SixJList;
316 
317  static unordered_map<unsigned long int,double> SixJList;
318  static unordered_map<long long unsigned int,double> NineJList;
319  static unordered_map<long long unsigned int,double> MoshList;
320 
321 };
322 
323 
324 
325 
326 #endif
327 
Definition: ModelSpace.hh:23
Definition: ReadWrite.cc:3012
Definition: IMSRGProfiler.hh:17
Definition: ModelSpace.hh:78
Definition: ModelSpace.hh:143
Definition: ModelSpace.hh:157
Definition: ModelSpace.hh:51