Unexpected results with adding carbon atoms. Is there a failure of describing potential?

Dear LAMMPS users/developpers,

I confined water molecules between two parallel graphene sheets and hereby tried to calculate the pressure tensor coefficients and viscosity as well (files attached for more details).
In the case of no carbon sheets, everything works well. However, when I add graphene sheets, unrealistic results for viscosity starts to appear. The viscosity should be 1 Pa.s for this case. But I get 2.30E7 Pa.s which is too strange. I personally believe something wrong with calculating of pressure coefficient in the presence of carbon atoms.

So, I have tried:

  1. I ignored C-C interactions:
    neigh_modify exclude type 3 3
    but the viscosity result was too low (0.0007).

  2. Since the carbon atoms were frozen here, there are no inconsistency should be arisen between tersoff and Lj potential. However, I get 357 Pa.s for using tersoff potential (too different compares to Lj potential).

  3. when I exclude C-C interactions both in tersoff and Lj potential, the results look similar.

Is there something wrong with using potential and hereby need to use some extended potentials (like Reaxff) here?

Thanks for your comments,

solvate.in (3.0 KB)
solvate.data (128.2 KB)

It seems that what you are measuring is the viscosity of the water + cnt system. Shouln’t you use only the pressure coefficients for the water alone instead of the pressure coefficients for the entire system?

Btw I know that this method is valid for calculating the viscosity of bulk fluids, but are you sure that it is the right method to measure the viscosity of confined fluids?



Nanoconfined water is a very different substance from bulk water (for example, its dielectric constant is about 40 times smaller). You should expect neither methods nor material properties to carry over from bulk to nanoconfinement simulations, in general.

1 Like

Additionally, assuming you are interested in shear viscosity of the system, you can simply push one graphene wall in one direction and the other graphene wall the opposite way and see what velocity gradient is developed in the water. That way you have a concrete (in silico) observation against which you can benchmark other techniques.

1 Like

Thanks all for fruitful comments.

Motivated with @simongravelle comment regarding the pressure coefficient, I fixed this fault and now attached the modified input file here. But again the output result for viscosity is 0.001 Pa.s which is too low compared to reported value by other literature (0.9 Pa.s).

By this way, for the second attempt, I try to use reaxff potential here. So, I added below commands (input file also attached):

pair_style	reaxff  NULL
pair_coeff      * * ffield.reax.FC H O C
fix      re all qeq/reaxff 1 0.0 10.0 1e-6 reaxff

and unfortunately it was failed to converge. So, any comments concerning the convergence and using of reaxff potential will be highly appreciated.

WARNING: Fix qeq/reaxff CG convergence failed after 200 iterations at step 11 (src/REAXFF/fix_qeq_reaxff.cpp:769)

solvate.in (3.1 KB)
solvate.data (128.3 KB)

Wait, bulk water viscosity is 0.001 Pa.s (at ambient temperature and pressure), what literature reports a viscosity of 0.9 Pa.s for water? In which system and in what conditions ?

If your simulation fails with a simple water model, don’t use reaxff, your are making things way too complicated and more likely to fail even more. For instance reaxff requires a very short timestep, like 0.1fs, but you are using 1fs.


1 Like

For example, the following literature reports the viscosity of water confined between two graphene sheets separated with h=7.5 A gap is 0.9 Pa.s (figure 1). However, I got 0.001 Pa.s with spc/e water model for the similar structure (input file’s been attached with my previous comment).


Oh yes you are right. In this case, to validate your protocol, I would start by trying to reproduce their results in the “easy” limit, i.e. for large slit pore with h > 2 nm. In this limit it is very clear that you should get something like 1 mPa.s (the value will depend a little bit on the water model, but is should still be somewhere between 0.3 mPa.s for the less good water models like tip3p and 1mPa.s for the best water model like tip4p/2005).

1 Like