My Project
Loading...
Searching...
No Matches
fxn_main.h
1void Setup_RandomNumber(RandomNumber& Random, size_t SIZE);
2void Copy_Atom_data_to_device(size_t NumberOfComponents, Atoms* device_System, Atoms* System);
3void Update_Components_for_framework(Components& SystemComponents);
4
5void Setup_Temporary_Atoms_Structure(Atoms& TempMol, Atoms* System);
6
7void Setup_Box_Temperature_Pressure(Units& Constants, Components& SystemComponents, Boxsize& device_Box);
8
9void Prepare_ForceField(ForceField& FF, ForceField& device_FF, PseudoAtomDefinitions PseudoAtom);
10
11void Prepare_Widom(WidomStruct& Widom, Boxsize Box, Simulations& Sims, Components SystemComponents, Atoms* System, Move_Statistics MoveStats);
12
13template<typename T>
14T* CUDA_allocate_array(size_t N, T InitVal)
15{
16 T* device_x;
17 cudaMalloc(&device_x, N * sizeof(T)); checkCUDAError("Error allocating Malloc");
18 T array[N];
19 for(size_t i = 0; i < N; i++) array[i] = InitVal;
20 cudaMemcpy(device_x, array, N * sizeof(T), cudaMemcpyHostToDevice);
21 //cudaMemset(device_x, (T) InitVal, N * sizeof(T));
22 return device_x;
23}
24
25template<typename T>
26T* CUDA_copy_allocate_array(T* x, size_t N)
27{
28 T* device_x;
29 cudaMalloc(&device_x, N * sizeof(T)); checkCUDAError("Error allocating Malloc");
30 cudaMemcpy(device_x, x, N * sizeof(T), cudaMemcpyHostToDevice); checkCUDAError("double Error Memcpy");
31 return device_x;
32}
33
34void Prepare_TempSystem_On_Host(Atoms& TempSystem)
35{
36 size_t Allocate_size = 1024;
37 TempSystem.pos =(double3*) malloc(Allocate_size*sizeof(double3));
38 TempSystem.scale = (double*) malloc(Allocate_size*sizeof(double));
39 TempSystem.charge = (double*) malloc(Allocate_size*sizeof(double));
40 TempSystem.scaleCoul = (double*) malloc(Allocate_size*sizeof(double));
41 TempSystem.Type = (size_t*) malloc(Allocate_size*sizeof(size_t));
42 TempSystem.MolID = (size_t*) malloc(Allocate_size*sizeof(size_t));
43 TempSystem.size = 0;
44 TempSystem.Allocate_size = Allocate_size;
45}
46
47inline void Setup_RandomNumber(RandomNumber& Random, size_t SIZE)
48{
49 Random.randomsize = SIZE; Random.offset = 0;
50 Random.host_random = (double3*) malloc(Random.randomsize * sizeof(double3));
51 for (size_t i = 0; i < Random.randomsize; i++)
52 {
53 Random.host_random[i].x = Get_Uniform_Random();
54 Random.host_random[i].y = Get_Uniform_Random();
55 Random.host_random[i].z = Get_Uniform_Random();
56 }
57 //Add some padding to match the sequence of the previous code, pad it up to 1 million numbers//
58 for (size_t i = Random.randomsize * 3; i < 1000000; i++) Get_Uniform_Random();
59
60 cudaMalloc(&Random.device_random, Random.randomsize * sizeof(double3));
61 cudaMemcpy(Random.device_random, Random.host_random, Random.randomsize * sizeof(double3), cudaMemcpyHostToDevice);
62}
63
64inline void Copy_Atom_data_to_device(size_t NumberOfComponents, Atoms* device_System, Atoms* System)
65{
66 size_t required_size = 0;
67 for(size_t i = 0; i < NumberOfComponents; i++)
68 {
69 if(i == 0){ required_size = System[i].size;} else { required_size = System[i].Allocate_size; }
70 device_System[i].pos = CUDA_copy_allocate_array(System[i].pos, required_size);
71 device_System[i].scale = CUDA_copy_allocate_array(System[i].scale, required_size);
72 device_System[i].charge = CUDA_copy_allocate_array(System[i].charge, required_size);
73 device_System[i].scaleCoul = CUDA_copy_allocate_array(System[i].scaleCoul, required_size);
74 device_System[i].Type = CUDA_copy_allocate_array(System[i].Type, required_size);
75 device_System[i].MolID = CUDA_copy_allocate_array(System[i].MolID, required_size);
76 device_System[i].size = System[i].size;
77 device_System[i].Molsize = System[i].Molsize;
78 device_System[i].Allocate_size = System[i].Allocate_size;
79 }
80}
81
82inline void Update_Components_for_framework(Components& SystemComponents)
83{
84 //Fill in Non-Important values for the framework components//
85 //Other values should be filled in CheckFrameworkCIF function//
86 for(size_t i = 0; i < SystemComponents.NComponents.y; i++)
87 {
88 SystemComponents.MolFraction.push_back(1.0);
89 SystemComponents.IdealRosenbluthWeight.push_back(1.0);
90 SystemComponents.FugacityCoeff.push_back(1.0);
91 SystemComponents.Tc.push_back(0.0); //Tc for framework is set to zero
92 SystemComponents.Pc.push_back(0.0); //Pc for framework is set to zero
93 SystemComponents.Accentric.push_back(0.0); //Accentric factor for framework is set to zero
94 //Zhao's note: for now, assume the framework is rigid//
95 SystemComponents.rigid.push_back(true);
96 SystemComponents.hasfractionalMolecule.push_back(false); //No fractional molecule for the framework//
97 SystemComponents.NumberOfCreateMolecules.push_back(0); //Create zero molecules for the framework//
98 LAMBDA lambda;
99 lambda.newBin = 0; lambda.delta = static_cast<double>(1.0/(lambda.binsize)); lambda.WangLandauScalingFactor = 0.0; //Zhao's note: in raspa3, delta is 1/(nbin - 1)
100 lambda.FractionalMoleculeID = 0;
101 SystemComponents.Lambda.push_back(lambda);
102
103 TMMC tmmc;
104 SystemComponents.Tmmc.push_back(tmmc); //Just use default values for tmmc for the framework, it will do nothing//
105 }
106 //Add PseudoAtoms from the Framework to the total PseudoAtoms array//
107 //SystemComponents.UpdatePseudoAtoms(INSERTION, 0);
108}
109
110inline void Setup_Temporary_Atoms_Structure(Atoms& TempMol, Atoms* System)
111{
112 //Set up MolArrays//
113 size_t Allocate_size_Temporary=1024; //Assign 1024 empty slots for the temporary structures//
114 //OLD//
115 TempMol.pos = CUDA_allocate_array<double3> (Allocate_size_Temporary, {0.0, 0.0, 0.0});
116 TempMol.scale = CUDA_allocate_array<double> (Allocate_size_Temporary, 0.0);
117 TempMol.charge = CUDA_allocate_array<double> (Allocate_size_Temporary, 0.0);
118 TempMol.scaleCoul = CUDA_allocate_array<double> (Allocate_size_Temporary, 0.0);
119 TempMol.Type = CUDA_allocate_array<size_t> (Allocate_size_Temporary, 0.0);
120 TempMol.MolID = CUDA_allocate_array<size_t> (Allocate_size_Temporary, 0.0);
121 TempMol.size = 0;
122 TempMol.Molsize = 0;
123 TempMol.Allocate_size = Allocate_size_Temporary;
124}
125
126inline void Setup_Box_Temperature_Pressure(Units& Constants, Components& SystemComponents, Boxsize& device_Box)
127{
128 SystemComponents.Beta = 1.0/(Constants.BoltzmannConstant/(Constants.MassUnit*pow(Constants.LengthUnit,2)/pow(Constants.TimeUnit,2))*SystemComponents.Temperature);
129 //Convert pressure from pascal
130 device_Box.Pressure/=(Constants.MassUnit/(Constants.LengthUnit*pow(Constants.TimeUnit,2)));
131 printf("------------------- SIMULATION BOX PARAMETERS -----------------\n");
132 printf("Pressure: %.5f\n", device_Box.Pressure);
133 printf("Box Volume: %.5f\n", device_Box.Volume);
134 printf("Box Beta: %.5f\n", SystemComponents.Beta);
135 printf("Box Temperature: %.5f\n", SystemComponents.Temperature);
136 printf("---------------------------------------------------------------\n");
137}
138
139inline void Prepare_ForceField(ForceField& FF, ForceField& device_FF, PseudoAtomDefinitions PseudoAtom)
140{
141 // COPY DATA TO DEVICE POINTER //
142 //device_FF.FFParams = CUDA_copy_allocate_array(FF.FFParams, 5);
143 device_FF.OverlapCriteria = FF.OverlapCriteria;
144 device_FF.CutOffVDW = FF.CutOffVDW;
145 device_FF.CutOffCoul = FF.CutOffCoul;
146 //device_FF.Prefactor = FF.Prefactor;
147 //device_FF.Alpha = FF.Alpha;
148
149 device_FF.epsilon = CUDA_copy_allocate_array(FF.epsilon, FF.size*FF.size);
150 device_FF.sigma = CUDA_copy_allocate_array(FF.sigma, FF.size*FF.size);
151 device_FF.z = CUDA_copy_allocate_array(FF.z, FF.size*FF.size);
152 device_FF.shift = CUDA_copy_allocate_array(FF.shift, FF.size*FF.size);
153 device_FF.FFType = CUDA_copy_allocate_array(FF.FFType, FF.size*FF.size);
154 device_FF.noCharges = FF.noCharges;
155 device_FF.size = FF.size;
156 device_FF.VDWRealBias = FF.VDWRealBias;
157 //Formulate Component statistics on the host
158 //ForceFieldParser(FF, PseudoAtom);
159 //PseudoAtomParser(FF, PseudoAtom);
160}
161
162inline void InitializeMaxTranslationRotation(Components& SystemComponents)
163{
164 for(size_t i = 0; i < SystemComponents.NComponents.x; i++)
165 {
166 double3 MaxTranslation = {1.0, 1.0, 1.0};
167 double3 MaxRotation = {30.0/(180/3.1415), 30.0/(180/3.1415), 30.0/(180/3.1415)};
168 SystemComponents.MaxTranslation.push_back(MaxTranslation);
169 SystemComponents.MaxRotation.push_back(MaxRotation);
170 SystemComponents.MaxSpecialRotation.push_back(MaxRotation);
171 }
172}
173
174inline void Prepare_Widom(WidomStruct& Widom, Boxsize Box, Simulations& Sims, Components& SystemComponents, Atoms* System)
175{
176 //Zhao's note: NumberWidomTrials is for first bead. NumberWidomTrialsOrientations is for the rest, here we consider single component, not mixture //
177
178 size_t MaxTrialsize = max(Widom.NumberWidomTrials, Widom.NumberWidomTrialsOrientations*(SystemComponents.Moleculesize[1]-1));
179
180 //Zhao's note: The previous way yields a size for blocksum that can be smaller than the number of kpoints
181 //This is a problem when you need to do parallel Ewald summation for the whole system//
182 //Might be good to add a flag or so//
183 //size_t MaxResultsize = MaxTrialsize*(System[0].Allocate_size+System[1].Allocate_size);
184 size_t MaxAllocatesize = max(System[0].Allocate_size, System[1].Allocate_size);
185 size_t MaxResultsize = MaxTrialsize * SystemComponents.Total_Components * MaxAllocatesize * 5; //For Volume move, it really needs a lot of blocks//
186
187
188 printf("----------------- MEMORY ALLOCAION STATUS -----------------\n");
189 //Compare Allocate sizes//
190 printf("System allocate_sizes are: %zu, %zu\n", System[0].Allocate_size, System[1].Allocate_size);
191 printf("Component allocate_sizes are: %zu, %zu\n", SystemComponents.Allocate_size[0], SystemComponents.Allocate_size[1]);
192
193 SystemComponents.flag = (bool*)malloc(MaxTrialsize * sizeof(bool));
194 cudaMallocHost(&Sims.device_flag, MaxTrialsize * sizeof(bool));
195
196 cudaMallocHost(&Sims.Blocksum, (MaxResultsize/DEFAULTTHREAD + 1)*sizeof(double));
197
198 cudaMallocManaged(&Sims.ExcludeList, 10 * sizeof(int2));
199 for(size_t i = 0; i < 10; i++) Sims.ExcludeList[i] = {-1, -1}; //Initialize with negative # so that nothing is ignored//
200 //cudaMalloc(&Sims.Blocksum, (MaxResultsize/DEFAULTTHREAD + 1)*sizeof(double));
201
202 printf("Allocated Blocksum size: %zu\n", (MaxResultsize/DEFAULTTHREAD + 1));
203
204 //cudaMalloc(&Sims.Blocksum, (MaxResultsize/DEFAULTTHREAD + 1)*sizeof(double));
205 Sims.Nblocks = MaxResultsize/DEFAULTTHREAD + 1;
206
207 printf("Allocated %zu doubles for Blocksums\n", MaxResultsize/DEFAULTTHREAD + 1);
208
209 for(size_t i = 0; i < SystemComponents.NComponents.x; i++)
210 {
211 double3 MaxTranslation = {Box.Cell[0]*0.1, Box.Cell[4]*0.1, Box.Cell[8]*0.1};
212 double3 MaxRotation = {30.0/(180/3.1415), 30.0/(180/3.1415), 30.0/(180/3.1415)};
213 SystemComponents.MaxTranslation[i] =MaxTranslation;
214 SystemComponents.MaxRotation[i] =MaxRotation;
215 SystemComponents.MaxSpecialRotation[i]=MaxRotation;
216 }
217 Sims.start_position = 0;
218 //Sims.Nblocks = 0;
219 Sims.TotalAtoms = 0;
220 Sims.AcceptedFlag = false;
221
222 Widom.WidomFirstBeadAllocatesize = MaxResultsize/DEFAULTTHREAD;
223 printf("------------------------------------------------------------\n");
224}
225
226inline void Allocate_Copy_Ewald_Vector(Boxsize& device_Box, Components SystemComponents)
227{
228 printf("****** Allocating Ewald WaveVectors (INITIAL STAGE ONLY) ******\n");
229 //Zhao's note: This only works if the box size is not changed, eik_xy might not be useful if box size is not changed//
230 size_t eikx_size = SystemComponents.eik_x.size() * 2;
231 size_t eiky_size = SystemComponents.eik_y.size() * 2; //added times 2 for box volume move//
232 size_t eikz_size = SystemComponents.eik_z.size() * 2;
233 printf("Allocated %zu %zu %zu space for eikxyz\n", eikx_size, eiky_size, eikz_size);
234 //size_t eikxy_size = SystemComponents.eik_xy.size();
235 size_t AdsorbateEiksize = SystemComponents.AdsorbateEik.size() * 2; //added times 2 for box volume move//
236 cudaMalloc(&device_Box.eik_x, eikx_size * sizeof(Complex));
237 cudaMalloc(&device_Box.eik_y, eiky_size * sizeof(Complex));
238 cudaMalloc(&device_Box.eik_z, eikz_size * sizeof(Complex));
239 //cudaMalloc(&device_Box.eik_xy, eikxy_size * sizeof(Complex));
240 cudaMalloc(&device_Box.AdsorbateEik, AdsorbateEiksize * sizeof(Complex));
241 cudaMalloc(&device_Box.tempEik, AdsorbateEiksize * sizeof(Complex));
242 cudaMalloc(&device_Box.FrameworkEik, AdsorbateEiksize * sizeof(Complex));
243
244 Complex AdsorbateEik[AdsorbateEiksize]; //Temporary Complex struct on the host//
245 Complex FrameworkEik[AdsorbateEiksize];
246 //for(size_t i = 0; i < SystemComponents.AdsorbateEik.size(); i++)
247 for(size_t i = 0; i < AdsorbateEiksize; i++)
248 {
249 if(i < SystemComponents.AdsorbateEik.size())
250 {
251 AdsorbateEik[i].real = SystemComponents.AdsorbateEik[i].real();
252 AdsorbateEik[i].imag = SystemComponents.AdsorbateEik[i].imag();
253
254 FrameworkEik[i].real = SystemComponents.FrameworkEik[i].real();
255 FrameworkEik[i].imag = SystemComponents.FrameworkEik[i].imag();
256 }
257 else
258 {
259 AdsorbateEik[i].real = 0.0; AdsorbateEik[i].imag = 0.0;
260 FrameworkEik[i].real = 0.0; FrameworkEik[i].imag = 0.0;
261 }
262 if(i < 10) printf("Wave Vector %zu is %.5f %.5f\n", i, AdsorbateEik[i].real, AdsorbateEik[i].imag);
263 }
264 cudaMemcpy(device_Box.AdsorbateEik, AdsorbateEik, AdsorbateEiksize * sizeof(Complex), cudaMemcpyHostToDevice); checkCUDAError("error copying Complex");
265 cudaMemcpy(device_Box.FrameworkEik, FrameworkEik, AdsorbateEiksize * sizeof(Complex), cudaMemcpyHostToDevice); checkCUDAError("error copying Complex");
266 printf("****** DONE Allocating Ewald WaveVectors (INITIAL STAGE ONLY) ******\n");
267}
268
269inline void Check_Simulation_Energy(Boxsize& Box, Atoms* System, ForceField FF, ForceField device_FF, Components& SystemComponents, int SIMULATIONSTAGE, size_t Numsim, Simulations& Sim)
270{
271 std::string STAGE;
272 switch(SIMULATIONSTAGE)
273 {
274 case INITIAL:
275 { STAGE = "INITIAL"; break;}
276 case CREATEMOL:
277 { STAGE = "CREATE_MOLECULE"; break;}
278 case FINAL:
279 { STAGE = "FINAL"; break;}
280 }
281 printf("======================== CALCULATING %s STAGE ENERGY ========================\n", STAGE.c_str());
282 MoveEnergy ENERGY;
283 Atoms device_System[SystemComponents.Total_Components];
284 cudaMemcpy(device_System, Sim.d_a, SystemComponents.Total_Components * sizeof(Atoms), cudaMemcpyDeviceToHost);
285 cudaMemcpy(Box.Cell, Sim.Box.Cell, 9 * sizeof(double), cudaMemcpyDeviceToHost);
286 cudaMemcpy(Box.InverseCell, Sim.Box.InverseCell, 9 * sizeof(double), cudaMemcpyDeviceToHost);
287 //Update every value that can be changed during a volume move//
288 Box.Volume = Sim.Box.Volume;
289 Box.ReciprocalCutOff = Sim.Box.ReciprocalCutOff;
290 Box.Cubic = Sim.Box.Cubic;
291 Box.kmax = Sim.Box.kmax;
292
293 MoveEnergy GPU_Energy;
294
295 double start = omp_get_wtime();
296 VDWReal_Total_CPU(Box, System, device_System, FF, SystemComponents, ENERGY);
297 double end = omp_get_wtime();
298 double CPUSerialTime = end - start;
299 start = omp_get_wtime();
300 double* xxx; xxx = (double*) malloc(sizeof(double)*2);
301 double* device_xxx; device_xxx = CUDA_copy_allocate_array(xxx, 2);
302 //Zhao's note: if the serial GPU energy test is too slow, comment it out//
303 //one_thread_GPU_test<<<1,1>>>(Sim.Box, Sim.d_a, device_FF, device_xxx);
304 cudaMemcpy(xxx, device_xxx, sizeof(double), cudaMemcpyDeviceToHost);
305 end = omp_get_wtime();
306 cudaDeviceSynchronize();
307
308 double SerialGPUTime = end - start;
309 //For total energy, divide the parallelization into several parts//
310 //For framework, every thread treats the interaction between one framework atom with an adsorbate molecule//
311 //For adsorbate/adsorbate, every thread treats one adsorbate molecule with an adsorbate molecule//
312 start = omp_get_wtime();
313 size_t Host_threads = 0;
314 size_t Guest_threads = 0;
315 size_t NFrameworkAtomsPerThread = 4;
316 size_t NAdsorbate = 0;
317 for(size_t i = 1; i < SystemComponents.Total_Components; i++) NAdsorbate += SystemComponents.NumberOfMolecule_for_Component[i];
318 Host_threads = SystemComponents.Moleculesize[0] / NFrameworkAtomsPerThread; //Per adsorbate molecule//
319 if(SystemComponents.Moleculesize[0] % NFrameworkAtomsPerThread != 0) Host_threads ++;
320 Host_threads *= NAdsorbate; //Total = Host_thread_per_molecule * number of Adsorbate molecule
321 Guest_threads = NAdsorbate * (NAdsorbate-1)/2;
322 if(Host_threads + Guest_threads > 0)
323 {
324 bool ConsiderHostHost = false;
325 bool UseOffset = false;
326 GPU_Energy += Total_VDW_Coulomb_Energy(Sim, device_FF, NAdsorbate, Host_threads, Guest_threads, NFrameworkAtomsPerThread, ConsiderHostHost, UseOffset);
327 }
328 end = omp_get_wtime();
329
330 //Do Parallel Total Ewald//
331 double CPUEwaldTime = 0.0;
332 double GPUEwaldTime = 0.0;
333
334 if(!device_FF.noCharges)
335 {
336 cudaDeviceSynchronize();
337 double EwStart = omp_get_wtime();
338 CPU_GPU_EwaldTotalEnergy(Box, Sim.Box, System, Sim.d_a, FF, device_FF, SystemComponents, ENERGY);
339 ENERGY.GGEwaldE -= SystemComponents.FrameworkEwald;
340
341 double EwEnd = omp_get_wtime();
342 //Zhao's note: if it is in the initial stage, calculate the intra and self exclusion energy for ewald summation//
343 if(SIMULATIONSTAGE == INITIAL) Calculate_Exclusion_Energy_Rigid(Box, System, FF, SystemComponents);
344 CPUEwaldTime = EwEnd - EwStart;
345
346 cudaDeviceSynchronize();
347 //Zhao's note: if doing initial energy, initialize and copy host Ewald to device//
348 if(SIMULATIONSTAGE == INITIAL) Allocate_Copy_Ewald_Vector(Sim.Box, SystemComponents);
349 Check_WaveVector_CPUGPU(Sim.Box, SystemComponents); //Check WaveVector on the CPU and GPU//
350 cudaDeviceSynchronize();
351 EwStart = omp_get_wtime();
352 bool UseOffset = false;
353 GPU_Energy += Ewald_TotalEnergy(Sim, SystemComponents, UseOffset);
354 GPU_Energy.GGEwaldE -= SystemComponents.FrameworkEwald;
355 cudaDeviceSynchronize();
356 EwEnd = omp_get_wtime();
357 GPUEwaldTime = EwEnd - EwStart;
358 }
359 printf("Ewald Summation (total energy) on the CPU took %.5f secs\n", CPUEwaldTime);
360 printf("Ewald Summation (total energy) on the GPU took %.5f secs\n", GPUEwaldTime);
361 //Calculate Tail Correction Energy//
362 ENERGY.TailE = TotalTailCorrection(SystemComponents, FF.size, Sim.Box.Volume);
363 if(SystemComponents.UseDNNforHostGuest) ENERGY.DNN_E = DNN_Prediction_Total(SystemComponents, Sim);
364 if(SystemComponents.UseDNNforHostGuest) double Correction = ENERGY.DNN_Correction();
365
366 if(SIMULATIONSTAGE == INITIAL) SystemComponents.Initial_Energy = ENERGY;
367 else if(SIMULATIONSTAGE == CREATEMOL)
368 {
369 SystemComponents.CreateMol_Energy = ENERGY;
370 }
371 else
372 {
373 SystemComponents.Final_Energy = ENERGY;
374 }
375 printf("====================== DONE CALCULATING %s STAGE ENERGY ======================\n", STAGE.c_str());
376}
377
378inline void Copy_AtomData_from_Device(Atoms* System, Atoms* Host_System, Atoms* d_a, Components& SystemComponents)
379{
380 cudaMemcpy(System, d_a, SystemComponents.Total_Components * sizeof(Atoms), cudaMemcpyDeviceToHost);
381 //printval<<<1,1>>>(System[0]);
382 //printvald_a<<<1,1>>>(d_a);
383 for(size_t ijk=0; ijk < SystemComponents.Total_Components; ijk++)
384 {
385 // if the host allocate_size is different from the device, allocate more space on the host
386 size_t current_allocated_size = System[ijk].Allocate_size;
387 if(current_allocated_size != Host_System[ijk].Allocate_size) //Need to update host
388 {
389 Host_System[ijk].pos = (double3*) malloc(System[ijk].Allocate_size*sizeof(double3));
390 Host_System[ijk].scale = (double*) malloc(System[ijk].Allocate_size*sizeof(double));
391 Host_System[ijk].charge = (double*) malloc(System[ijk].Allocate_size*sizeof(double));
392 Host_System[ijk].scaleCoul = (double*) malloc(System[ijk].Allocate_size*sizeof(double));
393 Host_System[ijk].Type = (size_t*) malloc(System[ijk].Allocate_size*sizeof(size_t));
394 Host_System[ijk].MolID = (size_t*) malloc(System[ijk].Allocate_size*sizeof(size_t));
395 Host_System[ijk].Allocate_size = System[ijk].Allocate_size;
396 }
397 Host_System[ijk].size = System[ijk].size; //Zhao's note: no matter what, the size (not allocated size) needs to be updated
398
399 cudaMemcpy(Host_System[ijk].pos, System[ijk].pos, sizeof(double3)*System[ijk].Allocate_size, cudaMemcpyDeviceToHost);
400 cudaMemcpy(Host_System[ijk].scale, System[ijk].scale, sizeof(double)*System[ijk].Allocate_size, cudaMemcpyDeviceToHost);
401 cudaMemcpy(Host_System[ijk].charge, System[ijk].charge, sizeof(double)*System[ijk].Allocate_size, cudaMemcpyDeviceToHost);
402 cudaMemcpy(Host_System[ijk].scaleCoul, System[ijk].scaleCoul, sizeof(double)*System[ijk].Allocate_size, cudaMemcpyDeviceToHost);
403 cudaMemcpy(Host_System[ijk].Type, System[ijk].Type, sizeof(size_t)*System[ijk].Allocate_size, cudaMemcpyDeviceToHost);
404 cudaMemcpy(Host_System[ijk].MolID, System[ijk].MolID, sizeof(size_t)*System[ijk].Allocate_size, cudaMemcpyDeviceToHost);
405 Host_System[ijk].size = System[ijk].size;
406 }
407}
408
409inline void PRINT_ENERGY_AT_STAGE(Components& SystemComponents, int stage, Units& Constants)
410{
411 std::string stage_name;
412 MoveEnergy E;
413 switch(stage)
414 {
415 case INITIAL: {stage_name = "INITIAL STAGE"; E = SystemComponents.Initial_Energy; break;}
416 case CREATEMOL: {stage_name = "CREATE MOLECULE STAGE"; E = SystemComponents.CreateMol_Energy; break;}
417 case FINAL: {stage_name = "FINAL STAGE"; E = SystemComponents.Final_Energy; break;}
418 case CREATEMOL_DELTA: {stage_name = "RUNNING DELTA_E (CREATE MOLECULE - INITIAL)"; E = SystemComponents.CreateMoldeltaE; break;}
419 case DELTA: {stage_name = "RUNNING DELTA_E (FINAL - CREATE MOLECULE)"; E = SystemComponents.deltaE; break;}
420 case CREATEMOL_DELTA_CHECK: {stage_name = "CHECK DELTA_E (CREATE MOLECULE - INITIAL)"; E = SystemComponents.CreateMol_Energy - SystemComponents.Initial_Energy; break;}
421 case DELTA_CHECK: {stage_name = "CHECK DELTA_E (FINAL - CREATE MOLECULE)"; E = SystemComponents.Final_Energy - SystemComponents.CreateMol_Energy; break;}
422 case DRIFT: {stage_name = "ENERGY DRIFT"; E = SystemComponents.CreateMol_Energy + SystemComponents.deltaE - SystemComponents.Final_Energy; break;}
423 case AVERAGE: {stage_name = "PRODUCTION PHASE AVERAGE ENERGY"; E = SystemComponents.AverageEnergy;break;}
424 case AVERAGE_ERR: {stage_name = "PRODUCTION PHASE AVERAGE ENERGY ERRORBAR"; E = SystemComponents.AverageEnergy_Errorbar; break;}
425 }
426 printf(" *** %s *** \n", stage_name.c_str());
427 printf("========================================================================\n");
428 printf("VDW [Host-Host]: %.5f (%.5f [K])\n", E.HHVDW, E.HHVDW * Constants.energy_to_kelvin);
429 printf("VDW [Host-Guest]: %.5f (%.5f [K])\n", E.HGVDW, E.HGVDW * Constants.energy_to_kelvin);
430 printf("VDW [Guest-Guest]: %.5f (%.5f [K])\n", E.GGVDW, E.GGVDW * Constants.energy_to_kelvin);
431 printf("Real Coulomb [Host-Host]: %.5f (%.5f [K])\n", E.HHReal, E.HHReal * Constants.energy_to_kelvin);
432 printf("Real Coulomb [Host-Guest]: %.5f (%.5f [K])\n", E.HGReal, E.HGReal * Constants.energy_to_kelvin);
433 printf("Real Coulomb [Guest-Guest]: %.5f (%.5f [K])\n", E.GGReal, E.GGReal * Constants.energy_to_kelvin);
434 printf("Ewald [Host-Host]: %.5f (%.5f [K])\n", E.HHEwaldE, E.HHEwaldE * Constants.energy_to_kelvin);
435 printf("Ewald [Host-Guest]: %.5f (%.5f [K])\n", E.HGEwaldE, E.HGEwaldE * Constants.energy_to_kelvin);
436 printf("Ewald [Guest-Guest]: %.5f (%.5f [K])\n", E.GGEwaldE, E.GGEwaldE * Constants.energy_to_kelvin);
437 printf("DNN Energy: %.5f (%.5f [K])\n", E.DNN_E, E.DNN_E * Constants.energy_to_kelvin);
438 if(SystemComponents.UseDNNforHostGuest)
439 {
440 printf(" --> Stored Classical Host-Guest Interactions: \n");
441 printf(" VDW: %.5f (%.5f [K])\n", E.storedHGVDW, E.storedHGVDW * Constants.energy_to_kelvin);
442 printf(" Real Coulomb: %.5f (%.5f [K])\n", E.storedHGReal, E.storedHGReal * Constants.energy_to_kelvin);
443 printf(" Ewald: %.5f (%.5f [K])\n", E.storedHGEwaldE, E.storedHGEwaldE * Constants.energy_to_kelvin);
444 printf(" Total: %.5f (%.5f [K])\n", E.storedHGVDW + E.storedHGReal + E.storedHGEwaldE, (E.storedHGVDW + E.storedHGReal + E.storedHGEwaldE) * Constants.energy_to_kelvin);
445 printf(" --> DNN - Classical: %.5f (%.5f [K])\n", E.DNN_E - (E.storedHGVDW + E.storedHGReal + E.storedHGEwaldE), (E.DNN_E - (E.storedHGVDW + E.storedHGReal + E.storedHGEwaldE)) * Constants.energy_to_kelvin);
446 }
447 printf("Tail Correction Energy: %.5f (%.5f [K])\n", E.TailE, E.TailE * Constants.energy_to_kelvin);
448 printf("Total Energy: %.5f (%.5f [K])\n", E.total(), E.total() * Constants.energy_to_kelvin);
449 printf("========================================================================\n");
450}
451inline void ENERGY_SUMMARY(std::vector<Components>& SystemComponents, Units& Constants)
452{
453 size_t NumberOfSimulations = SystemComponents.size();
454 for(size_t i = 0; i < NumberOfSimulations; i++)
455 {
456 printf("======================== ENERGY SUMMARY (Simulation %zu) =========================\n", i);
457 PRINT_ENERGY_AT_STAGE(SystemComponents[i], INITIAL, Constants);
458 PRINT_ENERGY_AT_STAGE(SystemComponents[i], CREATEMOL, Constants);
459 PRINT_ENERGY_AT_STAGE(SystemComponents[i], CREATEMOL_DELTA, Constants);
460 PRINT_ENERGY_AT_STAGE(SystemComponents[i], CREATEMOL_DELTA_CHECK, Constants);
461 PRINT_ENERGY_AT_STAGE(SystemComponents[i], FINAL, Constants);
462 PRINT_ENERGY_AT_STAGE(SystemComponents[i], DELTA, Constants);
463 PRINT_ENERGY_AT_STAGE(SystemComponents[i], DELTA_CHECK, Constants);
464 PRINT_ENERGY_AT_STAGE(SystemComponents[i], DRIFT, Constants);
465 printf("================================================================================\n");
466 printf("======================== PRODUCTION PHASE AVERAGE ENERGIES (Simulation %zu) =========================\n", i);
467 PRINT_ENERGY_AT_STAGE(SystemComponents[i], AVERAGE, Constants);
468 PRINT_ENERGY_AT_STAGE(SystemComponents[i], AVERAGE_ERR, Constants);
469 printf("================================================================================\n");
470
471 printf("DNN Rejection Summary:\nTranslation+Rotation: %zu\nReinsertion: %zu\nInsertion: %zu\nDeletion: %zu\nSingleSwap: %zu\n", SystemComponents[i].TranslationRotationDNNReject, SystemComponents[i].ReinsertionDNNReject, SystemComponents[i].InsertionDNNReject, SystemComponents[i].DeletionDNNReject, SystemComponents[i].SingleSwapDNNReject);
472 printf("DNN Drift Summary:\nTranslation+Rotation: %.5f\nReinsertion: %.5f\nInsertion: %.5f\nDeletion: %.5f\nSingleSwap: %.5f\n", SystemComponents[i].SingleMoveDNNDrift, SystemComponents[i].ReinsertionDNNDrift, SystemComponents[i].InsertionDNNDrift, SystemComponents[i].DeletionDNNDrift, SystemComponents[i].SingleSwapDNNDrift);
473 }
474}
475
476static inline void Write_Lambda(size_t Cycle, Components SystemComponents, size_t SystemIndex)
477{
478 std::ofstream textrestartFile{};
479 std::filesystem::path cwd = std::filesystem::current_path();
480
481 std::string dirname="Lambda/System_" + std::to_string(SystemIndex) + "/";
482 std::string fname = dirname + "/" + "Lambda_Histogram.data";
483
484 std::filesystem::path directoryName = cwd /dirname;
485 std::filesystem::path fileName = cwd /fname;
486 std::filesystem::create_directories(directoryName);
487 textrestartFile = std::ofstream(fileName, std::ios::out);
488 for(size_t i = 0; i < SystemComponents.Total_Components; i++)
489 {
490 if(SystemComponents.hasfractionalMolecule[i])
491 {
492 textrestartFile << "Component " << i << ": " << SystemComponents.MoleculeName[i] << '\n';
493 textrestartFile << "BIN SIZE : " << SystemComponents.Lambda[i].binsize << '\n';
494 textrestartFile << "BIN WIDTH: " << SystemComponents.Lambda[i].delta << '\n';
495 textrestartFile << "WL SCALING FACTOR: " << SystemComponents.Lambda[i].WangLandauScalingFactor << '\n';
496 textrestartFile << "FRACTIONAL MOLECULE ID: " << SystemComponents.Lambda[i].FractionalMoleculeID << '\n';
497 textrestartFile << "CURRENT BIN: " << SystemComponents.Lambda[i].currentBin << '\n';
498 textrestartFile << "BINS: ";
499 for(size_t j = 0; j < SystemComponents.Lambda[i].binsize; j++)
500 textrestartFile << j << " ";
501 textrestartFile << "\nHistogram: ";
502 for(size_t j = 0; j < SystemComponents.Lambda[i].binsize; j++)
503 textrestartFile << SystemComponents.Lambda[i].Histogram[j] << " ";
504 textrestartFile << "\nBIAS FACTOR: ";
505 for(size_t j = 0; j < SystemComponents.Lambda[i].binsize; j++)
506 textrestartFile << SystemComponents.Lambda[i].biasFactor[j] << " ";
507 }
508 }
509 textrestartFile.close();
510}
511
512static inline void Write_TMMC(size_t Cycle, Components SystemComponents, size_t SystemIndex)
513{
514 std::ofstream textTMMCFile{};
515 std::filesystem::path cwd = std::filesystem::current_path();
516
517 std::string dirname="TMMC/System_" + std::to_string(SystemIndex) + "/";
518 std::string fname = dirname + "/" + "TMMC_Statistics.data";
519
520 std::filesystem::path directoryName = cwd /dirname;
521 std::filesystem::path fileName = cwd /fname;
522 std::filesystem::create_directories(directoryName);
523 textTMMCFile = std::ofstream(fileName, std::ios::out);
524 for(size_t i = 0; i < SystemComponents.Total_Components; i++)
525 {
526 if(SystemComponents.Tmmc[i].DoTMMC)
527 {
528 textTMMCFile << "Component " << i << ": " << SystemComponents.MoleculeName[i] << " -> Updated " << SystemComponents.Tmmc[i].TMUpdateTimes << " times \n";
529 textTMMCFile << "Min Macrostate : " << SystemComponents.Tmmc[i].MinMacrostate << '\n';
530 textTMMCFile << "Max Macrostate : " << SystemComponents.Tmmc[i].MaxMacrostate << '\n';
531 textTMMCFile << "Wang-Landau Factor : " << SystemComponents.Tmmc[i].WLFactor << '\n';
532 textTMMCFile << "N NMol Bin CM[-1] CM[0] CM[1] WLBias ln_g TMBias lnpi Forward_lnpi Reverse_lnpi Histogram" << '\n';
533 for(size_t j = 0; j < SystemComponents.Tmmc[i].Histogram.size(); j++)
534 {
535 size_t N = j / SystemComponents.Tmmc[i].nbinPerMacrostate;
536 size_t bin = j % SystemComponents.Tmmc[i].nbinPerMacrostate;
537 textTMMCFile << j << " " << N << " " << bin << " ";
538 textTMMCFile /*<< std::setprecision (15)*/ << SystemComponents.Tmmc[i].CMatrix[j].x << " " ;
539 textTMMCFile /*<< std::setprecision (15)*/ << SystemComponents.Tmmc[i].CMatrix[j].y << " " ;
540 textTMMCFile /*<< std::setprecision (15)*/ << SystemComponents.Tmmc[i].CMatrix[j].z << " " ;
541 textTMMCFile /*<< std::setprecision (15)*/ << SystemComponents.Tmmc[i].WLBias[j] << " " ;
542 textTMMCFile /*<< std::setprecision (15)*/ << SystemComponents.Tmmc[i].ln_g[j] << " ";
543 textTMMCFile /*<< std::setprecision (15)*/ << SystemComponents.Tmmc[i].TMBias[j] << " " ;
544 textTMMCFile /*<< std::setprecision (15)*/ << SystemComponents.Tmmc[i].lnpi[j] << " ";
545 textTMMCFile /*<< std::setprecision (15)*/ << SystemComponents.Tmmc[i].forward_lnpi[j] << " " ;
546 textTMMCFile /*<< std::setprecision (15)*/ << SystemComponents.Tmmc[i].reverse_lnpi[j] << " " ;
547 textTMMCFile << SystemComponents.Tmmc[i].Histogram[j] << '\n';
548 }
549 }
550 }
551 textTMMCFile.close();
552}
553
554inline void GenerateSummaryAtEnd(int Cycle, std::vector<Components>& SystemComponents, Simulations*& Sims, ForceField& FF, std::vector<Boxsize>& Box, PseudoAtomDefinitions& PseudoAtom)
555{
556
557 size_t NumberOfSimulations = SystemComponents.size();
558 for(size_t i = 0; i < NumberOfSimulations; i++)
559 {
560 printf("System %zu\n", i);
561 Atoms device_System[SystemComponents[i].Total_Components];
562 Copy_AtomData_from_Device(device_System, SystemComponents[i].HostSystem, Sims[i].d_a, SystemComponents[i]);
563
564 Write_Lambda(Cycle, SystemComponents[i], i);
565 Write_TMMC(Cycle, SystemComponents[i], i);
566 //Print Number of Pseudo Atoms//
567 for(size_t j = 0; j < SystemComponents[i].NumberOfPseudoAtoms.size(); j++) printf("PseudoAtom Type: %s[%zu], #: %zu\n", PseudoAtom.Name[j].c_str(), j, SystemComponents[i].NumberOfPseudoAtoms[j]);
568 }
569}
570
571inline void prepare_MixtureStats(Components& SystemComponents)
572{
573 double tot = 0.0;
574 printf("================= MOL FRACTIONS =================\n");
575 for(size_t j = 0; j < SystemComponents.Total_Components; j++)
576 {
577 SystemComponents.Moves[j].IdentitySwap_Total_TO.resize(SystemComponents.Total_Components, 0);
578 SystemComponents.Moves[j].IdentitySwap_Acc_TO.resize(SystemComponents.Total_Components, 0);
579 if(j != 0) tot += SystemComponents.MolFraction[j];
580 }
581 //Prepare MolFraction for adsorbate components//
582 for(size_t j = 1; j < SystemComponents.Total_Components; j++)
583 {
584 SystemComponents.MolFraction[j] /= tot;
585 printf("Component [%zu] (%s), Mol Fraction: %.5f\n", j, SystemComponents.MoleculeName[j].c_str(), SystemComponents.MolFraction[j]);
586 }
587 printf("=================================================\n");
588}
Atoms
Definition data_struct.h:746
Boxsize
Definition data_struct.h:819
Complex
Definition data_struct.h:46
Components
Definition data_struct.h:843
ForceField
Definition data_struct.h:794
LAMBDA
Definition data_struct.h:73
MoveEnergy
Definition data_struct.h:615
Move_Statistics
Definition data_struct.h:416
PseudoAtomDefinitions
Definition data_struct.h:767
RandomNumber
Definition data_struct.h:1053
Simulations
Definition data_struct.h:1027
TMMC
Definition data_struct.h:106
Units
Definition data_struct.h:52
WidomStruct
Definition data_struct.h:1044
Generated by doxygen 1.9.8