My Project
Loading...
Searching...
No Matches
mc_utilities.h
1
2// FUNCATIONS RELATED TO ALLOCATING MORE SPACE ON THE GPU //
4
5__global__ void AllocateMoreSpace_CopyToTemp(Atoms* d_a, Atoms temp, size_t Space, size_t SelectedComponent)
6{
7 for(size_t i = 0; i < Space; i++)
8 {
9 temp.pos[i] = d_a[SelectedComponent].pos[i];
10 temp.scale[i] = d_a[SelectedComponent].scale[i];
11 temp.charge[i] = d_a[SelectedComponent].charge[i];
12 temp.scaleCoul[i] = d_a[SelectedComponent].scaleCoul[i];
13 temp.Type[i] = d_a[SelectedComponent].Type[i];
14 temp.MolID[i] = d_a[SelectedComponent].MolID[i];
15 }
16}
17
18__global__ void AllocateMoreSpace_CopyBack(Atoms* d_a, Atoms temp, size_t Space, size_t Newspace, size_t SelectedComponent)
19{
20 d_a[SelectedComponent].Allocate_size = Newspace;
21 for(size_t i = 0; i < Space; i++)
22 {
23 d_a[SelectedComponent].pos[i] = temp.pos[i];
24 d_a[SelectedComponent].scale[i] = temp.scale[i];
25 d_a[SelectedComponent].charge[i] = temp.charge[i];
26 d_a[SelectedComponent].scaleCoul[i] = temp.scaleCoul[i];
27 d_a[SelectedComponent].Type[i] = temp.Type[i];
28 d_a[SelectedComponent].MolID[i] = temp.MolID[i];
29 /*
30 System_comp.pos[i] = temp.pos[i];
31 System_comp.scale[i] = temp.scale[i];
32 System_comp.charge[i] = temp.charge[i];
33 System_comp.scaleCoul[i] = temp.scaleCoul[i];
34 System_comp.Type[i] = temp.Type[i];
35 System_comp.MolID[i] = temp.MolID[i];
36 */
37 }
38 //test the new allocation//
39 //d_a[SelectedComponent].pos[Newspace-1].x = 0.0;
40 /*if(d_a[SelectedComponent].pos[Newspace-1].x = 0.0)
41 {
42 }*/
43}
44
45void AllocateMoreSpace(Atoms*& d_a, size_t SelectedComponent, Components& SystemComponents)
46{
47 printf("Allocating more space on device\n");
48 Atoms temp; // allocate a struct on the device for copying data.
49 //Atoms tempSystem[SystemComponents.Total_Components];
50 size_t Copysize=SystemComponents.Allocate_size[SelectedComponent];
51 size_t Morespace = 1024;
52 size_t Newspace = Copysize+Morespace;
53 //Allocate space on the temporary struct
54 cudaMalloc(&temp.pos, Copysize * sizeof(double3));
55 cudaMalloc(&temp.scale, Copysize * sizeof(double));
56 cudaMalloc(&temp.charge, Copysize * sizeof(double));
57 cudaMalloc(&temp.scaleCoul, Copysize * sizeof(double));
58 cudaMalloc(&temp.Type, Copysize * sizeof(size_t));
59 cudaMalloc(&temp.MolID, Copysize * sizeof(size_t));
60 // Copy data to temp
61 AllocateMoreSpace_CopyToTemp<<<1,1>>>(d_a, temp, Copysize, SelectedComponent);
62 // Allocate more space on the device pointers
63
64 Atoms System[SystemComponents.Total_Components]; cudaMemcpy(System, d_a, SystemComponents.Total_Components * sizeof(Atoms), cudaMemcpyDeviceToHost);
65
66 cudaMalloc(&System[SelectedComponent].pos, Newspace * sizeof(double3));
67 cudaMalloc(&System[SelectedComponent].scale, Newspace * sizeof(double));
68 cudaMalloc(&System[SelectedComponent].charge, Newspace * sizeof(double));
69 cudaMalloc(&System[SelectedComponent].scaleCoul, Newspace * sizeof(double));
70 cudaMalloc(&System[SelectedComponent].Type, Newspace * sizeof(size_t));
71 cudaMalloc(&System[SelectedComponent].MolID, Newspace * sizeof(size_t));
72 // Copy data from temp back to the new pointers
73 AllocateMoreSpace_CopyBack<<<1,1>>>(d_a, temp, Copysize, Newspace, SelectedComponent);
74 //System[SelectedComponent].Allocate_size = Newspace;
75}
76
78// FUNCTIONS RELATED TO ACCEPTING A MOVE //
80
81static inline void Update_NumberOfMolecules(Components& SystemComponents, Atoms*& d_a, size_t SelectedComponent, int MoveType)
82{
83 size_t Molsize = SystemComponents.Moleculesize[SelectedComponent]; //change in atom number counts
84 int NumMol = -1; //default: deletion; Insertion: +1, Deletion: -1, size_t is never negative
85
86 switch(MoveType)
87 {
88 case INSERTION: case SINGLE_INSERTION: case CBCF_INSERTION:
89 {
90 NumMol = 1; break;
91 }
92 case DELETION: case SINGLE_DELETION: case CBCF_DELETION:
93 {
94 NumMol = -1; break;
95 }
96 }
97 //Update Components
98 SystemComponents.NumberOfMolecule_for_Component[SelectedComponent] += NumMol;
99 SystemComponents.TotalNumberOfMolecules += NumMol;
100
101 // check System size //
102 size_t size = SystemComponents.Moleculesize[SelectedComponent] * SystemComponents.NumberOfMolecule_for_Component[SelectedComponent];
103
104 SystemComponents.UpdatePseudoAtoms(MoveType, SelectedComponent);
105
106 if(size > SystemComponents.Allocate_size[SelectedComponent])
107 {
108 AllocateMoreSpace(d_a, SelectedComponent, SystemComponents);
109 throw std::runtime_error("Need to allocate more space, not implemented\n");
110 }
111}
112
113__global__ void Update_SINGLE_INSERTION_data(Atoms* d_a, Atoms New, size_t SelectedComponent)
114{
115 //Assuming single thread//
116 size_t Molsize = d_a[SelectedComponent].Molsize;
117 size_t UpdateLocation = d_a[SelectedComponent].size;
118 for(size_t j = 0; j < Molsize; j++)
119 {
120 d_a[SelectedComponent].pos[UpdateLocation+j] = New.pos[j];
121 d_a[SelectedComponent].scale[UpdateLocation+j] = New.scale[j];
122 d_a[SelectedComponent].charge[UpdateLocation+j] = New.charge[j];
123 d_a[SelectedComponent].scaleCoul[UpdateLocation+j] = New.scaleCoul[j];
124 d_a[SelectedComponent].Type[UpdateLocation+j] = New.Type[j];
125 d_a[SelectedComponent].MolID[UpdateLocation+j] = New.MolID[j];
126 }
127 d_a[SelectedComponent].size += Molsize;
128}
129
130__global__ void Update_deletion_data(Atoms* d_a, size_t SelectedComponent, size_t UpdateLocation, int Moleculesize, size_t LastLocation)
131{
132 size_t i = blockIdx.x * blockDim.x + threadIdx.x;
133
134 //UpdateLocation should be the molecule that needs to be deleted
135 //Then move the atom at the last position to the location of the deleted molecule
136 //**Zhao's note** MolID of the deleted molecule should not be changed
137 //**Zhao's note** if Molecule deleted is the last molecule, then nothing is copied, just change the size later.
138 if(UpdateLocation != LastLocation)
139 {
140 for(size_t i = 0; i < Moleculesize; i++)
141 {
142 d_a[SelectedComponent].pos[UpdateLocation+i] = d_a[SelectedComponent].pos[LastLocation+i];
143 d_a[SelectedComponent].scale[UpdateLocation+i] = d_a[SelectedComponent].scale[LastLocation+i];
144 d_a[SelectedComponent].charge[UpdateLocation+i] = d_a[SelectedComponent].charge[LastLocation+i];
145 d_a[SelectedComponent].scaleCoul[UpdateLocation+i] = d_a[SelectedComponent].scaleCoul[LastLocation+i];
146 d_a[SelectedComponent].Type[UpdateLocation+i] = d_a[SelectedComponent].Type[LastLocation+i];
147 }
148 }
149 //Zhao's note: the single values in d_a and System are pointing to different locations//
150 //d_a is just device (cannot access on host), while System is shared (accessible on host), need to update d_a values here
151 //there are two of these values: size and Allocate_size
152 if(i==0)
153 {
154 d_a[SelectedComponent].size -= Moleculesize; //Zhao's special note: AllData.size doesn't work... So single values are painful, need to consider pointers for single values
155 }
156}
157
158__global__ void Update_insertion_data(Atoms* d_a, Atoms Mol, Atoms NewMol, size_t SelectedTrial, size_t SelectedComponent, size_t UpdateLocation, int Moleculesize)
159{
160 size_t i = blockIdx.x * blockDim.x + threadIdx.x;
161 //UpdateLocation should be the last position of the dataset
162 //Need to check if Allocate_size is smaller than size
163 if(Moleculesize == 1) //Only first bead is inserted, first bead data is stored in NewMol
164 {
165 d_a[SelectedComponent].pos[UpdateLocation] = NewMol.pos[SelectedTrial];
166 d_a[SelectedComponent].scale[UpdateLocation] = NewMol.scale[SelectedTrial];
167 d_a[SelectedComponent].charge[UpdateLocation] = NewMol.charge[SelectedTrial];
168 d_a[SelectedComponent].scaleCoul[UpdateLocation] = NewMol.scaleCoul[SelectedTrial];
169 d_a[SelectedComponent].Type[UpdateLocation] = NewMol.Type[SelectedTrial];
170 d_a[SelectedComponent].MolID[UpdateLocation] = NewMol.MolID[SelectedTrial];
171 }
172 else //Multiple beads: first bead + trial orientations
173 {
174 //Update the first bead, first bead data stored in position 0 of Mol //
175 d_a[SelectedComponent].pos[UpdateLocation] = Mol.pos[0];
176 d_a[SelectedComponent].scale[UpdateLocation] = Mol.scale[0];
177 d_a[SelectedComponent].charge[UpdateLocation] = Mol.charge[0];
178 d_a[SelectedComponent].scaleCoul[UpdateLocation] = Mol.scaleCoul[0];
179 d_a[SelectedComponent].Type[UpdateLocation] = Mol.Type[0];
180 d_a[SelectedComponent].MolID[UpdateLocation] = Mol.MolID[0];
181
182 size_t chainsize = Moleculesize - 1; // For trial orientations //
183 for(size_t j = 0; j < chainsize; j++) //Update the selected orientations//
184 {
185 size_t selectsize = SelectedTrial*chainsize+j;
186 d_a[SelectedComponent].pos[UpdateLocation+j+1] = NewMol.pos[selectsize];
187 d_a[SelectedComponent].scale[UpdateLocation+j+1] = NewMol.scale[selectsize];
188 d_a[SelectedComponent].charge[UpdateLocation+j+1] = NewMol.charge[selectsize];
189 d_a[SelectedComponent].scaleCoul[UpdateLocation+j+1] = NewMol.scaleCoul[selectsize];
190 d_a[SelectedComponent].Type[UpdateLocation+j+1] = NewMol.Type[selectsize];
191 d_a[SelectedComponent].MolID[UpdateLocation+j+1] = NewMol.MolID[selectsize];
192 }
193 }
194 //Zhao's note: the single values in d_a and System are pointing to different locations//
195 //d_a is just device (cannot access on host), while System is shared (accessible on host), need to update d_a values here
196 //there are two of these values: size and Allocate_size
197 if(i==0)
198 {
199 d_a[SelectedComponent].size += Moleculesize;
200 /*
201 for(size_t j = 0; j < Moleculesize; j++)
202 printf("xyz: %.5f %.5f %.5f, scale/charge/scaleCoul: %.5f %.5f %.5f, Type: %lu, MolID: %lu\n", d_a[SelectedComponent].pos[UpdateLocation+j].x, d_a[SelectedComponent].pos[UpdateLocation+j].y, d_a[SelectedComponent].pos[UpdateLocation+j].z, d_a[SelectedComponent].scale[UpdateLocation+j], d_a[SelectedComponent].charge[UpdateLocation+j], d_a[SelectedComponent].scaleCoul[UpdateLocation+j], d_a[SelectedComponent].Type[UpdateLocation+j], d_a[SelectedComponent].MolID[UpdateLocation+j]);
203 */
204 }
205}
206
207__global__ void update_translation_position(Atoms* d_a, Atoms NewMol, size_t start_position, size_t SelectedComponent)
208{
209 size_t i = blockIdx.x * blockDim.x + threadIdx.x;
210 d_a[SelectedComponent].pos[start_position+i] = NewMol.pos[i];
211 d_a[SelectedComponent].scale[start_position+i] = NewMol.scale[i];
212 d_a[SelectedComponent].charge[start_position+i] = NewMol.charge[i];
213 d_a[SelectedComponent].scaleCoul[start_position+i] = NewMol.scaleCoul[i];
214}
215
217// GET PREFACTOR FOR INSERTION/DELETION MOVES //
219
220static inline double GetPrefactor(Components& SystemComponents, Simulations& Sims, size_t SelectedComponent, int MoveType)
221{
222 double MolFraction = SystemComponents.MolFraction[SelectedComponent];
223 double FugacityCoefficient = SystemComponents.FugacityCoeff[SelectedComponent];
224 double NumberOfMolecules = static_cast<double>(SystemComponents.NumberOfMolecule_for_Component[SelectedComponent]);
225
226 //If component has fractional molecule, subtract the number of molecules by 1.//
227 if(SystemComponents.hasfractionalMolecule[SelectedComponent]){NumberOfMolecules-=1.0;}
228 if(NumberOfMolecules < 0.0) NumberOfMolecules = 0.0;
229
230 double preFactor = 0.0;
231 //Hard to generalize the prefactors, when you consider identity swap
232 //Certainly you can use insertion/deletion for identity swap, but you may lose accuracies
233 //since you are multiplying then dividing by the same variables
234 switch(MoveType)
235 {
236 case INSERTION: case SINGLE_INSERTION:
237 {
238 preFactor = SystemComponents.Beta * MolFraction * Sims.Box.Pressure * FugacityCoefficient * Sims.Box.Volume / (1.0+NumberOfMolecules);
239 break;
240 }
241 case DELETION: case SINGLE_DELETION:
242 {
243 preFactor = (NumberOfMolecules) / (SystemComponents.Beta * MolFraction * Sims.Box.Pressure * FugacityCoefficient * Sims.Box.Volume);
244 break;
245 }
246 case IDENTITY_SWAP:
247 {
248 throw std::runtime_error("Sorry, but no IDENTITY_SWAP OPTION for now. That move uses its own function\n");
249 }
250 case TRANSLATION: case ROTATION: case SPECIAL_ROTATION:
251 {
252 preFactor = 1.0;
253 }
254 }
255 return preFactor;
256}
257
259// ACCEPTION OF MOVES //
261
262static inline void AcceptInsertion(Components& SystemComponents, Simulations& Sims, size_t SelectedComponent, size_t SelectedTrial, bool noCharges, int MoveType)
263{
264 size_t UpdateLocation = SystemComponents.Moleculesize[SelectedComponent] * SystemComponents.NumberOfMolecule_for_Component[SelectedComponent];
265 //printf("AccInsertion, SelectedTrial: %zu, UpdateLocation: %zu\n", SelectedTrial, UpdateLocation);
266 //Zhao's note: here needs more consideration: need to update after implementing polyatomic molecule
267 if(MoveType == INSERTION)
268 {
269 Update_insertion_data<<<1,1>>>(Sims.d_a, Sims.Old, Sims.New, SelectedTrial, SelectedComponent, UpdateLocation, (int) SystemComponents.Moleculesize[SelectedComponent]);
270 }
271 else if(MoveType == SINGLE_INSERTION)
272 {
273 Update_SINGLE_INSERTION_data<<<1,1>>>(Sims.d_a, Sims.New, SelectedComponent);
274 }
275 Update_NumberOfMolecules(SystemComponents, Sims.d_a, SelectedComponent, INSERTION); //true = Insertion//
276 if(!noCharges && SystemComponents.hasPartialCharge[SelectedComponent])
277 {
278 Update_Ewald_Vector(Sims.Box, false, SystemComponents, SelectedComponent);
279 }
280}
281
282static inline void AcceptDeletion(Components& SystemComponents, Simulations& Sims, size_t SelectedComponent, size_t UpdateLocation, size_t SelectedMol, bool noCharges)
283{
284 size_t LastMolecule = SystemComponents.NumberOfMolecule_for_Component[SelectedComponent]-1;
285 size_t LastLocation = LastMolecule*SystemComponents.Moleculesize[SelectedComponent];
286 Update_deletion_data<<<1,1>>>(Sims.d_a, SelectedComponent, UpdateLocation, (int) SystemComponents.Moleculesize[SelectedComponent], LastLocation);
287
288 Update_NumberOfMolecules(SystemComponents, Sims.d_a, SelectedComponent, DELETION); //false = Deletion//
289 if(!noCharges && SystemComponents.hasPartialCharge[SelectedComponent])
290 {
291 Update_Ewald_Vector(Sims.Box, false, SystemComponents, SelectedComponent);
292 }
293 //Zhao's note: the last molecule can be the fractional molecule, (fractional molecule ID is stored on the host), we need to update it as well (at least check it)//
294 //The function below will only be processed if the system has a fractional molecule and the transfered molecule is NOT the fractional one //
295 if((SystemComponents.hasfractionalMolecule[SelectedComponent])&&(LastMolecule == SystemComponents.Lambda[SelectedComponent].FractionalMoleculeID))
296 {
297 //Since the fractional molecule is moved to the place of the selected deleted molecule, update fractional molecule ID on host
298 SystemComponents.Lambda[SelectedComponent].FractionalMoleculeID = SelectedMol;
299 }
300}
301
303// PREPARING NEW (TRIAL) LOCATIONS/ORIENTATIONS //
305
306__device__ void Rotate_Quaternions(double3 &Vec, double3 RANDOM)
307{
308 //https://stackoverflow.com/questions/31600717/how-to-generate-a-random-quaternion-quickly
309 //https://stackoverflow.com/questions/38978441/creating-uniform-random-quaternion-and-multiplication-of-two-quaternions
310 //James J. Kuffner 2004//
311 //Zhao's note: this one needs three random numbers//
312 const double u = RANDOM.x;
313 const double v = RANDOM.y;
314 const double w = RANDOM.z;
315 const double pi=3.14159265358979323846;//Zhao's note: consider using M_PI
316 //Sadly, no double4 type available//
317 const double q0 = sqrt(1-u) * std::sin(2*pi*v);
318 const double q1 = sqrt(1-u) * std::cos(2*pi*v);
319 const double q2 = sqrt(u) * std::sin(2*pi*w);
320 const double q3 = sqrt(u) * std::cos(2*pi*w);
321
322 double rot[3*3];
323 const double a01=q0*q1; const double a02=q0*q2; const double a03=q0*q3;
324 const double a11=q1*q1; const double a12=q1*q2; const double a13=q1*q3;
325 const double a22=q2*q2; const double a23=q2*q3; const double a33=q3*q3;
326
327 rot[0]=1.0-2.0*(a22+a33);
328 rot[1]=2.0*(a12-a03);
329 rot[2]=2.0*(a13+a02);
330 rot[3]=2.0*(a12+a03);
331 rot[4]=1.0-2.0*(a11+a33);
332 rot[5]=2.0*(a23-a01);
333 rot[6]=2.0*(a13-a02);
334 rot[7]=2.0*(a23+a01);
335 rot[8]=1.0-2.0*(a11+a22);
336 const double r=Vec.x*rot[0*3+0] + Vec.y*rot[0*3+1] + Vec.z*rot[0*3+2];
337 const double s=Vec.x*rot[1*3+0] + Vec.y*rot[1*3+1] + Vec.z*rot[1*3+2];
338 const double c=Vec.x*rot[2*3+0] + Vec.y*rot[2*3+1] + Vec.z*rot[2*3+2];
339 Vec={r, s, c};
340}
341
342__device__ void RotationAroundAxis(double3* pos, size_t i, double theta, double3 Axis)
343{
344 double w,s,c,rot[3*3];
345
346 c=cos(theta);
347 w=1.0-c;
348 s=sin(theta);
349
350 rot[0*3+0]=(Axis.x)*(Axis.x)*w+c;
351 rot[0*3+1]=(Axis.x)*(Axis.y)*w+(Axis.z)*s;
352 rot[0*3+2]=(Axis.x)*(Axis.z)*w-(Axis.y)*s;
353 rot[1*3+0]=(Axis.x)*(Axis.y)*w-(Axis.z)*s;
354 rot[1*3+1]=(Axis.y)*(Axis.y)*w+c;
355 rot[1*3+2]=(Axis.y)*(Axis.z)*w+(Axis.x)*s;
356 rot[2*3+0]=(Axis.x)*(Axis.z)*w+(Axis.y)*s;
357 rot[2*3+1]=(Axis.y)*(Axis.z)*w-(Axis.x)*s;
358 rot[2*3+2]=(Axis.z)*(Axis.z)*w+c;
359
360 w=pos[i].x*rot[0*3+0]+pos[i].y*rot[0*3+1]+pos[i].z*rot[0*3+2];
361 s=pos[i].x*rot[1*3+0]+pos[i].y*rot[1*3+1]+pos[i].z*rot[1*3+2];
362 c=pos[i].x*rot[2*3+0]+pos[i].y*rot[2*3+1]+pos[i].z*rot[2*3+2];
363 pos[i].x=w;
364 pos[i].y=s;
365 pos[i].z=c;
366}
367
368__global__ void get_new_position(Simulations& Sim, ForceField FF, size_t start_position, size_t SelectedComponent, double3 MaxChange, double3* RANDOM, size_t index, int MoveType)
369{
370 const size_t i = blockIdx.x * blockDim.x + threadIdx.x;
371 const size_t real_pos = start_position + i;
372
373 const double3 pos = Sim.d_a[SelectedComponent].pos[real_pos];
374 const double scale = Sim.d_a[SelectedComponent].scale[real_pos];
375 const double charge = Sim.d_a[SelectedComponent].charge[real_pos];
376 const double scaleCoul = Sim.d_a[SelectedComponent].scaleCoul[real_pos];
377 const size_t Type = Sim.d_a[SelectedComponent].Type[real_pos];
378 const size_t MolID = Sim.d_a[SelectedComponent].MolID[real_pos];
379
380 switch(MoveType)
381 {
382 case TRANSLATION://TRANSLATION:
383 {
384 Sim.New.pos[i] = pos + MaxChange * 2.0 * (RANDOM[index] - 0.5);
385 Sim.New.scale[i] = scale;
386 Sim.New.charge[i] = charge;
387 Sim.New.scaleCoul[i] = scaleCoul;
388 Sim.New.Type[i] = Type;
389 Sim.New.MolID[i] = MolID;
390
391 Sim.Old.pos[i] = pos;
392 Sim.Old.scale[i] = scale;
393 Sim.Old.charge[i] = charge;
394 Sim.Old.scaleCoul[i] = scaleCoul;
395 Sim.Old.Type[i] = Type;
396 Sim.Old.MolID[i] = MolID;
397 break;
398 }
399 case ROTATION://ROTATION:
400 {
401 Sim.New.pos[i] = pos - Sim.d_a[SelectedComponent].pos[start_position];
402 const double3 Angle = MaxChange * 2.0 * (RANDOM[index] - 0.5);
403
404 RotationAroundAxis(Sim.New.pos, i, Angle.x, {1.0, 0.0, 0.0}); //X
405 RotationAroundAxis(Sim.New.pos, i, Angle.y, {0.0, 1.0, 0.0}); //Y
406 RotationAroundAxis(Sim.New.pos, i, Angle.z, {0.0, 0.0, 1.0}); //Z
407
408 Sim.New.pos[i] += Sim.d_a[SelectedComponent].pos[start_position];
409
410 Sim.New.scale[i] = scale;
411 Sim.New.charge[i] = charge;
412 Sim.New.scaleCoul[i] = scaleCoul;
413 Sim.New.Type[i] = Type;
414 Sim.New.MolID[i] = MolID;
415
416 Sim.Old.pos[i] = pos;
417 Sim.Old.scale[i] = scale;
418 Sim.Old.charge[i] = charge;
419 Sim.Old.scaleCoul[i] = scaleCoul;
420 Sim.Old.Type[i] = Type;
421 Sim.Old.MolID[i] = MolID;
422 break;
423 }
424 case SINGLE_INSERTION:
425 {
426 //First ROTATION using QUATERNIONS//
427 //Then TRANSLATION//
428 double3 BoxLength = {Sim.Box.Cell[0], Sim.Box.Cell[4], Sim.Box.Cell[8]};
429 double3 NEW_COM = BoxLength * RANDOM[index];
430 if(i == 0) Sim.New.pos[0] = NEW_COM;
431 if(i > 0)
432 {
433 double3 Vec = pos - Sim.d_a[SelectedComponent].pos[start_position];
434 Rotate_Quaternions(Vec, RANDOM[index + 1]);
435 Sim.New.pos[i] = Vec + NEW_COM;
436 }
437 Sim.New.scale[i] = scale;
438 Sim.New.charge[i] = charge;
439 Sim.New.scaleCoul[i] = scaleCoul;
440 Sim.New.Type[i] = Type;
441 Sim.New.MolID[i] = Sim.d_a[SelectedComponent].size / Sim.d_a[SelectedComponent].Molsize;
442 break;
443 }
444 case SINGLE_DELETION: //Just Copy the old positions//
445 {
446 Sim.Old.pos[i] = pos;
447 Sim.Old.scale[i] = scale;
448 Sim.Old.charge[i] = charge;
449 Sim.Old.scaleCoul[i] = scaleCoul;
450 Sim.Old.Type[i] = Type;
451 Sim.Old.MolID[i] = MolID;
452 }
453
454 case SPECIAL_ROTATION:
455 {
456 // Rotation around axis (first/second atom) in the molecule //
457 Sim.New.pos[i] = pos - Sim.d_a[SelectedComponent].pos[start_position];
458
459 const double3 Angle = MaxChange * 2.0 * (RANDOM[index] - 0.5);
460 //Take direction, normalize//
461 double3 Axis;
462 //if(i == 0)
463 //{
464 Axis = Sim.d_a[SelectedComponent].pos[start_position + 1] - Sim.d_a[SelectedComponent].pos[start_position];
465 double norm = sqrt(dot(Axis, Axis));
466 Axis *= 1.0/norm;
467 //printf("Atom %lu, DistVec: %.5f %.5f %.5f\n", i, Sim.New.pos[i].x, Sim.New.pos[i].y, Sim.New.pos[i].z);
468 //}
469 if(i != 0 && i != 1) //Do not rotate the end points of the vector
470 RotationAroundAxis(Sim.New.pos, i, 3.0 * Angle.x, Axis);
471 Sim.New.pos[i] += Sim.d_a[SelectedComponent].pos[start_position];
472
473 Sim.New.scale[i] = scale;
474 Sim.New.charge[i] = charge;
475 Sim.New.scaleCoul[i] = scaleCoul;
476 Sim.New.Type[i] = Type;
477 Sim.New.MolID[i] = MolID;
478
479 Sim.Old.pos[i] = pos;
480 Sim.Old.scale[i] = scale;
481 Sim.Old.charge[i] = charge;
482 Sim.Old.scaleCoul[i] = scaleCoul;
483 Sim.Old.Type[i] = Type;
484 Sim.Old.MolID[i] = MolID;
485 //printf("Atom %lu, Angle: %.5f, Axis: %.5f %.5f %.5f, xyz: %.5f %.5f %.5f\n", i, Angle.x, Axis.x, Axis.y, Axis.z, Sim.New.pos[i].x, Sim.New.pos[i].y, Sim.New.pos[i].z);
486 }
487 }
488 Sim.device_flag[i] = false;
489}
490
492// OPTIMIZING THE ACCEPTANCE (MAINTAIN AROUND 0.5) FOR TRANSLATION/ROTATION MOVES //
494
495static inline void Update_Max_Translation(Components& SystemComponents, size_t Comp)
496{
497 if(SystemComponents.Moves[Comp].TranslationTotal == 0) return;
498 SystemComponents.Moves[Comp].TranslationAccRatio = static_cast<double>(SystemComponents.Moves[Comp].TranslationAccepted)/SystemComponents.Moves[Comp].TranslationTotal;
499 //printf("AccRatio is %.10f\n", MoveStats.TranslationAccRatio);
500 if(SystemComponents.Moves[Comp].TranslationAccRatio > 0.5)
501 {
502 SystemComponents.MaxTranslation[Comp] *= 1.05;
503 }
504 else
505 {
506 SystemComponents.MaxTranslation[Comp] *= 0.95;
507 }
508 if(SystemComponents.MaxTranslation[Comp].x < 0.01) SystemComponents.MaxTranslation[Comp].x = 0.01;
509 if(SystemComponents.MaxTranslation[Comp].y < 0.01) SystemComponents.MaxTranslation[Comp].y = 0.01;
510 if(SystemComponents.MaxTranslation[Comp].z < 0.01) SystemComponents.MaxTranslation[Comp].z = 0.01;
511
512 if(SystemComponents.MaxTranslation[Comp].x > 5.0) SystemComponents.MaxTranslation[Comp].x = 5.0;
513 if(SystemComponents.MaxTranslation[Comp].y > 5.0) SystemComponents.MaxTranslation[Comp].y = 5.0;
514 if(SystemComponents.MaxTranslation[Comp].z > 5.0) SystemComponents.MaxTranslation[Comp].z = 5.0;
515 SystemComponents.Moves[Comp].TranslationAccepted = 0;
516 SystemComponents.Moves[Comp].TranslationTotal = 0;
517}
518
519static inline void Update_Max_Rotation(Components& SystemComponents, size_t Comp)
520{
521 if(SystemComponents.Moves[Comp].RotationTotal == 0) return;
522 SystemComponents.Moves[Comp].RotationAccRatio = static_cast<double>(SystemComponents.Moves[Comp].RotationAccepted)/SystemComponents.Moves[Comp].RotationTotal;
523 //printf("AccRatio is %.10f\n", MoveStats.RotationAccRatio);
524 if(SystemComponents.Moves[Comp].RotationAccRatio > 0.5)
525 {
526 SystemComponents.MaxRotation[Comp] *= 1.05;
527 }
528 else
529 {
530 SystemComponents.MaxRotation[Comp] *= 0.95;
531 }
532 if(SystemComponents.MaxRotation[Comp].x < 0.01) SystemComponents.MaxRotation[Comp].x = 0.01;
533 if(SystemComponents.MaxRotation[Comp].y < 0.01) SystemComponents.MaxRotation[Comp].y = 0.01;
534 if(SystemComponents.MaxRotation[Comp].z < 0.01) SystemComponents.MaxRotation[Comp].z = 0.01;
535
536 if(SystemComponents.MaxRotation[Comp].x > 3.14) SystemComponents.MaxRotation[Comp].x = 3.14;
537 if(SystemComponents.MaxRotation[Comp].y > 3.14) SystemComponents.MaxRotation[Comp].y = 3.14;
538 if(SystemComponents.MaxRotation[Comp].z > 3.14) SystemComponents.MaxRotation[Comp].z = 3.14;
539 SystemComponents.Moves[Comp].RotationAccepted = 0;
540 SystemComponents.Moves[Comp].RotationTotal = 0;
541}
542
543static inline void Update_Max_SpecialRotation(Components& SystemComponents, size_t Comp)
544{
545 if(SystemComponents.Moves[Comp].SpecialRotationTotal == 0) return;
546 SystemComponents.Moves[Comp].SpecialRotationAccRatio = static_cast<double>(SystemComponents.Moves[Comp].SpecialRotationAccepted)/SystemComponents.Moves[Comp].SpecialRotationTotal;
547 //printf("AccRatio is %.10f\n", MoveStats.RotationAccRatio);
548 if(SystemComponents.Moves[Comp].SpecialRotationAccRatio > 0.5)
549 {
550 SystemComponents.MaxSpecialRotation[Comp] *= 1.05;
551 }
552 else
553 {
554 SystemComponents.MaxSpecialRotation[Comp] *= 0.95;
555 }
556 if(SystemComponents.MaxSpecialRotation[Comp].x < 0.01) SystemComponents.MaxSpecialRotation[Comp].x = 0.01;
557 if(SystemComponents.MaxSpecialRotation[Comp].y < 0.01) SystemComponents.MaxSpecialRotation[Comp].y = 0.01;
558 if(SystemComponents.MaxSpecialRotation[Comp].z < 0.01) SystemComponents.MaxSpecialRotation[Comp].z = 0.01;
559
560 if(SystemComponents.MaxSpecialRotation[Comp].x > 3.14) SystemComponents.MaxSpecialRotation[Comp].x = 3.14;
561 if(SystemComponents.MaxSpecialRotation[Comp].y > 3.14) SystemComponents.MaxSpecialRotation[Comp].y = 3.14;
562 if(SystemComponents.MaxSpecialRotation[Comp].z > 3.14) SystemComponents.MaxSpecialRotation[Comp].z = 3.14;
563 SystemComponents.Moves[Comp].SpecialRotationAccepted = 0;
564 SystemComponents.Moves[Comp].SpecialRotationTotal = 0;
565}
Atoms
Definition data_struct.h:746
Components
Definition data_struct.h:843
ForceField
Definition data_struct.h:794
Simulations
Definition data_struct.h:1027
Generated by doxygen 1.9.8