My Project
Loading...
Searching...
No Matches
write_data.h
1#include <iostream>
2#include <filesystem>
3#include <fstream>
4#include <vector>
5#include <string>
6#include <iomanip>
7
8
9static inline void WriteBox_LAMMPS(Atoms* System, Components SystemComponents, ForceField FF, Boxsize Box, std::ofstream& textrestartFile, std::vector<std::string> AtomNames)
10{
11 size_t NumberOfAtoms = 0;
12 for(size_t i = 0; i < SystemComponents.NComponents.x; i++)
13 {
14 NumberOfAtoms += SystemComponents.NumberOfMolecule_for_Component[i] * SystemComponents.Moleculesize[i];
15 printf("Printing: Component: %zu [ %s ], NumMol: %zu, Molsize: %zu\n", i, SystemComponents.MoleculeName[i].c_str(), SystemComponents.NumberOfMolecule_for_Component[i], SystemComponents.Moleculesize[i]);
16 }
17 textrestartFile << "# LAMMPS data file written by WriteLammpsdata function in RASPA(written by Zhao Li)" << '\n';
18 textrestartFile << NumberOfAtoms << " atoms" << '\n';
19 textrestartFile << 0 << " bonds" << '\n';
20 textrestartFile << 0 << " angles" << '\n';
21 textrestartFile << 0 << " dihedrals" << '\n';
22 textrestartFile << 0 << " impropers" << '\n';
23 textrestartFile << 0 << " " << Box.Cell[0] << " xlo xhi" << '\n';
24 textrestartFile << 0 << " " << Box.Cell[4] << " ylo yhi" << '\n';
25 textrestartFile << 0 << " " << Box.Cell[8] << " zlo zhi" << '\n';
26 textrestartFile << Box.Cell[3] << " " << Box.Cell[6] << " " << Box.Cell[7] << " xy xz yz" << '\n' << '\n';
27 textrestartFile << FF.size << " atom types" << '\n';
28 textrestartFile << 0 << " bond types" << '\n';
29 textrestartFile << 0 << " angle types" << '\n';
30 textrestartFile << 0 << " dihedral types" << '\n';
31 textrestartFile << 0 << " improper types" << '\n' << '\n';
32 //textrestartFile << std::puts("%zu bonds\n", NumberOfAtoms);
33 textrestartFile << "Masses" << '\n' << '\n';
34 for(size_t i = 0; i < FF.size; i++)
35 {
36 double mass = SystemComponents.PseudoAtoms.mass[i];
37 textrestartFile << i+1 << " " << mass << " # " << AtomNames[i] << '\n';
38 }
39 textrestartFile << '\n' << "Pair Coeffs" << '\n' << '\n';
40 for(size_t i = 0; i < FF.size; i++)
41 {
42 const size_t row = i*FF.size+i;
43 textrestartFile << i+1 << " " << FF.epsilon[row]*0.00239006 << " " << FF.sigma[row] << " # " << AtomNames[i] << '\n';
44 }
45}
46
47static inline void WriteAtoms_LAMMPS(Atoms* System, Components SystemComponents, Boxsize Box, std::ofstream& textrestartFile, std::vector<std::string> AtomNames)
48{
49 textrestartFile << '\n' << "Atoms" << '\n' << '\n';
50 size_t Atomcount=0; size_t Molcount=0;
51 for(size_t i = 0; i < SystemComponents.NComponents.x; i++)
52 {
53 Atoms Data = System[i];
54 printf("Component %zu, Molsize: %zu\n", i, Data.Molsize);
55 double3 Wrap_shift;
56 for(size_t j = 0; j < Data.size; j++)
57 {
58 //Wrap To Box//
59 double3 pos = Data.pos[j];
60 if(j % Data.Molsize == 0)
61 {
62 WrapInBox(pos, Box.Cell, Box.InverseCell, Box.Cubic);
63 Wrap_shift = pos - Data.pos[j];
64 }
65 else
66 {
67 pos += Wrap_shift;
68 }
69 textrestartFile << Atomcount+1 << " " << Molcount+Data.MolID[j]+1 << " " << Data.Type[j]+1 << " " << Data.charge[j] << " " << pos.x << " " << pos.y << " " << pos.z << " # " << SystemComponents.MoleculeName[i] << " " << AtomNames[Data.Type[j]] << '\n';
70 Atomcount++;
71 }
72 Molcount+=SystemComponents.NumberOfMolecule_for_Component[i];
73 }
74}
75
76static inline void WriteAtoms_Restart(Atoms* System, Components SystemComponents, std::ofstream& textrestartFile, std::vector<std::string> AtomNames)
77{
78 textrestartFile << "Reactions: 0" << "\n";
79 //size_t Atomcount=0;
80 size_t prevMol = 0;
81 for(size_t i = SystemComponents.NComponents.y; i < SystemComponents.NComponents.x; i++)
82 {
83 Atoms Data = System[i];
84 textrestartFile << '\n' << "Component: " << i - SystemComponents.NComponents.y << " Adsorbate " << SystemComponents.NumberOfMolecule_for_Component[i] << " molecules of " << SystemComponents.MoleculeName[i] << '\n';
85 textrestartFile << "------------------------------------------------------------------------" << "\n";
86 size_t molsize = SystemComponents.Moleculesize[i];
87 //First write positions//
88 for(size_t j = 0; j < Data.size; j++)
89 textrestartFile << "Adsorbate-atom-position:" << " " << prevMol + Data.MolID[j] << " " << j - Data.MolID[j]*molsize << " " << std::setprecision (15) << Data.pos[j].x << " " << std::setprecision (15) << Data.pos[j].y << " " << std::setprecision (15) << Data.pos[j].z << '\n';
90 //Then write velocities (Zhao's note: Not implemented yet)//
91 for(size_t j = 0; j < Data.size; j++)
92 textrestartFile << "Adsorbate-atom-velocity:" << " " << prevMol + Data.MolID[j] << " " << j - Data.MolID[j]*molsize << " " << 0.0 << " " << 0.0 << " " << 0.0 << '\n';
93 //Then write force (Zhao's note: Not implemented yet)//
94 for(size_t j = 0; j < Data.size; j++)
95 textrestartFile << "Adsorbate-atom-force:" << " " << prevMol + Data.MolID[j] << " " << j - Data.MolID[j]*molsize << " " << 0.0 << " " << 0.0 << " " << 0.0 << '\n';
96 //Then write charge//
97 for(size_t j = 0; j < Data.size; j++)
98 textrestartFile << "Adsorbate-atom-charge:" << " " << prevMol + Data.MolID[j] << " " << j - Data.MolID[j]*molsize << " " << Data.charge[j] << '\n';
99 //Then write scale//
100 for(size_t j = 0; j < Data.size; j++)
101 textrestartFile << "Adsorbate-atom-scaling:" << " " << prevMol + Data.MolID[j] << " " << j - Data.MolID[j]*molsize << " " << Data.scale[j] << '\n';
102 //Finally write fixed (Zhao's note: Not implemented yet)//
103 for(size_t j = 0; j < Data.size; j++)
104 textrestartFile << "Adsorbate-atom-fixed:" << " " << prevMol + Data.MolID[j] << " " << j - Data.MolID[j]*molsize << " " << "0 0 0" << '\n';
105 prevMol += SystemComponents.NumberOfMolecule_for_Component[i];
106 }
107}
108
109static inline void WriteComponent_Restart(Atoms* System, Components SystemComponents, std::ofstream& textrestartFile, Boxsize& Box)
110{
111 textrestartFile << "Components: " << SystemComponents.Total_Components - SystemComponents.NumberOfFrameworks << " (Adsorbates " << SystemComponents.TotalNumberOfMolecules-SystemComponents.NumberOfFrameworks << ", Cations 0)" << "\n";
112 textrestartFile << "========================================================================\n";
113 for(size_t i = SystemComponents.NumberOfFrameworks; i < SystemComponents.Total_Components; i++)
114 {
115 textrestartFile << "Components 0 (" << SystemComponents.MoleculeName[i] << ") " << "\n";
116 int fracID = -1;
117 if(SystemComponents.hasfractionalMolecule[i])
118 {
119 fracID = SystemComponents.Lambda[i].FractionalMoleculeID;
120 textrestartFile << "Fractional-molecule-id component " << i-SystemComponents.NumberOfFrameworks << ": " << fracID << "\n";
121 textrestartFile << "Lambda-factors component " << i-SystemComponents.NumberOfFrameworks << ": " << SystemComponents.Lambda[i].WangLandauScalingFactor << "\n";
122 textrestartFile << "Number-of-biasing-factors component " << i-SystemComponents.NumberOfFrameworks << ": " << SystemComponents.Lambda[i].binsize << "\n";
123 textrestartFile << "Biasing-factors component " << i-SystemComponents.NumberOfFrameworks << ": ";
124 for(size_t j = 0; j < SystemComponents.Lambda[i].binsize; j++)
125 textrestartFile << SystemComponents.Lambda[i].biasFactor[j] << " ";
126 textrestartFile << "\n";
127 textrestartFile << "Maximum-CF-Lambda-change component " << i-SystemComponents.NumberOfFrameworks << ": 0.50000\n"; //Zhao's note: continuous lambda not implemented
128 textrestartFile << "Maximum-CBCF-Lambda-change component " << i-SystemComponents.NumberOfFrameworks << ": 0.50000\n"; //Zhao's note: continuous lambda not implemented
129 }
130 textrestartFile << "\n";
131 textrestartFile << "Maximum-translation-change component " << i-SystemComponents.NumberOfFrameworks << ": " << SystemComponents.MaxTranslation[i].x << " " << SystemComponents.MaxTranslation[i].y << " " << SystemComponents.MaxTranslation[i].z << "\n";
132 textrestartFile << "Maximum-translation-in-plane-change component " << i-SystemComponents.NumberOfFrameworks << ": " << "0.000000,0.000000,0.000000" << "\n";
133 textrestartFile << "Maximum-rotation-change component " << i-SystemComponents.NumberOfFrameworks << ": " << SystemComponents.MaxRotation[i].x << " " << SystemComponents.MaxRotation[i].y << " " << SystemComponents.MaxRotation[i].z << "\n";
134 textrestartFile << "\n";
135 }
136}
137
138static inline void WriteCellInfo_Restart(Atoms* System, Components SystemComponents, std::ofstream& textrestartFile, Boxsize& Box)
139{
140 //size_t Atomcount=0;
141 textrestartFile << "Cell info:\n";
142 textrestartFile << "========================================================================\n";
143 textrestartFile << "number-of-unit-cells: 1 1 1\n"; //Zhao's note: not allowing for multiple unit cells for now//
144 //First write Unit Box sizes//
145 textrestartFile << "unit-cell-vector-a:" << " " << Box.Cell[0] << " " << Box.Cell[1] << " " << Box.Cell[2] << '\n';
146 textrestartFile << "unit-cell-vector-b:" << " " << Box.Cell[3] << " " << Box.Cell[4] << " " << Box.Cell[5] << '\n';
147 textrestartFile << "unit-cell-vector-c:" << " " << Box.Cell[6] << " " << Box.Cell[7] << " " << Box.Cell[8] << '\n';
148 textrestartFile << "\n";
149 //Then write Total Box sizes//
150 textrestartFile << "cell-vector-a:" << " " << Box.Cell[0] << " " << Box.Cell[1] << " " << Box.Cell[2] << '\n';
151 textrestartFile << "cell-vector-b:" << " " << Box.Cell[3] << " " << Box.Cell[4] << " " << Box.Cell[5] << '\n';
152 textrestartFile << "cell-vector-c:" << " " << Box.Cell[6] << " " << Box.Cell[7] << " " << Box.Cell[8] << '\n';
153 textrestartFile << "\n";
154 //Cell lengths//
155 textrestartFile << "cell-lengths:" << " " << Box.Cell[0] << " " << Box.Cell[4] << " " << Box.Cell[8] << " " << "\n";
156 textrestartFile << "cell-angles:" << " " << 90.00 << " " << 90.00 << " " << 90.00 << " " << "\n";
157 textrestartFile << "\n\n";
158 //Maximum changes for MC-moves (Zhao's note: Not implemented)//
159 textrestartFile << "Maximum changes for MC-moves:" << "\n";
160 textrestartFile << "========================================================================" << "\n";
161 textrestartFile << "Maximum-volume-change: 0.006250" << "\n";
162 textrestartFile << "Maximum-Gibbs-volume-change: 0.025000" << "\n";
163 textrestartFile << "Maximum-box-shape-change: 0.100000 0.100000 0.100000, 0.100000 0.100000 0.100000, 0.100000 0.100000 0.100000" << "\n";
164 textrestartFile << "\n\n";
165 //Acceptance targets for MC-moves (Zhao's note: Not implemented)//
166 textrestartFile << "Acceptance targets for MC-moves:" << "\n";
167 textrestartFile << "========================================================================" << "\n";
168 textrestartFile << "Target-volume-change: 0.500000" << "\n";
169 textrestartFile << "Target-box-shape-change: 0.500000" << "\n";
170 textrestartFile << "Target-Gibbs-volume-change: 0.500000" << "\n";
171 textrestartFile << "\n\n";
172 //Write Component Data//
173 WriteComponent_Restart(System, SystemComponents, textrestartFile, Box);
174}
175
176static inline void create_movie_file(Atoms* System, Components& SystemComponents, Boxsize& HostBox, std::vector<std::string> AtomNames, size_t SystemIndex)
177{
178 //std::ofstream textrestartFile{};
179 std::string dirname="Movies/System_" + std::to_string(SystemIndex) + "/";
180 std::string fname = dirname + "/" + "result_" + std::to_string(SystemComponents.CURRENTCYCLE) + ".data";
181 std::filesystem::path cwd = std::filesystem::current_path();
182
183 std::filesystem::path directoryName = cwd /dirname;
184 std::filesystem::path fileName = cwd /fname;
185 std::filesystem::create_directories(directoryName);
186
187 std::ofstream textrestartFile = std::ofstream(fileName, std::ios::out);
188
189 WriteBox_LAMMPS(System, SystemComponents, SystemComponents.FF, HostBox, textrestartFile, AtomNames);
190 WriteAtoms_LAMMPS(System, SystemComponents, HostBox, textrestartFile, AtomNames);
191 textrestartFile.close();
192}
193
194static inline void copyFile(const std::string& sourcePath, const std::string& destinationPath)
195{
196 std::ifstream source(sourcePath, std::ios::binary);
197 std::ofstream destination(destinationPath, std::ios::binary);
198
199 if (!source.is_open() || !destination.is_open())
200 {
201 std::cerr << "Error opening file." << std::endl;
202 return; // false;
203 }
204
205 destination << source.rdbuf();
206
207 source.close();
208 destination.close();
209 //return true;
210}
211
212static inline void create_Restart_file(size_t Cycle, Atoms* System, Components SystemComponents, ForceField FF, Boxsize Box, std::vector<std::string> AtomNames, size_t SystemIndex)
213{
214 std::ofstream textrestartFile{};
215 std::string dirname="Restart/System_" + std::to_string(SystemIndex) + "/";
216 std::string fname = dirname + "/" + "restartfile";
217 std::string backup_fname = dirname + "/" + "previous_restartfile";
218 std::filesystem::path cwd = std::filesystem::current_path();
219
220 //Keep a previous copy of the current restart file//
221 copyFile(fname, backup_fname);
222
223 std::filesystem::path directoryName = cwd /dirname;
224 std::filesystem::path fileName = cwd /fname;
225 std::filesystem::create_directories(directoryName);
226 textrestartFile = std::ofstream(fileName, std::ios::out);
227 WriteCellInfo_Restart(System, SystemComponents, textrestartFile, Box);
228 WriteAtoms_Restart(System, SystemComponents, textrestartFile, AtomNames);
229 textrestartFile.close();
230}
231
232static inline void WriteAllData(Atoms* System, Components SystemComponents, std::ofstream& textrestartFile, std::vector<std::string> AtomNames, size_t i)
233{
234 //size_t Atomcount=0;
235 Atoms Data = System[i];
236 size_t molsize = SystemComponents.Moleculesize[i];
237 //First write positions//
238 //printf("Writing AllData for Component %zu, There are %zu atoms, molsize: %zu\n", i, Data.size, molsize);
239 textrestartFile << "x y z charge scale scaleCoul Type" << '\n';
240 for(size_t j = 0; j < Data.size; j++)
241 textrestartFile << " " << Data.MolID[j] << " " << j - Data.MolID[j]*molsize << " " << Data.pos[j].x << " " << Data.pos[j].y << " " << Data.pos[j].z << " " << Data.charge[j] << " " << Data.scale[j] << " " << Data.scaleCoul[j] << " " << Data.Type[j] << '\n';
242}
243
244static inline void Write_All_Adsorbate_data(size_t Cycle, Atoms* System, Components SystemComponents, ForceField FF, Boxsize Box, std::vector<std::string> AtomNames, size_t SystemIndex)
245{
246 //std::ofstream textrestartFile{};
247 std::filesystem::path cwd = std::filesystem::current_path();
248
249 for(size_t i = 0; i < SystemComponents.NComponents.x; i++)
250 {
251 std::string dirname="AllData/System_" + std::to_string(SystemIndex) + "/";
252 std::string fname = dirname + "/" + "Component_" + std::to_string(i) + ".data";
253
254 std::filesystem::path directoryName = cwd /dirname;
255 std::filesystem::path fileName = cwd /fname;
256 std::filesystem::create_directories(directoryName);
257 std::ofstream textrestartFile = std::ofstream(fileName, std::ios::out);
258 WriteAllData(System, SystemComponents, textrestartFile, AtomNames, i);
259 textrestartFile.close();
260 }
261}
262
263//Zhao's note on 112623, moved the following to fxn_main.h//
264/*
265static inline void Write_Lambda(size_t Cycle, Components SystemComponents, size_t SystemIndex)
266{
267 std::ofstream textrestartFile{};
268 std::filesystem::path cwd = std::filesystem::current_path();
269
270 std::string dirname="Lambda/System_" + std::to_string(SystemIndex) + "/";
271 std::string fname = dirname + "/" + "Lambda_Histogram.data";
272
273 std::filesystem::path directoryName = cwd /dirname;
274 std::filesystem::path fileName = cwd /fname;
275 std::filesystem::create_directories(directoryName);
276 textrestartFile = std::ofstream(fileName, std::ios::out);
277 for(size_t i = 0; i < SystemComponents.Total_Components; i++)
278 {
279 if(SystemComponents.hasfractionalMolecule[i])
280 {
281 textrestartFile << "Component " << i << ": " << SystemComponents.MoleculeName[i] << '\n';
282 textrestartFile << "BIN SIZE : " << SystemComponents.Lambda[i].binsize << '\n';
283 textrestartFile << "BIN WIDTH: " << SystemComponents.Lambda[i].delta << '\n';
284 textrestartFile << "WL SCALING FACTOR: " << SystemComponents.Lambda[i].WangLandauScalingFactor << '\n';
285 textrestartFile << "FRACTIONAL MOLECULE ID: " << SystemComponents.Lambda[i].FractionalMoleculeID << '\n';
286 textrestartFile << "CURRENT BIN: " << SystemComponents.Lambda[i].currentBin << '\n';
287 textrestartFile << "BINS: ";
288 for(size_t j = 0; j < SystemComponents.Lambda[i].binsize; j++)
289 textrestartFile << j << " ";
290 textrestartFile << "\nHistogram: ";
291 for(size_t j = 0; j < SystemComponents.Lambda[i].binsize; j++)
292 textrestartFile << SystemComponents.Lambda[i].Histogram[j] << " ";
293 textrestartFile << "\nBIAS FACTOR: ";
294 for(size_t j = 0; j < SystemComponents.Lambda[i].binsize; j++)
295 textrestartFile << SystemComponents.Lambda[i].biasFactor[j] << " ";
296 }
297 }
298 textrestartFile.close();
299}
300
301static inline void Write_TMMC(size_t Cycle, Components SystemComponents, size_t SystemIndex)
302{
303 std::ofstream textTMMCFile{};
304 std::filesystem::path cwd = std::filesystem::current_path();
305
306 std::string dirname="TMMC/System_" + std::to_string(SystemIndex) + "/";
307 std::string fname = dirname + "/" + "TMMC_Statistics.data";
308
309 std::filesystem::path directoryName = cwd /dirname;
310 std::filesystem::path fileName = cwd /fname;
311 std::filesystem::create_directories(directoryName);
312 textTMMCFile = std::ofstream(fileName, std::ios::out);
313 for(size_t i = 0; i < SystemComponents.Total_Components; i++)
314 {
315 if(SystemComponents.Tmmc[i].DoTMMC)
316 {
317 textTMMCFile << "Component " << i << ": " << SystemComponents.MoleculeName[i] << " -> Updated " << SystemComponents.Tmmc[i].TMUpdateTimes << " times \n";
318 textTMMCFile << "Min Macrostate : " << SystemComponents.Tmmc[i].MinMacrostate << '\n';
319 textTMMCFile << "Max Macrostate : " << SystemComponents.Tmmc[i].MaxMacrostate << '\n';
320 textTMMCFile << "Wang-Landau Factor : " << SystemComponents.Tmmc[i].WLFactor << '\n';
321 textTMMCFile << "N NMol Bin CM[-1] CM[0] CM[1] WLBias ln_g TMBias lnpi Forward_lnpi Reverse_lnpi Histogram" << '\n';
322 for(size_t j = 0; j < SystemComponents.Tmmc[i].Histogram.size(); j++)
323 {
324 size_t N = j / SystemComponents.Tmmc[i].nbinPerMacrostate;
325 size_t bin = j % SystemComponents.Tmmc[i].nbinPerMacrostate;
326 textTMMCFile << j << " " << N << " " << bin << " " << SystemComponents.Tmmc[i].CMatrix[j].x << " " << SystemComponents.Tmmc[i].CMatrix[j].y << " " << SystemComponents.Tmmc[i].CMatrix[j].z << " " << SystemComponents.Tmmc[i].WLBias[j] << " " << SystemComponents.Tmmc[i].ln_g[j] << " " << SystemComponents.Tmmc[i].TMBias[j] << " " << SystemComponents.Tmmc[i].lnpi[j] << " " << SystemComponents.Tmmc[i].forward_lnpi[j] << " " << SystemComponents.Tmmc[i].reverse_lnpi[j] << " " << SystemComponents.Tmmc[i].Histogram[j] << '\n';
327 }
328 }
329 }
330 textrestartFile.close();
331}
332*/
Atoms
Definition data_struct.h:746
Boxsize
Definition data_struct.h:819
Components
Definition data_struct.h:843
ForceField
Definition data_struct.h:794
Generated by doxygen 1.9.8