Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2025 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "AtomDistribution.h"
23 :
24 : namespace PLMD {
25 :
26 23 : std::unique_ptr<AtomDistribution> AtomDistribution::getAtomDistribution(std::string_view atomicDistr) {
27 23 : std::unique_ptr<AtomDistribution> distribution;
28 23 : if(atomicDistr == "line") {
29 1 : distribution = std::make_unique<theLine>();
30 22 : } else if (atomicDistr == "cube") {
31 1 : distribution = std::make_unique<uniformCube>();
32 21 : } else if (atomicDistr == "sphere") {
33 9 : distribution = std::make_unique<uniformSphere>();
34 12 : } else if (atomicDistr == "globs") {
35 1 : distribution = std::make_unique<twoGlobs>();
36 11 : } else if (atomicDistr == "sc") {
37 11 : distribution = std::make_unique<tiledSimpleCubic>();
38 : } else {
39 0 : plumed_error() << R"(The atomic distribution can be only "line", "cube", "sphere", "globs" and "sc", the input was ")"
40 0 : << atomicDistr <<'"';
41 : }
42 23 : return distribution;
43 0 : }
44 :
45 : class UniformSphericalVector {
46 : //double rminCub;
47 : double rCub;
48 :
49 : public:
50 : //assuming rmin=0
51 19 : explicit UniformSphericalVector(const double rmax):
52 19 : rCub (rmax*rmax*rmax/*-rminCub*/) {}
53 3800 : PLMD::Vector operator()(Random& rng) {
54 3800 : double rho = std::cbrt (/*rminCub + */rng.RandU01()*rCub);
55 3800 : double theta =std::acos (2.0*rng.RandU01() -1.0);
56 3800 : double phi = 2.0 * PLMD::pi * rng.RandU01();
57 : return Vector (
58 3800 : rho * sin (theta) * cos (phi),
59 3800 : rho * sin (theta) * sin (phi),
60 3800 : rho * cos (theta));
61 : }
62 : };
63 :
64 1 : void theLine::frame(std::vector<Vector>& posToUpdate,
65 : std::vector<double>& box,
66 : unsigned step,
67 : Random& rng) {
68 : auto nat = posToUpdate.size();
69 : UniformSphericalVector usv(0.5);
70 :
71 201 : for (unsigned i=0; i<nat; ++i) {
72 200 : posToUpdate[i] = Vector(i, 0, 0) + usv(rng);
73 : }
74 1 : box[0]=nat;
75 1 : box[1]=0.0;
76 1 : box[2]=0.0;
77 1 : box[3]=0.0;
78 1 : box[4]=1.75;
79 1 : box[5]=0.0;
80 1 : box[6]=0.0;
81 1 : box[7]=0.0;
82 1 : box[8]=1.75;
83 1 : }
84 :
85 17 : void uniformSphere::frame(std::vector<Vector>& posToUpdate,
86 : std::vector<double>& box,
87 : unsigned /*step*/,
88 : Random& rng) {
89 :
90 : //giving more or less a cubic udm of volume for each atom: V=nat
91 17 : const double rmax= std::cbrt ((3.0/(4.0*PLMD::pi)) * posToUpdate.size());
92 :
93 : UniformSphericalVector usv(rmax);
94 : auto s=posToUpdate.begin();
95 : auto e=posToUpdate.end();
96 : //I am using the iterators:this is slightly faster,
97 : // enough to overcome the cost of the vtable that I added
98 3417 : for (unsigned i=0; s!=e; ++s,++i) {
99 3400 : *s = usv (rng);
100 : }
101 :
102 17 : box[0]=2.0*rmax;
103 17 : box[1]=0.0;
104 17 : box[2]=0.0;
105 17 : box[3]=0.0;
106 17 : box[4]=2.0*rmax;
107 17 : box[5]=0.0;
108 17 : box[6]=0.0;
109 17 : box[7]=0.0;
110 17 : box[8]=2.0*rmax;
111 :
112 17 : }
113 :
114 1 : void twoGlobs::frame(std::vector<Vector>& posToUpdate,
115 : std::vector<double>& box,
116 : unsigned /*step*/,
117 : Random&rng) {
118 : //I am using two unigform spheres and 2V=n
119 1 : const double rmax= std::cbrt ((3.0/(8.0*PLMD::pi)) * posToUpdate.size());
120 :
121 : UniformSphericalVector usv(rmax);
122 1 : const std::array<const Vector,2> centers{
123 : PLMD::Vector{0.0,0.0,0.0},
124 : //so they do not overlap
125 : PLMD::Vector{2.0*rmax,2.0*rmax,2.0*rmax}
126 1 : };
127 1 : std::generate(posToUpdate.begin(),posToUpdate.end(),[&]() {
128 : //RandInt is only declared
129 : // return usv (rng) + centers[rng.RandInt(1)];
130 200 : return usv (rng) + centers[rng.RandU01()>0.5];
131 : });
132 :
133 1 : box[0]=4.0 *rmax;
134 1 : box[1]=0.0;
135 1 : box[2]=0.0;
136 1 : box[3]=0.0;
137 1 : box[4]=4.0 *rmax;
138 1 : box[5]=0.0;
139 1 : box[6]=0.0;
140 1 : box[7]=0.0;
141 1 : box[8]=4.0 *rmax;
142 1 : }
143 :
144 1 : void uniformCube::frame(std::vector<Vector>& posToUpdate,
145 : std::vector<double>& box,
146 : unsigned /*step*/,
147 : Random& rng) {
148 : //giving more or less a cubic udm of volume for each atom: V = nat
149 1 : const double rmax = std::cbrt(static_cast<double>(posToUpdate.size()));
150 :
151 :
152 :
153 : // std::generate(posToUpdate.begin(),posToUpdate.end(),[&]() {
154 : // return Vector (rndR(rng),rndR(rng),rndR(rng));
155 : // });
156 : auto s=posToUpdate.begin();
157 : auto e=posToUpdate.end();
158 : //I am using the iterators:this is slightly faster,
159 : // enough to overcome the cost of the vtable that I added
160 201 : for (unsigned i=0; s!=e; ++s,++i) {
161 200 : *s = Vector (rng.RandU01()*rmax,rng.RandU01()*rmax,rng.RandU01()*rmax);
162 : }
163 : //+0.05 to avoid overlap
164 1 : box[0]=rmax+0.05;
165 1 : box[1]=0.0;
166 1 : box[2]=0.0;
167 1 : box[3]=0.0;
168 1 : box[4]=rmax+0.05;
169 1 : box[5]=0.0;
170 1 : box[6]=0.0;
171 1 : box[7]=0.0;
172 1 : box[8]=rmax+0.05;
173 :
174 1 : }
175 :
176 15 : void tiledSimpleCubic::frame(std::vector<Vector>& posToUpdate,
177 : std::vector<double>& box,
178 : unsigned /*step*/,
179 : Random& rng) {
180 : //Tiling the space in this way will not tests 100% the pbc, but
181 : //I do not think that write a spacefilling curve, like Hilbert, Peano or Morton
182 : //could be a good idea, in this case
183 15 : const unsigned rmax = std::ceil(std::cbrt(static_cast<double>(posToUpdate.size())));
184 :
185 : auto s=posToUpdate.begin();
186 : auto e=posToUpdate.end();
187 : //I am using the iterators:this is slightly faster,
188 : // enough to overcome the cost of the vtable that I added
189 91 : for (unsigned k=0; k<rmax&&s!=e; ++k) {
190 460 : for (unsigned j=0; j<rmax&&s!=e; ++j) {
191 2330 : for (unsigned i=0; i<rmax&&s!=e; ++i) {
192 1946 : *s = Vector (i,j,k);
193 : ++s;
194 : }
195 : }
196 : }
197 15 : box[0]=rmax;
198 15 : box[1]=0.0;
199 15 : box[2]=0.0;
200 15 : box[3]=0.0;
201 15 : box[4]=rmax;
202 15 : box[5]=0.0;
203 15 : box[6]=0.0;
204 15 : box[7]=0.0;
205 15 : box[8]=rmax;
206 :
207 15 : }
208 :
209 0 : void fileTraj::rewind() {
210 0 : auto errormessage=parser.rewind();
211 0 : if (errormessage) {
212 : // A workarounf for not implemented rewind is to dump the trajectory in an xyz and then read that
213 0 : plumed_error()<<*errormessage;
214 : }
215 : //the extra false prevents an infinite loop in case of unexpected consecutice EOFs after a rewind
216 0 : step(false);
217 0 : }
218 :
219 : //read the next step
220 0 : void fileTraj::step(bool doRewind) {
221 0 : read=false;
222 0 : long long int mystep=0;
223 : double timeStep;
224 0 : std::optional<std::string> errormessage;
225 0 : if (masses.empty()) {
226 0 : errormessage=parser.readHeader(
227 : mystep,
228 : timeStep
229 : );
230 0 : if (errormessage) {
231 0 : plumed_error()<<*errormessage;
232 : }
233 0 : const size_t natoms = parser.nOfAtoms();
234 :
235 0 : masses.assign(natoms,0.0);
236 0 : charges.assign(natoms,0.0);
237 0 : coordinates.assign(natoms,Vector(0.0,0.0,0.0));
238 0 : cell.assign(9,0.0);
239 0 : errormessage=parser.readAtoms(1,
240 0 : dont_read_pbc,
241 : false,
242 : 0,
243 : 0,
244 : mystep,
245 : masses.data(),
246 : charges.data(),
247 : &coordinates[0][0],
248 : cell.data()
249 : );
250 : } else {
251 0 : errormessage=parser.readFrame(1,
252 0 : dont_read_pbc,
253 : false,
254 : 0,
255 : 0,
256 : mystep,
257 : timeStep,
258 : masses.data(),
259 : charges.data(),
260 : &coordinates[0][0],
261 : cell.data()
262 : );
263 : }
264 :
265 0 : if (errormessage) {
266 0 : if (*errormessage =="EOF" && doRewind) {
267 0 : rewind();
268 : } else {
269 0 : plumed_error()<<*errormessage;
270 : }
271 : }
272 0 : }
273 :
274 0 : void fileTraj::frame(std::vector<Vector>& posToUpdate,
275 : std::vector<double>& box,
276 : unsigned /*step*/,
277 : Random& /*rng*/) {
278 0 : if (read) {
279 0 : step();
280 : }
281 0 : read=true;
282 0 : std::copy(coordinates.begin(),coordinates.end(),posToUpdate.begin());
283 0 : std::copy(cell.begin(),cell.end(),box.begin());
284 0 : }
285 :
286 0 : fileTraj::fileTraj(std::string_view fmt,
287 : std::string_view fname,
288 : bool useMolfile,
289 0 : int command_line_natoms) {
290 0 : parser.init(fmt,
291 : fname,
292 : useMolfile,
293 : command_line_natoms);
294 0 : step();
295 0 : }
296 :
297 0 : bool fileTraj::overrideNat(unsigned& natoms) {
298 0 : natoms = masses.size();
299 0 : return true;
300 : }
301 :
302 12 : repliedTrajectory::repliedTrajectory(std::unique_ptr<AtomDistribution>&& d,
303 : const unsigned repeatX,
304 : const unsigned repeatY,
305 : const unsigned repeatZ,
306 12 : const unsigned nat):
307 12 : distribution(std::move(d)),
308 12 : rX(repeatX),
309 12 : rY(repeatY),
310 12 : rZ(repeatZ),
311 12 : coordinates(nat,Vector{})
312 12 : {}
313 :
314 12 : void repliedTrajectory::frame(std::vector<Vector>& posToUpdate, std::vector<double>& box,
315 : unsigned step,
316 : Random& rng) {
317 12 : distribution->frame(coordinates,box,step,rng);
318 :
319 : // repetitions
320 : auto p = posToUpdate.begin();
321 : const auto nat = coordinates.size();
322 :
323 : assert((rX*rY*rZ)*nat == posToUpdate.size());
324 :
325 12 : Vector boxX(box[0],box[1],box[2]);
326 12 : Vector boxY(box[3],box[4],box[5]);
327 12 : Vector boxZ(box[6],box[7],box[8]);
328 33 : for (unsigned x=0; x<rX; ++x) {
329 59 : for (unsigned y=0; y<rY; ++y) {
330 162 : for (unsigned z=0; z<rZ; ++z) {
331 24624 : for (unsigned i=0; i<nat; ++i) {
332 24500 : *p=coordinates[i]
333 : + x * boxX
334 : + y * boxY
335 : + z * boxZ;
336 : ++p;
337 : }
338 : }
339 : }
340 : }
341 12 : box[0]*=rX;
342 12 : box[1]*=rX;
343 12 : box[2]*=rX;
344 :
345 12 : box[3]*=rY;
346 12 : box[4]*=rY;
347 12 : box[5]*=rY;
348 :
349 12 : box[6]*=rZ;
350 12 : box[7]*=rZ;
351 12 : box[8]*=rZ;
352 12 : }
353 :
354 12 : bool repliedTrajectory::overrideNat(unsigned& natoms) {
355 12 : natoms = (rX*rY*rZ)*coordinates.size();
356 12 : return true;
357 : }
358 : } //namespace PLMD
|
