moldyn-loader.h 12 KB
Newer Older
Tobias WEBER's avatar
Tobias WEBER committed
1
2
3
4
5
/**
 * loads atom dynamics file
 * @author Tobias Weber <tweber@ill.fr>
 * @date Dec-2019
 * @license GPLv3, see 'LICENSE' file
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 *
 * ----------------------------------------------------------------------------
 * mag-core (part of the Takin software suite)
 * Copyright (C) 2018-2021  Tobias WEBER (Institut Laue-Langevin (ILL),
 *                          Grenoble, France).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3 of the License.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 * ----------------------------------------------------------------------------
Tobias WEBER's avatar
Tobias WEBER committed
24
25
26
27
28
29
30
31
32
33
34
 */

#ifndef __MOLDYN_H__
#define __MOLDYN_H__


#include <fstream>
#include <iostream>
#include <vector>
#include <string>

35
36
37
// for progress callback
#include <boost/signals2/signal.hpp>

38
39
#include "tlibs2/libs/file.h"
#include "tlibs2/libs/str.h"
Tobias WEBER's avatar
Tobias WEBER committed
40
41
42
43
44
45
46
47
48
49
50
51
52


template<class t_real, class t_vec>
class MolFrame
{
	public:
		MolFrame()
		{
			m_config.reserve(128);
		}


		void AddAtomConfig(std::vector<t_vec>&& config)
Tobias WEBER's avatar
Tobias WEBER committed
53
		{ m_config.emplace_back(std::move(config)); }
Tobias WEBER's avatar
Tobias WEBER committed
54
55

		void AddAtomConfig(const std::vector<t_vec>& config)
Tobias WEBER's avatar
Tobias WEBER committed
56
		{ m_config.push_back(config); }
Tobias WEBER's avatar
Tobias WEBER committed
57

Tobias WEBER's avatar
Tobias WEBER committed
58

Tobias WEBER's avatar
Tobias WEBER committed
59
		std::size_t GetNumAtomTypes() const
Tobias WEBER's avatar
Tobias WEBER committed
60
		{ return m_config.size(); }
Tobias WEBER's avatar
Tobias WEBER committed
61

Tobias WEBER's avatar
Tobias WEBER committed
62
63
64

		const std::vector<t_vec>& GetCoords(std::size_t idxType) const
		{ return m_config[idxType]; }
Tobias WEBER's avatar
Tobias WEBER committed
65

Tobias WEBER's avatar
Tobias WEBER committed
66

Tobias WEBER's avatar
Tobias WEBER committed
67
68
69
70
71
72
73
74
75
		/**
		 * get the atom coordinates
		 */
		const t_vec& GetAtomCoords(std::size_t idxType, std::size_t idxSubType) const
		{
			return m_config[idxType][idxSubType];
		}


Tobias WEBER's avatar
Tobias WEBER committed
76
77
78
79
80
81
82
83
84
		/**
		 * removes one atoms of type idxType and index idxSubType
		 */
		void RemoveAtom(std::size_t idxType, std::size_t idxSubType)
		{
			m_config[idxType].erase(m_config[idxType].begin() + idxSubType);
		}


Tobias WEBER's avatar
Tobias WEBER committed
85
86
87
88
89
90
91
92
93
		/**
		 * removes all atoms at index idx
		 */
		void RemoveAtoms(std::size_t idx)
		{
			m_config.erase(m_config.begin()+idx);
		}


Tobias WEBER's avatar
Tobias WEBER committed
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
	private:
		// atoms -> coordinates
		std::vector<std::vector<t_vec>> m_config;
};



template<class t_real, class t_vec>
class MolDyn
{
	public:
		MolDyn() : m_baseA(3), m_baseB(3), m_baseC(3)
		{
			m_frames.reserve(16384);
		}

		void SetBaseA(t_real x, t_real y, t_real z)
		{
			m_baseA[0] = x;
			m_baseA[1] = y;
			m_baseA[2] = z;
		}

		void SetBaseB(t_real x, t_real y, t_real z)
		{
			m_baseB[0] = x;
			m_baseB[1] = y;
			m_baseB[2] = z;
		}

		void SetBaseC(t_real x, t_real y, t_real z)
		{
			m_baseC[0] = x;
			m_baseC[1] = y;
			m_baseC[2] = z;
		}


Tobias WEBER's avatar
Tobias WEBER committed
132
133
134
135
136
137
138
		const t_vec& GetBaseA() const { return m_baseA; }
		const t_vec& GetBaseB() const { return m_baseB; }
		const t_vec& GetBaseC() const { return m_baseC; }


		std::size_t GetFrameCount() const
		{ return m_frames.size(); }
139

Tobias WEBER's avatar
Tobias WEBER committed
140
141
142
143
		const MolFrame<t_real, t_vec>& GetFrame(std::size_t frame) const
		{ return m_frames[frame]; }


Tobias WEBER's avatar
Tobias WEBER committed
144
		std::size_t GetNumAtomTypes() const
Tobias WEBER's avatar
Tobias WEBER committed
145
146
		{ return m_vecAtoms.size(); }

Tobias WEBER's avatar
Tobias WEBER committed
147
148
149
150
151
152
153
154
155
156
157
158
		std::size_t GetNumAtomsTotal() const
		{
			std::size_t num = 0;

			for(std::size_t numpertype : m_vecAtomNums)
				num += numpertype;

			return num;
		}

		const std::string& GetAtomName(std::size_t idxType) const
		{ return m_vecAtoms[idxType]; }
Tobias WEBER's avatar
Tobias WEBER committed
159

Tobias WEBER's avatar
Tobias WEBER committed
160
161
		unsigned int GetAtomNum(std::size_t idxType) const
		{ return m_vecAtomNums[idxType]; }
Tobias WEBER's avatar
Tobias WEBER committed
162

Tobias WEBER's avatar
Tobias WEBER committed
163

Tobias WEBER's avatar
Tobias WEBER committed
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
		void AddAtomType(const std::string& name, unsigned int number)
		{
			m_vecAtoms.push_back(name);
			m_vecAtomNums.push_back(number);
		}


		void AddFrame(MolFrame<t_real, t_vec>&& frame)
		{
			m_frames.emplace_back(std::move(frame));
		}

		void AddFrame(const MolFrame<t_real, t_vec>& frame)
		{
			m_frames.push_back(frame);
		}


Tobias WEBER's avatar
Tobias WEBER committed
182
183
184
185
186
187
188
189
190
		/**
		 * get atom coordinates for a specific frame
		 */
		const t_vec& GetAtomCoords(std::size_t idxType, std::size_t idxSubType, std::size_t iFrameIdx) const
		{
			return GetFrame(iFrameIdx).GetAtomCoords(idxType, idxSubType);
		}


Tobias WEBER's avatar
Tobias WEBER committed
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
		/**
		 * get atom coordinates for all frames
		 */
		std::vector<t_vec> GetAtomCoords(std::size_t idxType, std::size_t idxSubType) const
		{
			std::vector<t_vec> allcoords;
			allcoords.reserve(m_frames.size());

			if(!m_vecAtomNums[idxType])
				return allcoords;

			for(const MolFrame<t_real, t_vec>& frame : m_frames)
				allcoords.push_back(frame.GetAtomCoords(idxType, idxSubType));

			return allcoords;
		}


Tobias WEBER's avatar
Tobias WEBER committed
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
		/**
		 * removes one atoms of type idxType and index idxSubType
		 */
		void RemoveAtom(std::size_t idxType, std::size_t idxSubType)
		{
			if(!m_vecAtomNums[idxType])
				return;

			for(MolFrame<t_real, t_vec>& frame : m_frames)
				frame.RemoveAtom(idxType, idxSubType);

			--m_vecAtomNums[idxType];
		}


Tobias WEBER's avatar
Tobias WEBER committed
224
225
226
227
228
229
230
231
232
233
234
235
		/**
		 * removes all atoms at index idx
		 */
		void RemoveAtoms(std::size_t idx)
		{
			m_vecAtoms.erase(m_vecAtoms.begin() + idx);
			m_vecAtomNums.erase(m_vecAtomNums.begin() + idx);

			for(MolFrame<t_real, t_vec>& frame : m_frames)
				frame.RemoveAtoms(idx);
		}

Tobias WEBER's avatar
Tobias WEBER committed
236

Tobias WEBER's avatar
Tobias WEBER committed
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
		void Clear()
		{
			m_vecAtoms.clear();
			m_vecAtomNums.clear();
			m_frames.clear();

			m_sigLoadProgress.disconnect_all_slots();
			m_sigSaveProgress.disconnect_all_slots();
		}



		/**
		 * loading of files
		 */
Tobias WEBER's avatar
Tobias WEBER committed
252
253
254
255
256
257
258
		bool LoadFile(const std::string& filename, unsigned int frameskip = 0)
		{
			const std::string strDelim{" \t"};

			std::ifstream ifstr{filename};
			if(!ifstr)
			{
Tobias WEBER's avatar
Tobias WEBER committed
259
				std::cerr << "Cannot open \"" << filename << "\" for loading.";
Tobias WEBER's avatar
Tobias WEBER committed
260
261
262
				return 0;
			}

Tobias WEBER's avatar
Tobias WEBER committed
263

Tobias WEBER's avatar
Tobias WEBER committed
264
265
266
			std::size_t filesize = tl2::get_file_size(ifstr);
			std::cout << "File size: " << filesize / 1024 / 1024 << " MB." << std::endl;

Tobias WEBER's avatar
Tobias WEBER committed
267
268
269
			std::getline(ifstr, m_strSys);
			tl2::trim(m_strSys);
			std::cout << "System: " << m_strSys << std::endl;
Tobias WEBER's avatar
Tobias WEBER committed
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324



			std::string strScale;
			std::getline(ifstr, strScale);
			t_real scale = tl2::str_to_var<t_real>(strScale);
			std::cout << "scale: " << scale << std::endl;



			std::string strVecs1;
			std::string strVecs2;
			std::string strVecs3;
			std::getline(ifstr, strVecs1);
			std::getline(ifstr, strVecs2);
			std::getline(ifstr, strVecs3);

			std::vector<t_real> _vecBase1, _vecBase2, _vecBase3;
			tl2::get_tokens<t_real>(strVecs1, strDelim, _vecBase1);
			tl2::get_tokens<t_real>(strVecs2, strDelim, _vecBase2);
			tl2::get_tokens<t_real>(strVecs3, strDelim, _vecBase3);

			if(_vecBase1.size()!=3 || _vecBase2.size()!=3 || _vecBase3.size()!=3)
			{
				std::cerr << "Invalid base vectors." << std::endl;
				return 0;
			}


			SetBaseA(_vecBase1[0]*scale, _vecBase2[0]*scale, _vecBase3[0]*scale);
			SetBaseB(_vecBase1[1]*scale, _vecBase2[1]*scale, _vecBase3[1]*scale);
			SetBaseC(_vecBase1[2]*scale, _vecBase2[2]*scale, _vecBase3[2]*scale);


			std::string strAtoms;
			std::string strAtomNums;
			std::getline(ifstr, strAtoms);
			std::getline(ifstr, strAtomNums);
			tl2::get_tokens<std::string>(strAtoms, strDelim, m_vecAtoms);
			tl2::get_tokens<unsigned int>(strAtomNums, strDelim, m_vecAtomNums);

			if(m_vecAtoms.size() != m_vecAtomNums.size())
			{
				std::cerr << "Atom size mismatch." << std::endl;
				return 0;
			}

			for(std::size_t i=0; i<m_vecAtoms.size(); ++i)
			{
				std::cout << m_vecAtomNums[i] << " " << m_vecAtoms[i] << " atoms." << std::endl;
			}



			std::size_t iNumConfigs = 0;
Tobias WEBER's avatar
Tobias WEBER committed
325
			t_real percentage = 0;
Tobias WEBER's avatar
Tobias WEBER committed
326
327
328
329
330
331
332
333
334
335
336
			while(true)
			{
				std::string strConfig;
				std::getline(ifstr, strConfig);
				tl2::trim(strConfig);

				if(ifstr.eof())
					break;

				if(frameskip || iNumConfigs % 100)
				{
Tobias WEBER's avatar
Tobias WEBER committed
337
					std::cout << "\rReading " << strConfig << ". "
Tobias WEBER's avatar
Tobias WEBER committed
338
						<< static_cast<unsigned>(percentage) << " %.                ";
Tobias WEBER's avatar
Tobias WEBER committed
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
					std::cout.flush();
				}


				MolFrame<t_real, t_vec> frame;
				for(std::size_t iAtomType=0; iAtomType<m_vecAtoms.size(); ++iAtomType)
				{
					std::vector<t_vec> atomconf;

					for(std::size_t iAtom=0; iAtom<m_vecAtomNums[iAtomType]; ++iAtom)
					{
						std::string strCoords;
						std::getline(ifstr, strCoords);

						t_vec vecCoords;
						vecCoords.reserve(3);
						tl2::get_tokens<t_real>(strCoords, strDelim, vecCoords);


						if(vecCoords.size() != 3)
						{
							std::cerr << "Invalid coordinate." << std::endl;
							return 0;
						}

Tobias WEBER's avatar
Tobias WEBER committed
364
365
366
367
368
						// center cell (in rlu)
						vecCoords[0] -= 0.5;
						vecCoords[1] -= 0.5;
						vecCoords[2] -= 0.5;

Tobias WEBER's avatar
Tobias WEBER committed
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
						atomconf.emplace_back(std::move(vecCoords));
					}

					frame.AddAtomConfig(std::move(atomconf));
				}

				AddFrame(std::move(frame));
				++iNumConfigs;


				// skip frames
				for(unsigned int skipped=0; skipped<frameskip; ++skipped)
				{
					std::string strTmp;
					std::getline(ifstr, strTmp);
					//std::cout << "Skipping " << strTmp << "..." << "        ";

					for(std::size_t iAtomType=0; iAtomType<m_vecAtoms.size(); ++iAtomType)
					{
						for(std::size_t iAtom=0; iAtom<m_vecAtomNums[iAtomType]; ++iAtom)
							std::getline(ifstr, strTmp);
					}
				}

Tobias WEBER's avatar
Tobias WEBER committed
393
394

				std::size_t filepos = tl2::get_file_pos(ifstr);
Tobias WEBER's avatar
Tobias WEBER committed
395
				percentage = static_cast<t_real>(filepos*100) / static_cast<t_real>(filesize);
Tobias WEBER's avatar
Tobias WEBER committed
396
				if(m_sigLoadProgress.num_slots() && !*m_sigLoadProgress(percentage))
397
398
399
400
				{
					std::cerr << "\nLoading cancelled." << std::endl;
					return 0;
				}
Tobias WEBER's avatar
Tobias WEBER committed
401
			}
402

Tobias WEBER's avatar
Tobias WEBER committed
403
			std::cout << "\rRead " << iNumConfigs << " configurations. " << "                        " << std::endl;
Tobias WEBER's avatar
Tobias WEBER committed
404
405
406
407
			return 1;
		}


Tobias WEBER's avatar
Tobias WEBER committed
408
409
410
411
412

		/**
		 * saving of files
		 */
		bool SaveFile(const std::string& filename)
413
		{
Tobias WEBER's avatar
Tobias WEBER committed
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
			std::cout << m_sigSaveProgress.num_slots() << std::endl;
			std::ofstream ofstr{filename};
			if(!ofstr)
			{
				std::cerr << "Cannot open \"" << filename << "\" for saving.";
				return 0;
			}

			ofstr.precision(8);

			if(m_baseA.size() != 3)
				return 0;

			ofstr << m_strSys << "\n" << "1" << std::endl;
			for(int i=0; i<m_baseA.size(); ++i)
				ofstr << m_baseA[i] << " " << m_baseB[i] << " " << m_baseC[i] << std::endl;
430

Tobias WEBER's avatar
Tobias WEBER committed
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
			for(const std::string& strAtom : m_vecAtoms)
				ofstr << strAtom << " ";
			ofstr << std::endl;

			for(unsigned int numAtom : m_vecAtomNums)
				ofstr << numAtom << " ";
			ofstr << std::endl;

			// iterate frames
			t_real percentage = 0;
			for(std::size_t frame=0; frame<m_frames.size(); ++frame)
			{
				ofstr << "Config " << (frame+1) << "\n";
				const MolFrame<t_real, t_vec>& config = m_frames[frame];

				// iterate atom types
Tobias WEBER's avatar
Tobias WEBER committed
447
				for(std::size_t atomidx=0; atomidx<config.GetNumAtomTypes(); ++atomidx)
Tobias WEBER's avatar
Tobias WEBER committed
448
449
450
451
452
453
454
455
456
457
458
				{
					const auto& coords = config.GetCoords(atomidx);
					// iterate coordinates
					for(const auto& vec : coords)
					{
						if(vec.size() != 3)
							return 0;
						ofstr << vec[0]+0.5 << " " << vec[1]+0.5 << " " << vec[2]+0.5 << "\n";
					}
				}

Tobias WEBER's avatar
Tobias WEBER committed
459
				percentage = static_cast<t_real>((frame+1)*100)/static_cast<t_real>(m_frames.size());
Tobias WEBER's avatar
Tobias WEBER committed
460
461
462
463
464
465
466
467
468
				if(m_sigSaveProgress.num_slots() && !*m_sigSaveProgress(percentage))
				{
					std::cerr << "\nSaving cancelled." << std::endl;
					return 0;
				}

				if(frame % 100)
				{
					std::cout << "\rSaving configuration " << (frame+1) << " of " << m_frames.size() << ". "
Tobias WEBER's avatar
Tobias WEBER committed
469
						<< static_cast<unsigned>(percentage) << " %.                ";
Tobias WEBER's avatar
Tobias WEBER committed
470
471
472
473
474
475
476
					std::cout.flush();
				}
			}

			std::cout << "\rSaved " << m_frames.size() << " configurations. " << "                        " << std::endl;
			ofstr.flush();
			return 1;
477
478
479
		}


Tobias WEBER's avatar
Tobias WEBER committed
480
481
482
483
484
485
486
487
488

		template<class subscriber> void SubscribeToLoadProgress(const subscriber& subs)
		{ m_sigLoadProgress.connect(subs); }

		template<class subscriber> void SubscribeToSaveProgress(const subscriber& subs)
		{ m_sigSaveProgress.connect(subs); }


		template<class subscriber> void UnsubscribeFromLoadProgress(const subscriber* subs=nullptr)
Tobias WEBER's avatar
Tobias WEBER committed
489
		{
Tobias WEBER's avatar
Tobias WEBER committed
490
491
492
493
494
495
496
497
498
499
			// TODO
			//if(!subs)
				m_sigLoadProgress.disconnect_all_slots();
		}

		template<class subscriber> void UnsubscribeFromSaveProgress(const subscriber* subs=nullptr)
		{
			// TODO
			//if(!subs)
				m_sigSaveProgress.disconnect_all_slots();
Tobias WEBER's avatar
Tobias WEBER committed
500
501
502
		}


Tobias WEBER's avatar
Tobias WEBER committed
503
	private:
Tobias WEBER's avatar
Tobias WEBER committed
504
505
		std::string m_strSys;

Tobias WEBER's avatar
Tobias WEBER committed
506
507
508
509
510
511
512
513
		t_vec m_baseA;
		t_vec m_baseB;
		t_vec m_baseC;

		std::vector<std::string> m_vecAtoms;
		std::vector<unsigned int> m_vecAtomNums;

		std::vector<MolFrame<t_real, t_vec>> m_frames;
514
515

		boost::signals2::signal<bool (t_real)> m_sigLoadProgress;
Tobias WEBER's avatar
Tobias WEBER committed
516
		boost::signals2::signal<bool (t_real)> m_sigSaveProgress;
Tobias WEBER's avatar
Tobias WEBER committed
517
518
519
520
};


#endif