Coloumbic energy in GCMC calculation (CO2 in a box)

Dear all,

I am currently working on GCMC calculation of CO2 adosption in MOFs.

I started with a simple test with CO2 in a empty box and I played with several Coloumbic energy calculation method.

I tried three different setting for charge calculation.

(1)
pair_style lj/cut 12

(2)
pair_style lj/cut/coul/cut 12 12

(3)
pair_style lj/cut/coul/long 12.0
kspace_style ewald 1.0e-4
kspace_modify mesh 500 500 500 order 50

The first method did not include the electrostatic interactions, I expected the gas in the box would be slightly higher than this.21
The second method provide surprisingly high negative coloumbic energy and filled the box with hundreds of gas molecules. It means that the molecules are strongly attracting each other, right ?

The third method give high positive E_coul and high negative E_long, so it ends up with slightly negative value. The number of gas molecules is slightly less than the calculation without coloumbic interactions.

I am confused what is happening with the second and third method and not sure how to deal with that.

I have attached the files needed for the calculations.

I have also tried the three for CO2 in MOF. From my personal experience, the Coloumbic energy contributes to 5%-30% of the total uptake. However, the third method shows a lower numbers of uptake than the one without coloumbic interactions. And the second method shows a ten times higher uptake than the one without coloumbic interactions.

I am confused about the results. Could someone tell me what is the right way to describe coloumbic interactions.

Looking forward to your reply.

Many thanks in advance!

in.lammps (926 Bytes)

test.data (1.19 KB)

CO2.mol (397 Bytes)

I am glad to hear you are doing testing with pure gas first. I wish everyone did that before starting out using fix gcmc for the first time. Based on your description, it sounds like the code is working just fine. Different interaction models will give you different densities. If you want to know which model is best for CO2, read the literature.

Aidan

Dear Aidan,

Thank you very much for the prompt rely.

I am sorry I did not state my question clearly. I concerned about the coloumbic energy calculated in last two methods.

For the second method, I get results as below. The Coloumbic energy are surprisingly high.

Step nwater Temp PotEng TotEng Press Volume E_vdwl E_coul E_long E_pair E_bond E_angle E_dihed E_impro

0 1 300 452.44081 454.2293 -650.34957 27000 0 0 0 0 192.95904 259.48177 0 0

500 350 289.96154 -12653.509 -11746.837 1452.9113 27000 322.94723 -13428.897 0 -13105.95 192.95904 259.48177 0 0

1000 410 288.39993 -18056.818 -17000.289 3186.013 27000 524.98495 -19034.244 0 -18509.259 192.95904 259.48177 0 0

Similarly, in the third method where I include the ewald summation method. It shows

Step nwater Temp PotEng TotEng Press Volume E_vdwl E_coul E_long E_pair E_bond E_angle E_dihed E_impro

0 1 300 452.38191 454.17039 -650.17949 27000 0 11.564912 -11.623814 -0.058902223 192.95904 259.48177 0 0

WARNING: Using kspace solver on system with no charge (…/kspace.cpp:284)

500 156 303.8724 -194.77579 228.22623 392.05115 27000 -139.28908 1787.3222 -1842.8089 -194.77579 3.2728505e-25 1.6929547e-13 0 0

1000 262 298.37507 -515.08626 183.09227 509.69966 27000 -368.42725 2948.4402 -3095.0992 -515.08626 1.0408062e-24 4.4950865e-13 0 0

1500 279 315.36293 -611.91774 173.95318 502.91164 27000 -418.2962 3102.3809 -3296.0025 -611.91774 1.4593741e-24 6.24641

It shows that the molecules strongly repels at the real space and strongly attracting each other in reciprocal space.

Moreover, the density of the CO2 are overwhelming higher that expected. So I doubt there is something messed up with the coloumbic interactions. I guess there must be something wrong in calculating the Coloumbic energy. However, I checked the input file many times and I still do not know the reason.

Could you please give me some suggestions ?

Thank you very much for your time and patience.

Much appreciated,

Runhong Huang

Dear Aidan,

Thank you very much for the prompt rely.

I am sorry I did not state my question clearly. I concerned about the
coloumbic energy calculated in last two methods.

For the second method, I get results as below. The Coloumbic energy are
surprisingly high.

Step nwater Temp PotEng TotEng Press Volume E_vdwl E_coul E_long E_pair
E_bond E_angle E_dihed E_impro

       0 1 300 452.44081 454.2293 -650.34957
27000 0 0 0 0 192.95904
259.48177 0 0

     500 350 289.96154 -12653.509 -11746.837 1452.9113
27000 322.94723 -13428.897 0 -13105.95 192.95904
259.48177 0 0

    1000 410 288.39993 -18056.818 -17000.289 3186.013
27000 524.98495 -19034.244 0 -18509.259 192.95904
259.48177 0 0

Similarly, in the third method where I include the ewald summation method.
It shows

Step nwater Temp PotEng TotEng Press Volume E_vdwl E_coul E_long E_pair
E_bond E_angle E_dihed E_impro

       0 1 300 452.38191 454.17039 -650.17949
27000 0 11.564912 -11.623814 -0.058902223 192.95904
259.48177 0 0

WARNING: Using kspace solver on system with no charge (../kspace.cpp:284)

     500 156 303.8724 -194.77579 228.22623 392.05115
27000 -139.28908 1787.3222 -1842.8089 -194.77579 3.2728505e-25
1.6929547e-13 0 0

    1000 262 298.37507 -515.08626 183.09227 509.69966
27000 -368.42725 2948.4402 -3095.0992 -515.08626 1.0408062e-24
4.4950865e-13 0 0

    1500 279 315.36293 -611.91774 173.95318 502.91164
27000 -418.2962 3102.3809 -3296.0025 -611.91774 1.4593741e-24
6.24641

It shows that the molecules strongly repels at the real space and strongly
attracting each other in reciprocal space.

Moreover, the density of the CO2 are overwhelming higher that expected. So I
doubt there is something messed up with the coloumbic interactions. I guess
there must be something wrong in calculating the Coloumbic energy. However,
I checked the input file many times and I still do not know the reason.

this is an *extremely* unlikely hypothesis.

have you considered how long LAMMPS is around, especially in its
current C++ incarnation? don't you think, that people would have
noticed that something as fundamental as coulomb interactions would be
significantly wrong, specifically when using pair styles that are
commonly used? don't you think people will do some checks? so, the
gcmc process inside an MD code is a complex operation and thus small
oversights at that level are still possible, but questioning the
computation of coulomb interactions for standard pair styles is the
same as questioning at least half of the last 10+ year's worth of
publications.

please review aidan's comments more carefully and also follow his
advice to study the literature and possibly also some text book
explaining how ewald summation works.

axel.

Dear axel ,

Thank you very much for the email.

I am by no means questioning the accuracy and reliability of the lammps.

I am concerning that my settings is not correct so I get high Coloumbic energy. I am sorry I did not explain it well and it sounds like I am questioning the lammps.

I will try to read some literature, and try to get this done.

Thank you very much

Regards,

Runhong