Questions and answers regarding FEASST
42 views
Skip to first unread message
harold...@nist.gov
Nov 6, 2022, 9:14:20 AM11/6/22
to feasst
Most users contact me directly, but I thought I would include my answers to some questions that may be of general use (and may be soon incorported into the documentation..). The user questions are numbered, with my responses beginning with HWH.
1. How do you choose the value of the alpha parameter? I know it should be about 5/L, for SPC/E you have used 5.6/L and for RPM simulation it was 6.87/L. Is there a specific rule indicating how to tune this value?
HWH> There is no specific rule, and to my knowledge the tuning of this parameter in the grand canonical ensemble is an open question of research. For example, there have been some studies that mostly focus on a fixed density, and sometimes for a particular system (e.g., water or biomolecules). The "tolerance" and "tolerance_num_sites" arguments for Ewald use the same formula as LAMMPS.
2. When using the ChargeScreened class, the erfc_table_size argument is used. I don't quite understand the meaning of this argument. On what basis do you choose the size?
HWH> See https://pages.nist.gov/feasst/plugin/charge/doc/ChargeScreened.html . That argument can be set to 0 if you wish to disable tabulation of erfc. It is tabulated for optimization (as also done in a similar way in LAMMPS, etc).
3. For what purpose is RefPotential used in the tutorial for SPC/E? For what purpose is DCCB used here and why is the cut off value 0.9*3.165 (as I understand it is the value of the oxygen sigma parameter).
HWH> DCCB is briefly mentioned here: https://pages.nist.gov/feasst/plugin/monte_carlo/doc/TrialStage.html . The DCCB cutoff value is relatively arbitrary but can be optimized. In this case, DCCB essentially uses a hard sphere reference potential to quickly find cavities (with the help of a cell list, and considering only oxygen sites). I used 90% of the sigma value here because SPC/E h-bonds can lead to a surprising amount of overlap for the Oxygen LJ term.
4. Is the use of ChargeScreenedIntra necessary whenever we model molecules that have bonds and use the Ewald sum?
HWH> Yes, this term cancels the intramolecular O-H charge term computed in Fourier space.
5. Using dummy atoms containing only charge, is it necessary to apply any additional commands using Ewald sum?
HWH> A dummy atom requires you to modify the TrialGrowFile or perhaps simply use TrialTransfer. I happened to be working on adding an example TIP4P simulation and can follow up with you on this soon.
6. How do you choose the tolerance value when determining the last macrostate (for example, the maximum number of molecules)? In your paper (https://doi.org/10.1063/1.5123683) you used a value of 25.0, for example in J. R. Errington's paper (https://doi.org/10.1063/1.1572463) it was a value of 10.0. Is there any specific rule based on which to determine the value of this tolerance?
HWH> The variation you see in the last macrostate is mainly because we don't a priori know what N_max value is necessary. As an exercise, you can perform your analysis (say, VLE density calculation) with varying N_max. You may find that you can get away with less than 10, because probabilities e^-10 don't contribute very much in the GCE averages. We often go a bit higher in N than necessary just because we don't know when we start the simulation. Depending on your familiarity with the system and how quickly you want your results, sometimes its worth parallelizing into windows with a N_max that you know is too large and then simply throw out the largest density results that had prohibitively small acceptance and therefore did not converge well.
1. How do you choose the value of the alpha parameter? I know it should be about 5/L, for SPC/E you have used 5.6/L and for RPM simulation it was 6.87/L. Is there a specific rule indicating how to tune this value?
HWH> There is no specific rule, and to my knowledge the tuning of this parameter in the grand canonical ensemble is an open question of research. For example, there have been some studies that mostly focus on a fixed density, and sometimes for a particular system (e.g., water or biomolecules). The "tolerance" and "tolerance_num_sites" arguments for Ewald use the same formula as LAMMPS.
2. When using the ChargeScreened class, the erfc_table_size argument is used. I don't quite understand the meaning of this argument. On what basis do you choose the size?
HWH> See https://pages.nist.gov/feasst/plugin/charge/doc/ChargeScreened.html . That argument can be set to 0 if you wish to disable tabulation of erfc. It is tabulated for optimization (as also done in a similar way in LAMMPS, etc).
3. For what purpose is RefPotential used in the tutorial for SPC/E? For what purpose is DCCB used here and why is the cut off value 0.9*3.165 (as I understand it is the value of the oxygen sigma parameter).
HWH> DCCB is briefly mentioned here: https://pages.nist.gov/feasst/plugin/monte_carlo/doc/TrialStage.html . The DCCB cutoff value is relatively arbitrary but can be optimized. In this case, DCCB essentially uses a hard sphere reference potential to quickly find cavities (with the help of a cell list, and considering only oxygen sites). I used 90% of the sigma value here because SPC/E h-bonds can lead to a surprising amount of overlap for the Oxygen LJ term.
4. Is the use of ChargeScreenedIntra necessary whenever we model molecules that have bonds and use the Ewald sum?
HWH> Yes, this term cancels the intramolecular O-H charge term computed in Fourier space.
5. Using dummy atoms containing only charge, is it necessary to apply any additional commands using Ewald sum?
HWH> A dummy atom requires you to modify the TrialGrowFile or perhaps simply use TrialTransfer. I happened to be working on adding an example TIP4P simulation and can follow up with you on this soon.
6. How do you choose the tolerance value when determining the last macrostate (for example, the maximum number of molecules)? In your paper (https://doi.org/10.1063/1.5123683) you used a value of 25.0, for example in J. R. Errington's paper (https://doi.org/10.1063/1.1572463) it was a value of 10.0. Is there any specific rule based on which to determine the value of this tolerance?
HWH> The variation you see in the last macrostate is mainly because we don't a priori know what N_max value is necessary. As an exercise, you can perform your analysis (say, VLE density calculation) with varying N_max. You may find that you can get away with less than 10, because probabilities e^-10 don't contribute very much in the GCE averages. We often go a bit higher in N than necessary just because we don't know when we start the simulation. Depending on your familiarity with the system and how quickly you want your results, sometimes its worth parallelizing into windows with a N_max that you know is too large and then simply throw out the largest density results that had prohibitively small acceptance and therefore did not converge well.
7. [...] I have not seen any examples of introducing any [rigid porous material/]framework in the system, so I assume for now it is not possible. Are you planning to implement such a possibility? Are there any particular difficulties in implementing such functionality?
HWH> We're actively working to support more functionality for these kinds of simulations. If the framework is rigid, it can be added as a second particle type. The framework particle can be added as the argument "add_particles_of_type[i]" when specifying the Configuration. You'll probably want to center the framework positions at the origin and ensure it fits properly in the Domain. With flat histogram,you can specify a particle_type in https://pages.nist.gov/feasst/plugin/flat_histogram/doc/MacrostateNumParticles.html because the macrostate is now the number of fluid particles, as well as in the Trials. You can also define groups in Config, and use the "group[_index]" commands in Movie via FileXYZ to print only the fluid positions.
HWH> We're actively working to support more functionality for these kinds of simulations. If the framework is rigid, it can be added as a second particle type. The framework particle can be added as the argument "add_particles_of_type[i]" when specifying the Configuration. You'll probably want to center the framework positions at the origin and ensure it fits properly in the Domain. With flat histogram,you can specify a particle_type in https://pages.nist.gov/feasst/plugin/flat_histogram/doc/MacrostateNumParticles.html because the macrostate is now the number of fluid particles, as well as in the Trials. You can also define groups in Config, and use the "group[_index]" commands in Movie via FileXYZ to print only the fluid positions.
Reply all
Reply to author
Forward
0 new messages