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: Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
8 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
9 //
10 // File created by: Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
11 //////////////////////////////////////////////////////////////////////////////////////
12
13
14 #include "ReferencePoints.h"
15 #include "Utilities/string_utils.h"
16 #include "OhmmsData/AttributeSet.h"
17
18 namespace qmcplusplus
19 {
put(xmlNodePtr cur,ParticleSet & P,std::vector<ParticleSet * > & Pref)20 bool ReferencePoints::put(xmlNodePtr cur, ParticleSet& P, std::vector<ParticleSet*>& Pref)
21 {
22 app_log() << " Entering ReferencePoints::put" << std::endl;
23 bool succeeded = true;
24 put(P, Pref);
25 OhmmsAttributeSet ra;
26 std::string coord = "";
27 ra.add(coord, "coord");
28 ra.put(cur);
29 for (int i = 0; i < DIM; i++)
30 for (int d = 0; d < DIM; d++)
31 axes(d, i) = P.Lattice.a(i)[d];
32 Tensor_t crd;
33 if (coord == "cell")
34 {
35 coordinate = cellC;
36 crd = axes;
37 }
38 else if (coord == "cartesian")
39 {
40 coordinate = cartesianC;
41 for (int i = 0; i < DIM; i++)
42 for (int d = 0; d < DIM; d++)
43 if (d == i)
44 crd(i, i) = 1.0;
45 else
46 crd(d, i) = 0.0;
47 }
48 else
49 {
50 app_log() << std::endl;
51 app_log() << " Valid coordinates must be provided for element reference_points." << std::endl;
52 app_log() << " You provided: " << coord << std::endl;
53 app_log() << " Options are cell or cartesian." << std::endl;
54 app_log() << std::endl;
55 succeeded = false;
56 }
57 //read in the point contents
58 app_log() << " reading reference_points contents" << std::endl;
59 std::vector<std::string> lines = split(strip(XMLNodeString{cur}), "\n");
60 for (int i = 0; i < lines.size(); i++)
61 {
62 std::vector<std::string> tokens = split(strip(lines[i]), " ");
63 if (tokens.size() != DIM + 1)
64 {
65 app_log() << " reference point has 4 entries, given " << tokens.size() << ": " << lines[i] << std::endl;
66 succeeded = false;
67 }
68 else
69 {
70 Point rp;
71 for (int d = 0; d < DIM; d++)
72 {
73 rp[d] = string2real(tokens[d + 1]);
74 }
75 rp = dot(crd, rp);
76 points[tokens[0]] = rp;
77 }
78 }
79 return succeeded;
80 }
81
put(ParticleSet & P,std::vector<ParticleSet * > & Psets)82 bool ReferencePoints::put(ParticleSet& P, std::vector<ParticleSet*>& Psets)
83 {
84 //get axes and origin information from the ParticleSet
85 points["zero"] = 0 * P.Lattice.a(0);
86 points["a1"] = P.Lattice.a(0);
87 points["a2"] = P.Lattice.a(1);
88 points["a3"] = P.Lattice.a(2);
89 //points["center"]= .5*(P.Lattice.a(0)+P.Lattice.a(1)+P.Lattice.a(2))
90 //set points on face centers
91 points["f1p"] = points["zero"] + .5 * points["a1"];
92 points["f1m"] = points["zero"] - .5 * points["a1"];
93 points["f2p"] = points["zero"] + .5 * points["a2"];
94 points["f2m"] = points["zero"] - .5 * points["a2"];
95 points["f3p"] = points["zero"] + .5 * points["a3"];
96 points["f3m"] = points["zero"] - .5 * points["a3"];
97 //set points on cell corners
98 points["cmmm"] = points["zero"] + .5 * (-1 * points["a1"] - points["a2"] - points["a3"]);
99 points["cpmm"] = points["zero"] + .5 * (points["a1"] - points["a2"] - points["a3"]);
100 points["cmpm"] = points["zero"] + .5 * (-1 * points["a1"] + points["a2"] - points["a3"]);
101 points["cmmp"] = points["zero"] + .5 * (-1 * points["a1"] - points["a2"] + points["a3"]);
102 points["cmpp"] = points["zero"] + .5 * (-1 * points["a1"] + points["a2"] + points["a3"]);
103 points["cpmp"] = points["zero"] + .5 * (points["a1"] - points["a2"] + points["a3"]);
104 points["cppm"] = points["zero"] + .5 * (points["a1"] + points["a2"] - points["a3"]);
105 points["cppp"] = points["zero"] + .5 * (points["a1"] + points["a2"] + points["a3"]);
106 //get points from requested particle sets
107 int cshift = 1;
108 for (int i = 0; i < Psets.size(); i++)
109 {
110 ParticleSet& PS = *Psets[i];
111 for (int p = 0; p < PS.getTotalNum(); p++)
112 {
113 std::stringstream ss;
114 ss << p + cshift;
115 points[PS.getName() + ss.str()] = PS.R[p];
116 }
117 }
118 return true;
119 }
120
121
write_description(std::ostream & os,std::string & indent)122 void ReferencePoints::write_description(std::ostream& os, std::string& indent)
123 {
124 os << indent + "reference_points" << std::endl;
125 std::map<std::string, Point>::const_iterator it, end = points.end();
126 for (it = points.begin(); it != end; ++it)
127 {
128 os << indent + " " << it->first << ": " << it->second << std::endl;
129 }
130 os << indent + "end reference_points" << std::endl;
131 return;
132 }
133
save(std::vector<observable_helper * > & h5desc,hid_t gid) const134 void ReferencePoints::save(std::vector<observable_helper*>& h5desc, hid_t gid) const
135 {
136 observable_helper* oh = new observable_helper("reference_points");
137 oh->open(gid);
138 std::map<std::string, Point>::const_iterator it;
139 for (it = points.begin(); it != points.end(); ++it)
140 {
141 oh->addProperty(const_cast<Point&>(it->second), it->first);
142 }
143 h5desc.push_back(oh);
144 return;
145 }
146
147
148 } // namespace qmcplusplus
149