My Project
Loading...
Searching...
No Matches
mc_cbcfc.h
1__global__ void Prepare_LambdaChange(Atoms* d_a, Atoms Mol, Simulations Sims, ForceField FF, size_t start_position, size_t SelectedComponent, bool* device_flag)
2{
3 size_t i = blockIdx.x * blockDim.x + threadIdx.x;
4 size_t real_pos = start_position + i;
5
6 const double3 pos = d_a[SelectedComponent].pos[real_pos];
7 const double scale = d_a[SelectedComponent].scale[real_pos];
8 const double charge = d_a[SelectedComponent].charge[real_pos];
9 const double scaleCoul = d_a[SelectedComponent].scaleCoul[real_pos];
10 const size_t Type = d_a[SelectedComponent].Type[real_pos];
11 const size_t MolID = d_a[SelectedComponent].MolID[real_pos];
12
13 //Atoms Temp = Sims.Old; Atoms TempNEW = Sims.New;
14
15 Mol.pos[i] = pos;
16 Mol.scale[i] = scale; Mol.charge[i] = charge; Mol.scaleCoul[i] = scaleCoul; Mol.Type[i] = Type; Mol.MolID[i] = MolID;
17 device_flag[i] = false;
18}
19
20static inline MoveEnergy CBCF_LambdaChange(Components& SystemComponents, Simulations& Sims, ForceField& FF, RandomNumber& Random, WidomStruct& Widom, size_t SelectedMolInComponent, size_t SelectedComponent, double2 oldScale, double2 newScale, size_t& start_position, int MoveType, bool& SuccessConstruction)
21{
22 MoveEnergy tot;
23 //double result = 0.0;
24 size_t Atomsize = 0;
25 for(size_t ijk = 0; ijk < SystemComponents.Total_Components; ijk++)
26 {
27 Atomsize += SystemComponents.Moleculesize[ijk] * SystemComponents.NumberOfMolecule_for_Component[ijk];
28 }
29 size_t Molsize;
30 //Set up Old position and New position arrays
31 Molsize = SystemComponents.Moleculesize[SelectedComponent]; //Get the size of the selected Molecule
32 if(Molsize >= 1024)
33 {
34 throw std::runtime_error("Molecule size is greater than allocated size, Why so big?\n");
35 }
36 start_position = SelectedMolInComponent*SystemComponents.Moleculesize[SelectedComponent];
37 Prepare_LambdaChange<<<1, Molsize>>>(Sims.d_a, Sims.Old, Sims, FF, start_position, SelectedComponent, Sims.device_flag);
38
39 // Setup for the pairwise calculation //
40 size_t NHostAtom = 0; size_t NGuestAtom = 0;
41 for(size_t i = 0; i < SystemComponents.NComponents.y; i++)
42 NHostAtom += SystemComponents.Moleculesize[i] * SystemComponents.NumberOfMolecule_for_Component[i];
43 for(size_t i = SystemComponents.NComponents.y; i < SystemComponents.NComponents.x; i++)
44 NGuestAtom+= SystemComponents.Moleculesize[i] * SystemComponents.NumberOfMolecule_for_Component[i];
45
46
47 size_t HH_Nthread=0; size_t HH_Nblock=0; Setup_threadblock(NHostAtom * Molsize, &HH_Nblock, &HH_Nthread);
48 size_t HG_Nthread=0; size_t HG_Nblock=0; Setup_threadblock(NHostAtom * Molsize, &HG_Nblock, &HG_Nthread);
49 size_t GG_Nthread=0; size_t GG_Nblock=0; Setup_threadblock(NGuestAtom * Molsize, &GG_Nblock, &GG_Nthread);
50
51 size_t CrossTypeNthread = 0;
52 if(SelectedComponent < SystemComponents.NComponents.y) //Framework-Framework + Framework-Adsorbate//
53 {GG_Nthread = 0; GG_Nblock = 0; CrossTypeNthread = HH_Nthread; }
54 else //Framework-Adsorbate + Adsorbate-Adsorbate//
55 {HH_Nthread = 0; HH_Nblock = 0; CrossTypeNthread = GG_Nthread; }
56
57 size_t Nthread = std::max(CrossTypeNthread, HG_Nthread);
58 size_t Total_Nblock = HH_Nblock + HG_Nblock + GG_Nblock;
59
60 int3 NBlocks = {(int) HH_Nblock, (int) HG_Nblock, (int) GG_Nblock}; //x: HH_Nblock, y: HG_Nblock, z: GG_Nblock;
61 bool Do_New = true; bool Do_Old = true;
62
63 Calculate_Single_Body_Energy_SEPARATE_HostGuest_VDWReal_LambdaChange<<<Total_Nblock, Nthread, Nthread * 2 * sizeof(double)>>>(Sims.Box, Sims.d_a, Sims.Old, Sims.New, FF, Sims.Blocksum, SelectedComponent, Atomsize, Molsize, Sims.device_flag, NBlocks, Do_New, Do_Old, SystemComponents.NComponents, newScale);
64
65 cudaMemcpy(SystemComponents.flag, Sims.device_flag, sizeof(bool), cudaMemcpyDeviceToHost);
66
67 if(!SystemComponents.flag[0])
68 {
69 SuccessConstruction = true;
70 double BlockResult[Total_Nblock + Total_Nblock];
71 cudaMemcpy(BlockResult, Sims.Blocksum, 2 * Total_Nblock * sizeof(double), cudaMemcpyDeviceToHost);
72
73 //VDW Part and Real Part Coulomb//
74 for(size_t i = 0; i < HH_Nblock; i++)
75 {
76 tot.HHVDW += BlockResult[i];
77 tot.HHReal+= BlockResult[i + Total_Nblock];
78 //if(MoveType == SPECIAL_ROTATION) printf("HH Block %zu, VDW: %.5f, Real: %.5f\n", i, BlockResult[i], BlockResult[i + Total_Nblock]);
79 }
80 for(size_t i = HH_Nblock; i < HH_Nblock + HG_Nblock; i++)
81 {
82 tot.HGVDW += BlockResult[i];
83 tot.HGReal+= BlockResult[i + Total_Nblock];
84 //if(SystemComponents.CURRENTCYCLE == 25) printf("HG Block %zu, VDW: %.5f, Real: %.5f\n", i, BlockResult[i], BlockResult[i + Total_Nblock]);
85 }
86 for(size_t i = HH_Nblock + HG_Nblock; i < Total_Nblock; i++)
87 {
88 tot.GGVDW += BlockResult[i];
89 tot.GGReal+= BlockResult[i + Total_Nblock];
90 //printf("GG Block %zu, VDW: %.5f, Real: %.5f\n", i, BlockResult[i], BlockResult[i + Total_Nblock]);
91 }
92
93 if(!FF.noCharges && SystemComponents.hasPartialCharge[SelectedComponent])
94 {
95 //Zhao's note: since we changed it from using translation/rotation functions to its own, this needs to be changed as well//
96 double2 EwaldE = GPU_EwaldDifference_LambdaChange(Sims.Box, Sims.d_a, Sims.Old, FF, Sims.Blocksum, SystemComponents, SelectedComponent, oldScale, newScale, MoveType);
97 if(HH_Nblock == 0)
98 {
99 tot.GGEwaldE = EwaldE.x;
100 tot.HGEwaldE = EwaldE.y;
101 }
102 else
103 {
104 tot.HHEwaldE = EwaldE.x;
105 tot.HGEwaldE = EwaldE.y;
106 //printf("HHEwald: %.5f, HGEwald: %.5f\n", tot.HHEwaldE, tot.HGEwaldE);
107 }
108
109 }
110 tot.TailE = TailCorrectionDifference(SystemComponents, SelectedComponent, FF.size, Sims.Box.Volume, MoveType);
111 }
112 else
113 {
114 SuccessConstruction = false;
115 }
116 //printf("Cycle: %zu, Newscale: %.5f/%.5f, Oldscale: %.5f/%.5f, Construction: %s, CBCF Lambda E", SystemComponents.CURRENTCYCLE, newScale.x, newScale.y, oldScale.x, oldScale.y, SuccessConstruction ? "SUCCESS" : "BAD"); tot.print();
117 return tot;
118}
119
120__global__ void Revert_CBCF_Insertion(Atoms* d_a, size_t SelectedComponent, size_t start_position, double2 oldScale)
121{
122 size_t i = blockIdx.x * blockDim.x + threadIdx.x;
123 d_a[SelectedComponent].scale[start_position+i] = oldScale.x;
124 d_a[SelectedComponent].scaleCoul[start_position+i] = oldScale.y;
125}
126
128
129__global__ void Update_deletion_data_fractional(Atoms* d_a, size_t SelectedComponent, size_t UpdateLocation, int Moleculesize, size_t LastLocation)
130{
131 size_t i = blockIdx.x * blockDim.x + threadIdx.x;
132
133 //UpdateLocation should be the molecule that needs to be deleted
134 //Then move the atom at the last position to the location of the deleted molecule
135 //**Zhao's note** MolID of the deleted molecule should not be changed
136 //**Zhao's note** if Molecule deleted is the last molecule, then nothing is copied, just change the size later.
137 //**Zhao's note** for fractional molecule, since we still need the data of the deleted fractional molecule (for reversion if the move is rejected), store the old data at the LastLocation of d_a
138 if(UpdateLocation != LastLocation)
139 {
140 for(size_t i = 0; i < Moleculesize; i++)
141 {
142 double3 pos = d_a[SelectedComponent].pos[UpdateLocation+i];
143 double scale = d_a[SelectedComponent].scale[UpdateLocation+i];
144 double charge = d_a[SelectedComponent].charge[UpdateLocation+i];
145 double scaleCoul = d_a[SelectedComponent].scaleCoul[UpdateLocation+i]; double Type = d_a[SelectedComponent].Type[UpdateLocation+i];
146 d_a[SelectedComponent].pos[UpdateLocation+i] = d_a[SelectedComponent].pos[LastLocation+i];
147 d_a[SelectedComponent].scale[UpdateLocation+i] = d_a[SelectedComponent].scale[LastLocation+i];
148 d_a[SelectedComponent].charge[UpdateLocation+i] = d_a[SelectedComponent].charge[LastLocation+i];
149 d_a[SelectedComponent].scaleCoul[UpdateLocation+i] = d_a[SelectedComponent].scaleCoul[LastLocation+i];
150 d_a[SelectedComponent].Type[UpdateLocation+i] = d_a[SelectedComponent].Type[LastLocation+i];
151 //Put the data of the old fractional molecule to LastLocation//
152 d_a[SelectedComponent].pos[LastLocation+i] = pos;
153 d_a[SelectedComponent].scale[LastLocation+i] = scale;
154 d_a[SelectedComponent].charge[LastLocation+i] = charge; d_a[SelectedComponent].scaleCoul[LastLocation+i] = scaleCoul;
155 d_a[SelectedComponent].Type[LastLocation+i] = Type;
156 }
157 }
158 //Zhao's note: the single values in d_a and System are pointing to different locations//
159 //d_a is just device (cannot access on host), while System is shared (accessible on host), need to update d_a values here
160 //there are two of these values: size and Allocate_size
161 if(i==0)
162 {
163 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
164 }
165}
166
167__global__ void Revert_CBCF_Deletion(Atoms* d_a, Atoms NewMol, size_t SelectedComponent, size_t UpdateLocation, int Moleculesize, size_t LastLocation)
168{
169 //Zhao's note: this is literally the same function as Update_deletion_data_fractional
170 //UpdateLocation should be the molecule that needs to be deleted
171 //Then move the atom at the last position to the location of the deleted molecule
172 //**Zhao's note** MolID of the deleted molecule should not be changed
173 //**Zhao's note** if Molecule deleted is the last molecule, then nothing is copied, just change the size later.
174 //**Zhao's note** for fractional molecule, since we still need the data of the deleted fractional molecule (for reversion if the move is rejected), store the old data at the LastLocation of d_a
175 if(UpdateLocation != LastLocation)
176 {
177 for(size_t i = 0; i < Moleculesize; i++)
178 {
179 double3 pos = d_a[SelectedComponent].pos[UpdateLocation+i];
180 double scale = d_a[SelectedComponent].scale[UpdateLocation+i];
181 double charge = d_a[SelectedComponent].charge[UpdateLocation+i];
182 double scaleCoul = d_a[SelectedComponent].scaleCoul[UpdateLocation+i]; double Type = d_a[SelectedComponent].Type[UpdateLocation+i];
183 d_a[SelectedComponent].pos[UpdateLocation+i] = d_a[SelectedComponent].pos[LastLocation+i];
184 d_a[SelectedComponent].scale[UpdateLocation+i] = d_a[SelectedComponent].scale[LastLocation+i];
185 d_a[SelectedComponent].charge[UpdateLocation+i] = d_a[SelectedComponent].charge[LastLocation+i];
186 d_a[SelectedComponent].scaleCoul[UpdateLocation+i] = d_a[SelectedComponent].scaleCoul[LastLocation+i];
187 d_a[SelectedComponent].Type[UpdateLocation+i] = d_a[SelectedComponent].Type[LastLocation+i];
188 //Put the data of the old fractional molecule to LastLocation//
189 d_a[SelectedComponent].pos[LastLocation+i] = pos;
190 d_a[SelectedComponent].scale[LastLocation+i] = scale;
191 d_a[SelectedComponent].charge[LastLocation+i] = charge; d_a[SelectedComponent].scaleCoul[LastLocation+i] = scaleCoul;
192 d_a[SelectedComponent].Type[LastLocation+i] = Type;
193 }
194 }
195 //Zhao's note: the single values in d_a and System are pointing to different locations//
196 //d_a is just device (cannot access on host), while System is shared (accessible on host), need to update d_a values here
197 //there are two of these values: size and Allocate_size
198 size_t i = blockIdx.x * blockDim.x + threadIdx.x;
199 if(i==0)
200 {
201 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
202 }
203}
204
205__global__ void update_CBCF_scale(Atoms* d_a, size_t start_position, size_t SelectedComponent, double2 newScale)
206{
207 size_t i = blockIdx.x * blockDim.x + threadIdx.x;
208 d_a[SelectedComponent].scale[start_position+i] = newScale.x;
209 d_a[SelectedComponent].scaleCoul[start_position+i] = newScale.y;
210}
211
213
214MoveEnergy CBCFMove(Components& SystemComponents, Simulations& Sims, ForceField& FF, RandomNumber& Random, WidomStruct& Widom, size_t SelectedMolInComponent, size_t SelectedComponent);
215
216static inline MoveEnergy CBCFMove(Components& SystemComponents, Simulations& Sims, ForceField& FF, RandomNumber& Random, WidomStruct& Widom, size_t SelectedMolInComponent, size_t SelectedComponent)
217{
218 //Get Number of Molecules for this component (For updating TMMC)//
219 double NMol = SystemComponents.NumberOfMolecule_for_Component[SelectedComponent];
220 if(SystemComponents.hasfractionalMolecule[SelectedComponent]) NMol--;
221
222 int MoveType;
223 double TMMCPacc = 0.0;
224
225 SystemComponents.Moves[SelectedComponent].CBCFTotal ++;
226 //First draw a New Lambda, from the Bin//
227 LAMBDA lambda = SystemComponents.Lambda[SelectedComponent];
228 MoveEnergy final_energy;
229 size_t oldBin = lambda.currentBin;
230 size_t nbin = lambda.binsize;
231 double delta = lambda.delta;
232 double oldLambda = delta * static_cast<double>(oldBin);
233
234 //printf("Old Bin is %zu, Old Lambda is %.5f\n", oldBin, oldLambda);
235
236 MoveEnergy energy;
237
238 int SelectednewBin = selectNewBinTMMC(SystemComponents.Lambda[SelectedComponent]);
239 //int SelectednewBin = 1;
240 int Binchange = 1; //change for updating/getting bias for TMMC + CBCFC//
241 //Zhao's note: if the new bin == oldbin, the scale is the same, no need for a move, so select a new lambda
242 //Zhao's note: Do not delete the fractional molecule when there is only one molecule, the system needs to have at least one fractional molecule for each species
243 while(SelectednewBin == oldBin || (SelectednewBin < 0 && SystemComponents.NumberOfMolecule_for_Component[SelectedComponent] == 1))
244 {
245 if(SystemComponents.Tmmc[SelectedComponent].DoTMMC) //CBCFC + TMMC, can only move to adjacent bin//
246 {
247 //Binchange = 1;
248 //if(Get_Uniform_Random() < 0.5) SelectednewBin = -1;//Binchange = -1;
249 //SelectednewBin = Binchange + oldBin;
250 SelectednewBin = selectNewBinTMMC(SystemComponents.Lambda[SelectedComponent]);
251 //printf("in while, SelectednewBin: %d, oldBin: %zu\n", SelectednewBin, oldBin);
252 }
253 else //non-TMMC//
254 {
255 SelectednewBin = selectNewBin(SystemComponents.Lambda[SelectedComponent]);
256 }
257 }
258 Binchange = SelectednewBin - oldBin;
259 //printf("Binchange: %d, NewBin: %d, oldBin: %zu\n", Binchange, SelectednewBin, oldBin);
261 //INSERTION MOVE//
263 if(SelectednewBin > static_cast<int>(nbin))
264 {
265 MoveType = CBCF_INSERTION;
266
267 //return 0.0;
268 SystemComponents.Moves[SelectedComponent].CBCFInsertionTotal ++;
269 int newBin = SelectednewBin - static_cast<int>(nbin); //printf("INSERTION, newBin: %zu\n", newBin);
270 double newLambda = delta * static_cast<double>(newBin);
271 //Zhao's note: Set the Scaling factors for VDW and Coulombic Interaction depending on the mode//
272 double2 oldScale = SystemComponents.Lambda[SelectedComponent].SET_SCALE(oldLambda);
273 double2 InterScale= SystemComponents.Lambda[SelectedComponent].SET_SCALE(1.0);
274 double2 newScale = SystemComponents.Lambda[SelectedComponent].SET_SCALE(newLambda);
276 //First step: Increase the scale of the *current* fractional molecule to 1.0//
278 size_t start_position = 1;
279 bool SuccessConstruction = false;
280 energy = CBCF_LambdaChange(SystemComponents, Sims, FF, Random, Widom, SelectedMolInComponent, SelectedComponent, oldScale, InterScale, start_position, CBCF_INSERTION, SuccessConstruction);
281 if(!SuccessConstruction)
282 {
283 //If unsuccessful move (Overlap), Pacc = 0//
284 SystemComponents.Tmmc[SelectedComponent].UpdateCBCF(0.0, NMol, Binchange);
285 energy.zero();
286 return energy;
287 }
288 MoveEnergy first_step_energy = energy;
289 //Pretend the first step is accepted//
290 update_CBCF_scale<<<1,SystemComponents.Moleculesize[SelectedComponent]>>>(Sims.d_a, start_position, SelectedComponent, InterScale);
291 SuccessConstruction = false;
292 double Rosenbluth = 0.0;
293 size_t SelectedTrial = 0;
294 double preFactor = 0.0;
295 bool Accepted = false;
297 //Second step: Insertion of the fractional molecule//
299 MoveEnergy second_step_energy = Insertion_Body(SystemComponents, Sims, FF, Random, Widom, SelectedMolInComponent, SelectedComponent, Rosenbluth, SuccessConstruction, SelectedTrial, preFactor, true, newScale);
300 if(SuccessConstruction)
301 {
302 energy += second_step_energy;
303 //Account for the biasing terms of different lambdas//
304 double biasTerm = SystemComponents.Lambda[SelectedComponent].biasFactor[newBin] - SystemComponents.Lambda[SelectedComponent].biasFactor[oldBin];
305 preFactor *= std::exp(-SystemComponents.Beta * first_step_energy.total() + biasTerm);
306 double IdealRosen = SystemComponents.IdealRosenbluthWeight[SelectedComponent];
307 TMMCPacc = preFactor * Rosenbluth / IdealRosen; //Unbiased Acceptance//
308 //Apply the bias according to the macrostate//
309 SystemComponents.Tmmc[SelectedComponent].ApplyWLBiasCBCF(preFactor, NMol, Binchange);
310 SystemComponents.Tmmc[SelectedComponent].ApplyTMBiasCBCF(preFactor, NMol, Binchange);
311
312 if(Get_Uniform_Random() < preFactor * Rosenbluth / IdealRosen) Accepted = true;
313 SystemComponents.Tmmc[SelectedComponent].TreatAccOutofBoundCBCF(Accepted, NMol, Binchange);
314
315 //printf("Insertion E: %.5f, Acc: %s\n", energy.total(), Accepted ? "Accept" : "Reject");
316
317 if(Accepted)
318 { // accept the move
319 SystemComponents.Moves[SelectedComponent].CBCFInsertionAccepted ++;
320 SystemComponents.Moves[SelectedComponent].CBCFAccepted ++;
321 size_t UpdateLocation = SystemComponents.Moleculesize[SelectedComponent] * SystemComponents.NumberOfMolecule_for_Component[SelectedComponent];
322 //Zhao's note: here needs more consideration: need to update after implementing polyatomic molecule
323 Update_insertion_data<<<1,1>>>(Sims.d_a, Sims.Old, Sims.New, SelectedTrial, SelectedComponent, UpdateLocation, (int) SystemComponents.Moleculesize[SelectedComponent]);
324 Update_NumberOfMolecules(SystemComponents, Sims.d_a, SelectedComponent, CBCF_INSERTION);
325 //Update the ID of the fractional molecule on the host//
326 SystemComponents.Lambda[SelectedComponent].FractionalMoleculeID = SystemComponents.NumberOfMolecule_for_Component[SelectedComponent] - 1;
327 SystemComponents.Lambda[SelectedComponent].currentBin = newBin;
328 if(SystemComponents.Tmmc[SelectedComponent].DoTMMC)
329 SystemComponents.Tmmc[SelectedComponent].currentBin = newBin;
330 if(!FF.noCharges && SystemComponents.hasPartialCharge[SelectedComponent])
331 {
332 Update_Ewald_Vector(Sims.Box, false, SystemComponents, SelectedComponent);
333 }
334 final_energy = energy;
335 }
336 }
337 if(!Accepted) Revert_CBCF_Insertion<<<1, SystemComponents.Moleculesize[SelectedComponent]>>>(Sims.d_a, SelectedComponent, start_position, oldScale);
338 }
339 else if(SelectednewBin < 0) //Deletion Move//
340 {
341 MoveType = CBCF_DELETION;
342 //return 0.0; //for debugging
343 SystemComponents.Moves[SelectedComponent].CBCFDeletionTotal ++;
344 int newBin = SelectednewBin + nbin; //printf("DELETION, newBin: %zu\n", newBin);
345 double newLambda = delta * static_cast<double>(newBin);
346 //Zhao's note: Set the Scaling factors for VDW and Coulombic Interaction depending on the mode//
347 double2 oldScale = SystemComponents.Lambda[SelectedComponent].SET_SCALE(oldLambda);
348 double2 InterScale= SystemComponents.Lambda[SelectedComponent].SET_SCALE(1.0);
349 double2 newScale = SystemComponents.Lambda[SelectedComponent].SET_SCALE(newLambda);
351 //First step: delete the fractional molecule//
353 bool SuccessConstruction = false;
354 MoveEnergy energy;
355 double Rosenbluth = 0.0;
356 double preFactor = 0.0;
357 size_t UpdateLocation = SelectedMolInComponent * SystemComponents.Moleculesize[SelectedComponent];
358 energy = Deletion_Body(SystemComponents, Sims, FF, Random, Widom, SelectedMolInComponent, SelectedComponent, UpdateLocation, Rosenbluth, SuccessConstruction, preFactor, oldScale);
359
360 //printf("Deletion First Step E: "); energy.print();
361
362 //MoveEnergy first_step_energy = energy;
363 //Pretend the deletion move is accepted//
364 size_t LastMolecule = SystemComponents.NumberOfMolecule_for_Component[SelectedComponent]-1;
365 size_t LastLocation = LastMolecule*SystemComponents.Moleculesize[SelectedComponent];
366 Update_deletion_data_fractional<<<1,1>>>(Sims.d_a, SelectedComponent, UpdateLocation, (int) SystemComponents.Moleculesize[SelectedComponent], LastLocation);
367 //Record the old fractionalmolecule ID//
368 size_t OldFracMolID = SystemComponents.Lambda[SelectedComponent].FractionalMoleculeID;
369 Update_NumberOfMolecules(SystemComponents, Sims.d_a, SelectedComponent, CBCF_DELETION);
371 //Second step: randomly choose a new fractional molecule, perform the lambda move//
373 SelectedMolInComponent = static_cast<size_t>(Get_Uniform_Random() * static_cast<double>(SystemComponents.NumberOfMolecule_for_Component[SelectedComponent]));
374 SystemComponents.Lambda[SelectedComponent].FractionalMoleculeID = SelectedMolInComponent;
375
376 bool Accepted = false;
377 if(SuccessConstruction)
378 {
379 size_t start_position = 1;
380 MoveEnergy second_step_energy;
381
382 SuccessConstruction = false;
383 second_step_energy = CBCF_LambdaChange(SystemComponents, Sims, FF, Random, Widom, SelectedMolInComponent, SelectedComponent, InterScale, newScale, start_position, CBCF_DELETION, SuccessConstruction); //Fractional deletion, lambda change is the second step (that uses the temporary tempEik vector)
384
385 //Account for the biasing terms of different lambdas//
386 double biasTerm = SystemComponents.Lambda[SelectedComponent].biasFactor[newBin] - SystemComponents.Lambda[SelectedComponent].biasFactor[oldBin];
387 preFactor *= std::exp(-SystemComponents.Beta * second_step_energy.total() + biasTerm);
388 double IdealRosen = SystemComponents.IdealRosenbluthWeight[SelectedComponent];
389 TMMCPacc = preFactor * IdealRosen / Rosenbluth; //Unbiased Acceptance//
390 //Apply the bias according to the macrostate//
391 //SystemComponents.Tmmc[SelectedComponent].ApplyWLBiasCBCF(preFactor, NMol, Binchange);
392 SystemComponents.Tmmc[SelectedComponent].ApplyTMBiasCBCF(preFactor, NMol, Binchange);
393
394 if(Get_Uniform_Random() < preFactor * IdealRosen / Rosenbluth) Accepted = true;
395 SystemComponents.Tmmc[SelectedComponent].TreatAccOutofBoundCBCF(Accepted, NMol, Binchange);
396
397 if(Accepted)
398 {
399 SystemComponents.Moves[SelectedComponent].CBCFDeletionAccepted ++;
400 SystemComponents.Moves[SelectedComponent].CBCFAccepted ++;
401 //update_translation_position<<<1,SystemComponents.Moleculesize[SelectedComponent]>>>(Sims.d_a, Sims.New, start_position, SelectedComponent);
402 update_CBCF_scale<<<1,SystemComponents.Moleculesize[SelectedComponent]>>>(Sims.d_a, start_position, SelectedComponent, newScale);
403 SystemComponents.Lambda[SelectedComponent].currentBin = newBin;
404 if(SystemComponents.Tmmc[SelectedComponent].DoTMMC)
405 SystemComponents.Tmmc[SelectedComponent].currentBin = newBin;
406 if(!FF.noCharges && SystemComponents.hasPartialCharge[SelectedComponent])
407 {
408 Update_Ewald_Vector(Sims.Box, false, SystemComponents, SelectedComponent);
409 }
410 energy.take_negative();
411 energy += second_step_energy;
412 final_energy = energy;
413 }
414 }
415 if(!Accepted)
416 {
417 Revert_CBCF_Deletion<<<1,1>>>(Sims.d_a, Sims.New, SelectedComponent, UpdateLocation, (int) SystemComponents.Moleculesize[SelectedComponent], LastLocation);
418 SystemComponents.Lambda[SelectedComponent].FractionalMoleculeID = OldFracMolID;
419 SystemComponents.TotalNumberOfMolecules ++;
420 SystemComponents.NumberOfMolecule_for_Component[SelectedComponent] ++;
421 SystemComponents.UpdatePseudoAtoms(CBCF_INSERTION, SelectedComponent);
422 }
423 //if(SystemComponents.CURRENTCYCLE == 917) printf("Deletion E: %.5f, Acc: %s\n", energy.total(), Accepted ? "Accept": "Reject");
424 }
425 else //Lambda Move//
426 {
427 MoveType = CBCF_LAMBDACHANGE;
428 SystemComponents.Moves[SelectedComponent].CBCFLambdaTotal ++;
429 int newBin = static_cast<size_t>(SelectednewBin); //printf("LAMBDACHANGE, newBin: %zu\n", newBin);
430 double newLambda = delta * static_cast<double>(newBin);
431 //Zhao's note: Set the Scaling factors for VDW and Coulombic Interaction depending on the mode//
432 double2 oldScale = SystemComponents.Lambda[SelectedComponent].SET_SCALE(oldLambda);
433 double2 newScale = SystemComponents.Lambda[SelectedComponent].SET_SCALE(newLambda);
434 size_t start_position = 1;
435 bool SuccessConstruction = false;
436 energy = CBCF_LambdaChange(SystemComponents, Sims, FF, Random, Widom, SelectedMolInComponent, SelectedComponent, oldScale, newScale, start_position, CBCF_LAMBDACHANGE, SuccessConstruction);
437 if(!SuccessConstruction)
438 {
439 //If unsuccessful move (Overlap), Pacc = 0//
440 SystemComponents.Tmmc[SelectedComponent].UpdateCBCF(0.0, NMol, Binchange);
441 energy.zero();
442 return energy;
443 }
444 //Account for the biasing terms of different lambdas//
445 double biasTerm = SystemComponents.Lambda[SelectedComponent].biasFactor[newBin] - SystemComponents.Lambda[SelectedComponent].biasFactor[oldBin];
446 TMMCPacc = std::exp(-SystemComponents.Beta * energy.total() + biasTerm);
447 double preFactor = 1.0;
448 SystemComponents.Tmmc[SelectedComponent].ApplyTMBiasCBCF(preFactor, NMol, Binchange);
449
450 bool Accepted = false;
451 if (Get_Uniform_Random() < preFactor * std::exp(-SystemComponents.Beta * energy.total() + biasTerm)) Accepted = true;
452 SystemComponents.Tmmc[SelectedComponent].TreatAccOutofBoundCBCF(Accepted, NMol, Binchange);
453
454 if(Accepted)
455 {
456 SystemComponents.Moves[SelectedComponent].CBCFLambdaAccepted ++;
457 //update_translation_position<<<1,SystemComponents.Moleculesize[SelectedComponent]>>>(Sims.d_a, Sims.New, start_position, SelectedComponent);
458 update_CBCF_scale<<<1,SystemComponents.Moleculesize[SelectedComponent]>>>(Sims.d_a, start_position, SelectedComponent, newScale);
459 SystemComponents.Moves[SelectedComponent].CBCFAccepted ++;
460 SystemComponents.Lambda[SelectedComponent].currentBin = newBin;
461 if(SystemComponents.Tmmc[SelectedComponent].DoTMMC)
462 SystemComponents.Tmmc[SelectedComponent].currentBin = newBin;
463 if(!FF.noCharges && SystemComponents.hasPartialCharge[SelectedComponent])
464 {
465 Update_Ewald_Vector(Sims.Box, false, SystemComponents, SelectedComponent);
466 }
467 final_energy = energy;
468 }
469 }
470 SystemComponents.Tmmc[SelectedComponent].UpdateCBCF(TMMCPacc, NMol, Binchange);
471 return final_energy;
472}
Atoms
Definition data_struct.h:746
Components
Definition data_struct.h:843
ForceField
Definition data_struct.h:794
LAMBDA
Definition data_struct.h:73
MoveEnergy
Definition data_struct.h:615
RandomNumber
Definition data_struct.h:1053
Simulations
Definition data_struct.h:1027
WidomStruct
Definition data_struct.h:1044
Generated by doxygen 1.9.8