Use of FixGCMC and FIX rigid

Dear Lammps Dev Team,

I would like to perform adsorption of CO2 gas on a porous structure followed by NVT/NPT dynamics. So my intention is to use fix GCMC+ fix rigid/npt/small…

I have the following question:

Given the following statement in the user guide :

“Because molecule insertion does not work in combination with fix rigid, simultaneous use of fix rigid or fix rigid/small with this fix is not allowed.” – i.e. use of Fix GCMC and fix rigid or fix rigid/small together.

Is it only restricted to use of fix GCMC and fix rigid (or) fix rigid/small ??? Or for all variants of fix rigid i.e. fix rigid/nvt fix rigid/npt etc.

If yes fix poems is the alternative?

Cheers,

Ravi.

Dear Lammps Dev Team,

I would like to perform adsorption of CO2 gas on a porous structure followed
by NVT/NPT dynamics. So my intention is to use fix GCMC+ fix
rigid/npt/small…

I have the following question:

Given the following statement in the user guide :

“Because molecule insertion does not work in combination with fix rigid,
simultaneous use of fix rigid or fix rigid/small with this fix is not
allowed.” – i.e. use of Fix GCMC and fix rigid or fix rigid/small together.

all variants.

Is it only restricted to use of fix GCMC and fix rigid (or) fix
rigid/small ??? Or for all variants of fix rigid i.e. fix rigid/nvt fix
rigid/npt etc.

If yes fix poems is the alternative?

no.

axel.

Thanks for your reply...

Is there any alternative (i.e. CO2 as rigid molecule ) ... as SHAKE or RATTLE won't be able to do so..

Regards,
Ravi.

Thanks for your reply...

Is there any alternative (i.e. CO2 as rigid molecule ) ... as SHAKE or RATTLE won't be able to do so..

not without some serious C++ programming. e.g. theoretically linear
rigid molecules could be handled by SHAKE/RATTLE, as they can be
projected on a diatomic, propagated accordingly, and then projected
back.

axel.

fix shake is the only rigid constraint that currently works with fix gcmc. If fix shake can not handle a linear molecule, your only other solution is to use a flexible CO2 model.

Aidan

Thank you Axel and Aidan.

I have one more observation to share on the same problem i.e. adsorption of CO2 on a porous structure (no charge)

Please note that I ran only GCMC without dynamics:

a**. Scenario 1 :** Simulation wit out electrostatic interactions

GCMC_Charge.profile (9.49 KB)

GCMC_nocharge.profile (8.47 KB)

Your question is not very clear, but I think you are asking why you are getting different results when using lj/cut/coul/cut versus lj/cut. It’s hard to say without spending a lot of time examine both cases in detail (which you can always do yourself). In general, when you change the pair style, you expect results to change, unless you know for sure that they should not, which is clearly not the case here. The CO2 molecules will interact with each other differently in both cases, unless you are using a molecule template co2.txt that contains no charges. Another possibility is that the full_energy option is being turned on (check the log file) and the energy of a single molecule is different in both cases, due to intramolecular charge-charge interactions. You can verify this by inserting a single molecule in an otherwise empty box and looking at the energy in both cases.

Hi Aidan,

Thanks for your reply. My comments are as follows:

MY co2 molecule has the charge and I expect the change in behaviour with different pair styles (i.e. with and w/o columbic interactions).

But my concern is - in one case able to reach equilibrium in finite steps and in another case it adsorbs as many as you proceed with #of steps.

Full energy option is not turned on as I am not worried about long range interactions at this point (coul/cut).

Regards,

Ravi

Aidan,

Also, I strongly believe the system is getting over populated with the gas molecules when the charges are turned on. To confirm that I ran NPT simulation after 1000000 GCMC steps and observed the increase in volume of the simulation cell.

Regards,

Ravi.

Your description of the problem is still unclear, despite my pointing this out to you. This suggests you are more interesting in complaining than learning.

Aidan,

Also, I strongly believe the system is getting over populated with the gas molecules when the charges are turned on. To confirm that I ran NPT simulation after 1000000 GCMC steps and observed the increase in volume of the simulation cell.

Regards,

Ravi.

Aidan,

Also, I strongly believe the system is getting over populated with the gas molecules when the charges are turned on. To confirm that I ran NPT simulation after 1000000 GCMC steps and observed the increase in volume of the simulation cell.

Regards,

Ravi.

Hi Aidan,

Thank for your good suggestion!!. Earlier, my concern was all about my set up and interaction parameters between gas and solid . But that’s not the case …

My complaint is as follows :

I performed GCMC simulations on an empty box by varying the pressure of the gas (CO2) reservoir using LAMMPS (21 Mar 2016) version. EPM2 model was used for CO2 gas. I ran the simulations with and without considering the charge. I observed that GCMC results from LAMMPS are not correct if we consider Electrostatic Interactions. All the output results are attached in the excel. All the simulations were performed at fixed temp (T=350K)

LAMMPS-Complaint.xlsx (58.5 KB)

“I observed that GCMC results from LAMMPS are not correct if we consider Electrostatic Interactions.”

I took a closer look at this. There is nothing strange about this behavior. With short range electrostatic interactions, the CO2 molecules are able to form strong attractive interactions between C and O atoms on neighboring molecules. Hence, for the same chemical potential, you go from a vapor (no electrostatics) to a liquid (with electrostatics). In other words, a different molecular model gives a different phase. This is basically what I predicted was happening back on June 17:

“It’s hard to say without spending a lot of time examine both cases in detail (which you can always do yourself). In general, when you change the pair style, you expect results to change, unless you know for sure that they should not, which is clearly not the case here.”

Interestingly, I obtained a density of 1.7 g/cc which is close to the density of dry ice. Incorrect result, certainly. A bug in LAMMPS, absolutely not.

Aidan

Thanks for your remarks.

When I am performing all the simulations well above the critical temp (304.25 K), how it can predict the liquid phase except the model parameters are wrong which I don’t think so.

Dry ice can form at atmospheric pressure below ~ 200K and 5 atm which is not the case here.

Could you please explain why do you think it’s not a bug in lammps??

Regards,

Ravi.

Could you please explain why do you think it’s not a bug in lammps??

Yes I can.

A good general rule of thumb for any kind of computation is: If it has not
been tested, it is almost certainly wrong.

This rule is closely related but different from Axel's law of GIGO. If
fact, were GIGO not true, then testing would not be very effective.

A version of this rule more specific to your problem is: "Model parameter
values input to a code probably contain errors, unless those inputs have
been used with the same code to reproduce a non-trivial result"

This is not a criticism of your cognitive abilities, but rather is due to
the entropic nature of complex calculations. Assume computer programs and
inputs are generated in a manner that allows for the occurrence of random
errors. There are exponentially many ways to take a correct program or
input make it incorrect. Many, maybe even most, of these errors will be
detected by carefully scrutiny. However, some fraction of errors will go
undetected, so the chances of hitting upon the actual correct input remains
exponentially small, even after careful scrutiny.

The only way to overcome entropy is through testing, which eliminates large
chunks of error space, thanks to GIGO. Since LAMMPS has been subjected to a
lot of testing, and your inputs have not, it is natural to assume that any
errors are in your inputs, not in the code, all other things being equal.
So, the pertinent questions are:

1. Have you managed to reproduce any non-trivial result with LAMMPS using
these model parameters?
2. Is there any evidence that points specifically to an error in LAMMPS?

Aidan

Could you please explain why do you think it’s not a bug in lammps??

Yes I can.

A good general rule of thumb for any kind of computation is: If it has not
been tested, it is almost certainly wrong.

This rule is closely related but different from Axel's law of GIGO. If fact,

FWIW, it is not "my" law. GIGO existed long before i knew how to use a
computer. i've learned it from the "jargon file", which collects lots
of "hacker folklore" often going back to the times when men were real
men, computers were real computers and hackers were a very select
elite of capable programmers that looked down on mere users of
software and computers... oh, wait... :wink:

http://www.catb.org/jargon/html/G/GIGO.html

Yes, but you have certainly helped raise awareness of GIGO on this mailing list. As Wikipedia rightly points out, GIGO was already known to Charles Babbage:

On two occasions I have been asked, “Pray, Mr. Babbage, if you put into the machine wrong figures, will the right answers come out?” … I am not able rightly to apprehend the kind of confusion of ideas that could provoke such a question.

Aidan

Hello,

Yes !! With the same set of parameters, I ran the NPT simulations for given N and Pressure and Temperature (350K) and computed the density. Please see the results below.

So you were able to approximate the ideal gas law using NPT MD in LAMMPS? That is not a non-trivial result, it is a generic behavior of molecular models. So it does not prove that you model is correct. You will have to do better than that.

Regarding you comparison between NPT MD and GCMC, where you specify the same input pressure in both cases, you seem to have the expectation that the GCMC system will equilibrate to an ideal gas at the specified pressure. That is not what the documentation says, and it is not true in general. It is only true in the limit of sufficiently low pressure so that ideal gas law is a good approximation. What the pressure keyword provides is an alternate way to specify chemical potential as the pressure in a fictitious ideal gas. This equivalence is encoded in the following two lines of code:

zz = exp(betachemical_potential)/(pow(lambda,3.0));
if (pressure_flag) zz = pressure
fugacity_coeff*beta/force->nktv2p;

I will also mention that both your charged and uncharged models have a serious pathology. The pair energy is truncated discontinuously to zero at a distance of 8 angstroms. The original molecular model was almost certainly not based on such a short cut off. It may also have used an energy shift or a tail correction. If you are serious about reproducing CO2 properties, you will spend the time to learn about these issues by reading the literature and learn how to handle them in LAMMPS by reading the documentation and performing non-trivial tests.

Asking the mailing list to do your work for you is not a viable path to success.

Aidan