1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
8 //                    Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
9 //
10 // File created by: Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
11 //////////////////////////////////////////////////////////////////////////////////////
12 
13 
14 #ifndef QMCPLUSPLUS_ORBITAL_IMAGES_H
15 #define QMCPLUSPLUS_ORBITAL_IMAGES_H
16 
17 #include "QMCHamiltonians/OperatorBase.h"
18 #include "QMCWaveFunctions/SPOSet.h"
19 #include "QMCWaveFunctions/WaveFunctionFactory.h"
20 
21 namespace qmcplusplus
22 {
23 /** "Estimator" to produce files for orbital plotting.
24  *  Orbitals are evaluated on a uniform grid and written to ascii files.
25  *  Only format currently supported is xsf (XCrySDen).
26  *  Can print real, imag, abs, and abs^2 of each orbital.
27  *  This class should work with any SPOSet.
28  *  All work is performed by omp thread 0 of mpi task 0.
29  *
30  *  TO USE THIS YOU NEED TO KNOW THE NAME OF THE SPOSET(S)!!!
31  *    For example, using sposet_builder the names are prominently displayed:
32  *
33  *      <sposet_builder type="bspline" href="pwscf.h5" tilematrix="1 0 0 0 1 0 0 0 1" twistnum="0" source="ion0" version="0.10">
34  *         <sposet type="bspline" name="spo_ud" size="2" spindataset="0"/>
35  *      </sposet_builder>
36  *
37  *    In this case a single sposet named "spo_ud" exists.
38  *
39  *    If you are using Slater-Jastrow w/o sposet_builder
40  *    the sposets should be named updet and downdet.
41  *
42  *
43  *  To make xsf files, add xml similar to the following to <hamiltonian/>:
44  *
45  *    minimal working example: single sposet named spo_ud on a 20x20x20 grid
46  *
47  *      <estimator name="OrbitalImages" type="orbitalimages" ions="ion0">
48  *        <parameter name="sposets"> spo_ud </parameter>
49  *        <parameter name="grid"> 20 20 20 </parameter>
50  *      </estimator>
51  *
52  *    as above, but print real, imag, abs, and abs^2 of orbitals (default real)
53  *
54  *      <estimator name="OrbitalImages" type="orbitalimages" ions="ion0">
55  *        <parameter name="sposets"> spo_ud </parameter>
56  *        <parameter name="grid"> 20 20 20 </parameter>
57  *        <parameter name="value"> real imag abs abs2 </parameter>
58  *      </estimator>
59  *
60  *    up and down sposets named spo_u and spo_d w/ individual orbitals selected
61  *
62  *      <estimator name="OrbitalImages" type="orbitalimages" ions="ion0">
63  *        <parameter name="sposets"> spo_u spo_d </parameter>
64  *        <parameter name="spo_u">   13 24 37    </parameter>
65  *        <parameter name="spo_d">   10 18 29 41 </parameter>
66  *        <parameter name="grid">    20 20 20    </parameter>
67  *        <parameter name="value"> real imag abs abs2 </parameter>
68  *      </estimator>
69  *
70  *    user defined cell by cell center and axes (openbc's, subset of pbc cell)
71  *
72  *      <estimator name="OrbitalImages" type="orbitalimages" ions="ion0">
73  *        <parameter name="sposets"> spo_ud </parameter>
74  *        <parameter name="grid"> 20 20 20 </parameter>
75  *        <parameter name="center"> 2.5 2.5 2.5 </parameter>
76  *        <parameter name="cell">
77  *           5.0  0.0  0.0
78  *           0.0  5.0  0.0
79  *           0.0  0.0  5.0
80  *        </parameter>
81  *      </estimator>
82  *
83  *    user defined cell by cell corner and axes (default corner is 0 0 0)
84  *
85  *      <estimator name="OrbitalImages" type="orbitalimages" ions="ion0">
86  *        <parameter name="sposets"> spo_ud </parameter>
87  *        <parameter name="grid"> 20 20 20 </parameter>
88  *        <parameter name="corner"> 0.0 0.0 0.0 </parameter>
89  *        <parameter name="cell">
90  *           5.0  0.0  0.0
91  *           0.0  5.0  0.0
92  *           0.0  0.0  5.0
93  *        </parameter>
94  *      </estimator>
95  *
96  *    store only a batch of orbital values in memory
97  *    this can save on memory for very large grids
98  *    but will evaluate all orbitals #orbitals/batch_size times
99  *
100  *      <estimator name="OrbitalImages" type="orbitalimages" ions="ion0">
101  *        <parameter name="sposets"> spo_ud </parameter>
102  *        <parameter name="grid"> 200 200 200 </parameter>
103  *        <parameter name="batch_size"> 10    </parameter>
104  *      </estimator>
105  *
106  */
107 class OrbitalImages : public OperatorBase
108 {
109 public:
110   enum
111   {
112     DIM = OHMMS_DIM
113   };
114 
115   typedef SPOSet::ValueVector_t ValueVector_t;
116   typedef SPOSet::GradVector_t GradVector_t;
117   typedef ParticleSet::ParticleLayout_t Lattice_t;
118   typedef std::map<std::string, ParticleSet*> PSPool;
119 
120 
121   ///derivative types
122   enum derivative_types_enum
123   {
124     value_d = 0,
125     gradient_d,
126     laplacian_d
127   };
128 
129 
130   ///options for orbital value to write
131   enum value_types_enum
132   {
133     real_val = 0,
134     imag_val,
135     abs_val,
136     abs2_val
137   };
138 
139   ///options for orbital output file format
140   enum formats_enum
141   {
142     xsf = 0
143   };
144 
145   ///at put() ion particleset is obtained from ParticleSetPool
146   PSPool& psetpool;
147 
148   ///electron particleset
149   ParticleSet* Peln;
150 
151   ///ion particleset
152   ParticleSet* Pion;
153 
154   ///mpi communicator
155   Communicate* comm;
156 
157 
158   ///file format selection
159   formats_enum format;
160 
161   ///orbital value selections
162   std::vector<value_types_enum> value_types;
163 
164   ///write out derivatives in addition to values
165   bool derivatives;
166 
167   ///names of sposets to evaluate
168   std::vector<std::string> sposet_names;
169 
170   ///indices of orbitals within each sposet to evaluate
171   std::vector<std::vector<int>*> sposet_indices;
172 
173   ///sposets obtained by name from WaveFunctionFactory
174   std::vector<SPOSet*> sposets;
175 
176   ///evaluate points at grid cell centers instead of edges
177   bool center_grid;
178 
179   ///cell bounding the evaluation grid, default is simulation cell
180   Lattice_t cell;
181 
182   ///location of cell corner, positions in the cell are corner+uvec*cell
183   PosType corner;
184 
185   ///number of grid points in each direction (cell axis)
186   TinyVector<int, DIM> grid;
187 
188   ///stride to generate grid in arbitrary dimensions
189   TinyVector<int, DIM> gdims;
190 
191   ///total number of grid points
192   int npoints;
193 
194   ///number of orbitals to store in memory at a time (batch_size*npoints)
195   int batch_size;
196 
197   ///temporary vector to hold values of all orbitals at a single point
198   ValueVector_t spo_vtmp;
199 
200   ///temporary vector to hold gradients of all orbitals at a single point
201   GradVector_t spo_gtmp;
202 
203   ///temporary vector to hold laplacians of all orbitals at a single point
204   ValueVector_t spo_ltmp;
205 
206   ///temporary array to hold values of a batch of orbitals at all grid points
207   Array<ValueType, 2> batch_values;
208 
209   ///temporary array to hold gradients of a batch of orbitals at all grid points
210   Array<GradType, 2> batch_gradients;
211 
212   ///temporary array to hold laplacians of a batch of orbitals at all grid points
213   Array<ValueType, 2> batch_laplacians;
214 
215   ///temporary array to hold values of a single orbital at all grid points
216   std::vector<ValueType> orbital;
217 
218 
219   //constructor/destructor
220   OrbitalImages(ParticleSet& P, PSPool& PSP, Communicate* mpicomm, const WaveFunctionFactory& factory);
~OrbitalImages()221   ~OrbitalImages(){};
222 
223   //standard interface
224   OperatorBase* makeClone(ParticleSet& P, TrialWaveFunction& psi);
225 
226   ///read xml input
227   bool put(xmlNodePtr cur);
228 
229   ///hijack estimator evaluate to evaluate and write all orbitals
230   Return_t evaluate(ParticleSet& P);
231 
232   //optional standard interface
233   //void get_required_traces(TraceManager& tm);
234   //void setRandomGenerator(RandomGenerator_t* rng);
235 
236   //required for Collectables interface
addObservables(PropertySetType & plist,BufferType & olist)237   void addObservables(PropertySetType& plist, BufferType& olist) {}
registerCollectables(std::vector<observable_helper * > & h5desc,hid_t gid)238   void registerCollectables(std::vector<observable_helper*>& h5desc, hid_t gid) const {}
239 
240   //should be empty for Collectables interface
resetTargetParticleSet(ParticleSet & P)241   void resetTargetParticleSet(ParticleSet& P) {}
setObservables(PropertySetType & plist)242   void setObservables(PropertySetType& plist) {}
setParticlePropertyList(PropertySetType & plist,int offset)243   void setParticlePropertyList(PropertySetType& plist, int offset) {}
244 #if !defined(REMOVE_TRACEMANAGER)
checkout_scalar_arrays(TraceManager & tm)245   void checkout_scalar_arrays(TraceManager& tm) {}
collect_scalar_samples()246   void collect_scalar_samples() {}
delete_scalar_arrays()247   void delete_scalar_arrays() {}
248 #endif
249 
250   //obsolete?
get(std::ostream & os)251   bool get(std::ostream& os) const { return false; }
252 
253   //local functions
254   ///write brief report of configuration data
255   void report(const std::string& pad);
256 
257   ///write a single orbital to file
258   void write_orbital(const std::string& sponame,
259                      int index,
260                      std::vector<ValueType>& orb,
261                      value_types_enum value_type,
262                      derivative_types_enum derivative_type = value_d,
263                      int dimension                         = 0);
264 
265   ///write a single orbital to an xsf file
266   void write_orbital_xsf(const std::string& sponame,
267                          int index,
268                          std::vector<ValueType>& orb,
269                          value_types_enum value_type,
270                          derivative_types_enum derivative_type = value_d,
271                          int dimension                         = 0);
272 
273 private:
274   /// reference to the sposet_builder_factory
275   const WaveFunctionFactory& wf_factory_;
276 };
277 
278 } // namespace qmcplusplus
279 
280 #endif
281