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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
9 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
10 //
11 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
12 //////////////////////////////////////////////////////////////////////////////////////
13
14
15 /**@file ParticleSet.BC.cpp
16 * @brief definition of functions controlling Boundary Conditions
17 */
18 #include "Particle/ParticleSet.h"
19 #include "Particle/FastParticleOperators.h"
20 #include "Message/OpenMP.h"
21 #include "LongRange/StructFact.h"
22
23 namespace qmcplusplus
24 {
25 /** Creating StructureFactor
26 *
27 * Currently testing only 1 component for PBCs.
28 */
createSK()29 void ParticleSet::createSK()
30 {
31 //if(!sorted_ids && !reordered_ids)
32 //{
33 // //save ID and GroupID
34 // orgID=ID;
35 // orgGroupID=GroupID;
36 // if(groups()<1)
37 // {
38 // int nspecies=mySpecies.getTotalNum();
39 // std::vector<int> ppg(nspecies,0);
40 // for(int iat=0; iat<GroupID.size(); ++iat) ppg[GroupID[iat]]+=1;
41 // SubPtcl.resize(nspecies+1);
42 // SubPtcl[0]=0;
43 // for(int i=0; i<nspecies; ++i) SubPtcl[i+1]=SubPtcl[i]+ppg[i];
44 // int new_id=0;
45 // for(int i=0; i<nspecies; ++i)
46 // for(int iat=0; iat<GroupID.size(); ++iat) if(GroupID[iat]==i) orgID[new_id++]=ID[iat];
47 // bool grouped=true;
48 // for(int iat=0; iat<ID.size(); ++iat) grouped &= (orgID[iat]==ID[iat]);
49 // if(grouped)
50 // {
51 // app_log() << " ParticleSet is grouped. No need to reorder." << std::endl;
52 // }
53 // else
54 // {
55 // app_log() << " Need to reorder. Only R is swapped." << std::endl;
56 // ParticlePos_t oldR(R);
57 // for(int iat=0; iat<R.size(); ++iat) R[iat]=oldR[orgID[iat]];
58 // for(int i=0; i<groups(); ++i)
59 // for(int iat=first(i); iat<last(i); ++iat) GroupID[iat]=i;
60 // reordered_ids=true;
61 // }
62 // }//once group is set, nothing to be done
63 // sorted_ids=true;
64 //}
65 //int membersize= mySpecies.addAttribute("membersize");
66 //for(int ig=0; ig<mySpecies.size(); ++ig)
67 // SubPtcl[ig+1]=SubPtcl[ig]+mySpecies(membersize,ig);
68
69 if (Lattice.explicitly_defined)
70 convert2Cart(R); //make sure that R is in Cartesian coordinates
71
72 if (Lattice.SuperCellEnum != SUPERCELL_OPEN)
73 {
74 Lattice.SetLRCutoffs(Lattice.Rv);
75 LRBox = Lattice;
76 bool changed = false;
77 if (Lattice.SuperCellEnum == SUPERCELL_SLAB && Lattice.VacuumScale != 1.0)
78 {
79 LRBox.R(2, 0) *= Lattice.VacuumScale;
80 LRBox.R(2, 1) *= Lattice.VacuumScale;
81 LRBox.R(2, 2) *= Lattice.VacuumScale;
82 changed = true;
83 }
84 else if (Lattice.SuperCellEnum == SUPERCELL_WIRE && Lattice.VacuumScale != 1.0)
85 {
86 LRBox.R(1, 0) *= Lattice.VacuumScale;
87 LRBox.R(1, 1) *= Lattice.VacuumScale;
88 LRBox.R(1, 2) *= Lattice.VacuumScale;
89 LRBox.R(2, 0) *= Lattice.VacuumScale;
90 LRBox.R(2, 1) *= Lattice.VacuumScale;
91 LRBox.R(2, 2) *= Lattice.VacuumScale;
92 changed = true;
93 }
94 LRBox.reset();
95 LRBox.SetLRCutoffs(LRBox.Rv);
96 LRBox.printCutoffs(app_log());
97
98 if (changed)
99 {
100 app_summary() << " Simulation box changed by vacuum supercell conditions" << std::endl;
101 app_log() << "--------------------------------------- " << std::endl;
102 LRBox.print(app_log());
103 app_log() << "--------------------------------------- " << std::endl;
104 }
105
106 if (SK)
107 {
108 app_log() << "\n Structure Factor is reset by " << Lattice.LR_kc << std::endl;
109 SK->UpdateNewCell(*this, LRBox.LR_kc);
110 }
111 else
112 {
113 app_log() << "\n Creating Structure Factor for periodic systems " << LRBox.LR_kc << std::endl;
114 SK = std::make_unique<StructFact>(*this, LRBox.LR_kc);
115 }
116 //Lattice.print(app_log());
117 //This uses the copy constructor to avoid recomputing the data.
118 //SKOld = new StructFact(*SK);
119 }
120 //set the mass array
121 int beforemass = mySpecies.numAttributes();
122 int massind = mySpecies.addAttribute("mass");
123 if (beforemass == massind)
124 {
125 app_log() << " ParticleSet::createSK setting mass of " << getName() << " to 1.0" << std::endl;
126 for (int ig = 0; ig < mySpecies.getTotalNum(); ++ig)
127 mySpecies(massind, ig) = 1.0;
128 }
129 for (int iat = 0; iat < GroupID.size(); iat++)
130 Mass[iat] = mySpecies(massind, GroupID[iat]);
131
132 coordinates_->setAllParticlePos(R);
133 }
134
turnOnPerParticleSK()135 void ParticleSet::turnOnPerParticleSK()
136 {
137 if (SK)
138 SK->turnOnStorePerParticle(*this);
139 else
140 APP_ABORT(
141 "ParticleSet::turnOnPerParticleSK trying to turn on per particle storage in SK but SK has not been created.");
142 }
143
convert(const ParticlePos_t & pin,ParticlePos_t & pout)144 void ParticleSet::convert(const ParticlePos_t& pin, ParticlePos_t& pout)
145 {
146 if (pin.getUnit() == pout.getUnit())
147 {
148 pout = pin;
149 return;
150 }
151 if (pin.getUnit() == PosUnit::Lattice)
152 //convert to CartesianUnit
153 {
154 ConvertPosUnit<ParticlePos_t, Tensor_t, DIM>::apply(pin, Lattice.R, pout, 0, pin.size());
155 }
156 else
157 //convert to LatticeUnit
158 {
159 ConvertPosUnit<ParticlePos_t, Tensor_t, DIM>::apply(pin, Lattice.G, pout, 0, pin.size());
160 }
161 }
162
convert2Unit(const ParticlePos_t & pin,ParticlePos_t & pout)163 void ParticleSet::convert2Unit(const ParticlePos_t& pin, ParticlePos_t& pout)
164 {
165 pout.setUnit(PosUnit::Lattice);
166 if (pin.getUnit() == PosUnit::Lattice)
167 pout = pin;
168 else
169 ConvertPosUnit<ParticlePos_t, Tensor_t, DIM>::apply(pin, Lattice.G, pout, 0, pin.size());
170 }
171
convert2Cart(const ParticlePos_t & pin,ParticlePos_t & pout)172 void ParticleSet::convert2Cart(const ParticlePos_t& pin, ParticlePos_t& pout)
173 {
174 pout.setUnit(PosUnit::Cartesian);
175 if (pin.getUnit() == PosUnit::Cartesian)
176 pout = pin;
177 else
178 ConvertPosUnit<ParticlePos_t, Tensor_t, DIM>::apply(pin, Lattice.R, pout, 0, pin.size());
179 }
180
convert2Unit(ParticlePos_t & pinout)181 void ParticleSet::convert2Unit(ParticlePos_t& pinout)
182 {
183 if (pinout.getUnit() == PosUnit::Lattice)
184 return;
185 else
186 {
187 pinout.setUnit(PosUnit::Lattice);
188 ConvertPosUnit<ParticlePos_t, Tensor_t, DIM>::apply(pinout, Lattice.G, 0, pinout.size());
189 }
190 }
191
convert2Cart(ParticlePos_t & pinout)192 void ParticleSet::convert2Cart(ParticlePos_t& pinout)
193 {
194 if (pinout.getUnit() == PosUnit::Cartesian)
195 return;
196 else
197 {
198 pinout.setUnit(PosUnit::Cartesian);
199 ConvertPosUnit<ParticlePos_t, Tensor_t, DIM>::apply(pinout, Lattice.R, 0, pinout.size());
200 }
201 }
202
applyBC(const ParticlePos_t & pin,ParticlePos_t & pout)203 void ParticleSet::applyBC(const ParticlePos_t& pin, ParticlePos_t& pout) { applyBC(pin, pout, 0, pin.size()); }
204
applyBC(const ParticlePos_t & pin,ParticlePos_t & pout,int first,int last)205 void ParticleSet::applyBC(const ParticlePos_t& pin, ParticlePos_t& pout, int first, int last)
206 {
207 if (pin.getUnit() == PosUnit::Cartesian)
208 {
209 if (pout.getUnit() == PosUnit::Cartesian)
210 ApplyBConds<ParticlePos_t, Tensor_t, DIM>::Cart2Cart(pin, Lattice.G, Lattice.R, pout, first, last);
211 else if (pout.getUnit() == PosUnit::Lattice)
212 ApplyBConds<ParticlePos_t, Tensor_t, DIM>::Cart2Unit(pin, Lattice.G, pout, first, last);
213 else
214 throw std::runtime_error("Unknown unit conversion");
215 }
216 else if (pin.getUnit() == PosUnit::Lattice)
217 {
218 if (pout.getUnit() == PosUnit::Cartesian)
219 ApplyBConds<ParticlePos_t, Tensor_t, DIM>::Unit2Cart(pin, Lattice.R, pout, first, last);
220 else if (pout.getUnit() == PosUnit::Lattice)
221 ApplyBConds<ParticlePos_t, Tensor_t, DIM>::Unit2Unit(pin, pout, first, last);
222 else
223 throw std::runtime_error("Unknown unit conversion");
224 }
225 else
226 throw std::runtime_error("Unknown unit conversion");
227 }
228
applyBC(ParticlePos_t & pos)229 void ParticleSet::applyBC(ParticlePos_t& pos)
230 {
231 if (pos.getUnit() == PosUnit::Lattice)
232 {
233 ApplyBConds<ParticlePos_t, Tensor_t, DIM>::Unit2Unit(pos, 0, TotalNum);
234 }
235 else
236 {
237 ApplyBConds<ParticlePos_t, Tensor_t, DIM>::Cart2Cart(pos, Lattice.G, Lattice.R, 0, TotalNum);
238 }
239 }
240
applyMinimumImage(ParticlePos_t & pinout)241 void ParticleSet::applyMinimumImage(ParticlePos_t& pinout)
242 {
243 if (Lattice.SuperCellEnum == SUPERCELL_OPEN)
244 return;
245 for (int i = 0; i < pinout.size(); ++i)
246 Lattice.applyMinimumImage(pinout[i]);
247 }
248
convert2UnitInBox(const ParticlePos_t & pin,ParticlePos_t & pout)249 void ParticleSet::convert2UnitInBox(const ParticlePos_t& pin, ParticlePos_t& pout)
250 {
251 pout.setUnit(PosUnit::Lattice);
252 convert2Unit(pin, pout); // convert to crystalline unit
253 put2box(pout);
254 }
255
convert2CartInBox(const ParticlePos_t & pin,ParticlePos_t & pout)256 void ParticleSet::convert2CartInBox(const ParticlePos_t& pin, ParticlePos_t& pout)
257 {
258 convert2UnitInBox(pin, pout); // convert to crystalline unit
259 convert2Cart(pout);
260 }
261 } // namespace qmcplusplus
262