[lammps-users] GCMC/NVT with fix rigid command

Dear all,

I want to simulate density of methane molecules in bulk condition. For this purpose, I use GCMC/NVT simulation with the CVFF force field. I use Peng-Robinson EOS to calculate fugacity coefficient. The version of Lammps is 8 Apr 2021.

When I use fix NVT/rigid/small and GCMC, I face some errors. On the one hand, the GCMC/NVT simulation runs successfully without using fix rigid command. On the other hand, without using the fix rigid command in GCMC/NVT simulation, the output pressure is more value than the input pressure (desire pressure), and the density is very different from the NIST density.

For more investigation, I use fix rigid command. I run the code with a molecule command in NVT and NVE ensemble separately, and the simulation runs successfully.

When I use molecule commend for methane in GCMC simulation with fix rigid command, should I use the exact coordinate methane atom similar to bond and angle parameters in the force field? Because when I use the exact parameters, I get this error at the first step. (ERROR: Lost atoms: original 430 current 410 (src/thermo.cpp:425)).

In contrast, when the xyz coordinate is not as the same as the bond and angle parameters in the force field, the simulation continues for some steps, and then shows the same error. (ERROR on proc 0: Non-numeric positions - simulation unstable (src/nbin.cpp:125)).

bond_style harmonic

angle_style harmonic

bond_coeff 1 340.6 1.087

angle_coeff 1 39.5000 109.5000

1- When I use exact xyz coordinate for methane in molecule file (e.g., bond length: 1.087 A and angle=109.47) :

Coords

1 1.645 6.963 2.949

2 1.336 6.305 2.141

3 1.167 7.932 2.830

4 1.353 6.529 3.902

5 2.725 7.086 2.923

Step Temp Press Density PotEng KinEng v_iacc v_dacc v_tacc v_nC

0 305.8817 -688.79444 0.026640159 -0.69389448 115.79546 0 0 0 64

ERROR: Lost atoms: original 430 current 410 (src/thermo.cpp:425)

2- When I don’t use exact xyz coordinate for methane in the molecule file (e.g., bond length: 1.069-1.071 A and angle=109.4-109.5) ::

Coords

1 -0.370 0.900 0.000

2 0.700 0.900 0.000

3 -0.727 0.122 0.643

4 -0.727 0.731 -0.995

5 -0.727 1.845 0.351

Step Temp Press Density PotEng KinEng v_iacc v_dacc v_tacc v_nC

0 307.15115 -244.05364 0.026640159 -0.63825465 116.27602 0 0 0 64

25 210.52469 27.827673 0.036630218 -6.3652709 109.81847 0.79173838 0.52256532 0 88

50 222.49342 15.324639 0.036630218 -6.4317018 116.06186 0.79173838 0.52256532 0 88

75 240.52761 7.1352222 0.036630218 -7.4705819 125.46924 0.79173838 0.52256532 0 88

100 256.60519 16.567044 0.036630218 -7.379648 133.85598 0.79173838 0.52256532 0 88

125 270.21461 30.144297 0.036630218 -7.4233475 140.95523 0.79173838 0.52256532 0 88

150 280.66357 44.986566 0.036630218 -6.6454443 146.40584 0.79173838 0.52256532 0 88

175 289.98404 53.31387 0.036630218 -6.1965964 151.26778 0.79173838 0.52256532 0 88

200 299.3352 51.874975 0.036630218 -7.0024242 156.14574 0.79173838 0.52256532 0 88

225 376.38567 165.58168 0.037878976 -8.105518 203.07009 0.76117983 0.5573081 0 91

250 362.86073 192.8589 0.037878976 -9.3451126 195.77302 0.76117983 0.5573081 0 91.

7925 348.66461 120.08163 0.092408051 -103.27467 460.41121 0.59688882 0.5246215 0 222

7950 335.56181 118.66907 0.092408051 -103.84359 443.10899 0.59688882 0.5246215 0 222

7975 325.72888 108.39609 0.092408051 -105.97782 430.12461 0.59688882 0.5246215 0 222

8000 317.74374 101.8956 0.092408051 -107.85967 419.58025 0.59688882 0.5246215 0 222

8025 13147.427 5482.1281 0.094905566 -96.65286 17831.443 0.59783401 0.52275734 0 228

8050 26349.744 42803.876 0.094905566 3048.3808 35737.329 0.59783401 0.52275734 0 228

8075 nan nan 0.094905566 nan nan 0.59783401 0.52275734 0 228

8100 nan nan 0.094905566 nan nan 0.59783401 0.52275734 0 228

8125 nan nan 0.094905566 nan nan 0.59783401 0.52275734 0 228

8150 nan nan 0.094905566 nan nan 0.59783401 0.52275734 0 228

8175 nan nan 0.094905566 nan nan 0.59783401 0.52275734 0 228

8200 nan nan 0.094905566 nan nan 0.59783401 0.52275734 0 228

ERROR on proc 0: Non-numeric positions - simulation unstable (src/nbin.cpp:125)

I know that many simulations used methane all-atom with fix rigid comment and GCMC simulation in the articles with LAMMPS software, and I check the input file several times, ,but why is it impossible for me to do this? Is there any mistakes in the input file? I will be happy if anyone helps me.

I also tried different tfact and timesteps and got the error again. I also upload the LAMMPS input files on this mail.

Best regards

Saeed

log.lammps (47.1 KB)

CH4.txt (580 Bytes)

in.gcmc.ch4 (2.67 KB)

Saeed,

I am not much of an expert in running GCMC calculations, but I see multiple issues with your input deck where you are in conflict with very basic MD simulation best practices and common knowledge. I don’t think that any results from such a simulation setup you have can be trusted at all. Since GCMC is an advanced simulation technique, you should get the basics working properly first and then attempt to add GCMC to it later and carefully observe what are good settings there as well. Please note that MD simulations are as much a craft as they are science, so there are things where you just have to practice, experiment, and observe and check what you are doing in small steps thus and incrementally approach such a complex simulation setup in a fashion that will produce trustworthy results.

What is evident and very concerning is how much your temperature is fluctuating and “exploding” after during your run. That indicates problems long before it actually crashes. The crash is likely an indirect consequence and not something that has a direct cause.

  • first, your system is tiny which means that temperature will fluctuate a lot, you will have large(r) finite size effect errors and generally difficulties to get converged results except for the most simple and trivial systems. Just increasing the system after reading in the initial data with replicate 2 2 2 to 8x the size makes things much smoother. And that is still going to be a small (initial) system with only about 2500 atoms.

  • second, there is no equilibration. you start from an initial structure with molecules on regular lattice points, but random orientations. so you should first equilibrate that initial system to the target temperature and target conditions for starting the production run. that may take some time and that may include running for a while at elevated temperature to remove the “memory” of the regular initial configuration.

  • third, the choice of timestep seems problematic. for a system with hydrogen atoms and flexible bonds, a timestep between 0.25fs - 0.5fs would be the conservative to aggressive range. this can be increased to 1fs-2fs when using fix shake, but not with fix rigid since the time integration of the rotational degrees of freedom does not allow such a large timestep. you will get bad energy conservation with 1fs up to the point where the thermostat cannot hide it. the best practice in this case is to first run a simple system (without fix gcmc and without a thermostat) after equilibration to monitor energy conservation and adjust the timestep. keep in mind that you plan to do a “disruptive” operation afterwards, so it is probably a good idea to have a more conservative choice there, so that the disruptions will not kick the system off-balance.

  • fourth, since any GCMC step is a disruption, you should leave sufficient time between GCMC attempts so that there is sufficient time of relaxation after an insertion or deletion event and that the system mostly stays in equilibrium. be wary of how much using the thermostat can hide bad choices for GCMC settings. also the overall impact and amount of disruption relative to the whole system will (obviously) be much less when you simulate a larger system. notice how your temperature always jumps after an GCMC step. it probably would also be advisable to not only make the GCMC attempts less frequent but also make fewer attempts as well.

  • bonus items. your neighbor list skin seems very large. your performance should improve with a smaller value. also, for a larger system, the performance of PPPM is usually superior to direct Ewald summation (O(N log N) versus O(N^(3/2)).

since this is all about basic MD simulation technique and common sense, I suggest you also discuss this with your adviser/supervisor and/or senior/experienced colleagues to take advantage of their experience and knowledge. you won’t be getting that kind of support and advice from the mailing list (that is not its purpose). people usually will explain basic stuff to you a limited number of times until they tire of having to explain things that people trying to use advanced MD simulation techniques should already know well.

Axel.

1 Like