and when I use the stress/atom command, the result looks good.
However, I have to use centroid/stress/atom command in my modified structure and the bad thing, it doesn’t support kspace_style (my lammps version 29 Sep 2021). So I have to use lj/cut/tip4p/cut as follows:
pair_style lj/cut/tip4p/cut 2 1 1 1 0.1546 12.0
But there the result looks twice the expected value.
Is there any pair_modify options that more or less include long ranged Coulombic interactions?
This does not make any sense. Please look up in a text book or some other credible resource how Ewald summation works. Then you should realize that the real space calculation (which is what pair style does) is different between the two models for computing Coulomb, hence the different pair styles. Consequently, you cannot simply “fix up” a ‘coul/cut’ Coulomb after the fact. It is a different computation.
I have already commented on the futility of what you are trying to do here before, so no need to repeat it.
Why you have to use centroid/stress/atom instead of stress/atom? If that’s for calculating heat flux (as why it’s introduced in the first place), please note that virtual sites need special treatment (see, for example, J. Chem. Phys. 153, 194105 (2020), eq. (12)) and I don’t think the code in LAMMPS has considered this. It may be still possible to continue (use an explicit 4-site model, modify the source code to enable kspace centroid stress (the kspace part is same for stress/atom and centroid/stress/atom), calculate the extra correction term by yourself, etc.), but I strongly encourage you to discuss with someone more familiar with this topic first, and compare your results with some reference to make sure you have implemented it correctly.
Thanks for your fruitful comments. I appreciate it.
I use centroid/stress/atom to calculate viscosity of water inside the carbon nanotubes by using the Green-Kubo formula. I guess the water forms anisotropic media inside the nanotube, so I decided to use centroid/stress/atom instead of stress/atom.
Yes, they report a correction term for the virtual atom since it cannot be treated as a real atom as there is no independent momentum associated with it. And, actually Lammps does not include this term. So, for anisotropic media in which the normal stress is not valid, one needs to write a post-processing code to compute the extra term as well as modify the source code to enable kspace centroid stress.
All in all, since the rigid model freezes the high-frequency intramolecular vibrations and also lj/cut/tip4p/cut is a pair potential with spherical symmetry, it would be a good approximation to rely on the normal stress (stress/atom).
Unluckily I don’t know much about viscosity calculation. Here are some general comments I can say:
make sure that you really need centroid stress (e.g. has it been pointed out in the literature?), since you’d save a lot of effort if the plain stress/atom is enough;
compare results between different choices. For example, keep using stress/atom and compare between tip4p/cut and tip4p/long, or keep using tip4p/cut and compare between stress/atom and centroid/stress/atom; it’s possible that using a cut style introduces more error than using the plain stress/atom. You may also try different cutoff for tip4p/cut, 3-atom vs. 4-atom model, etc.
you may want to search literatures to see if a virtual site correction is required/available for your use case, or try to derive it by yourself. For typical tip4p models the correction is usually small on thermo conductivity, but I don’t know about the viscosity.
How large is your system and how many waters do you have in it?
Especially, if there is no water around it, you can just use a very, very, VERY large Coulomb cutoff and you will have a pretty accurate Coulomb. (especially, if you would use one of the smoothly force shifted styles). Using Kspace for a 1d-periodic system is not so accurate either, if there is vacuum around.
You may need to boost the neighbor list one and page parameters, but since your system will grow with the cutoff effectively only in one direction instead of 3, the cost of a large cutoff is going to be much less prohibitive than with a 3d-bulk system. LAMMPS should be able to handle a cutoff that is larger than the size of the box…
Why not just use a three site model like SPC/E or OPC3? Whatever gain in accuracy you get from a virtual site in TIP4P is going to be swamped by (1) the additional complexity in post-processing due to the virtual site (as you’re finding out yourself) (2) the error from simulating nanoconfined water with bulk-parametrized force fields (3) the loss of intramolecular motions with rigid water (4) uncertainty in parametrizing the carbon-water interactions, including (5) induced charges on the carbons from water due to the carbon nanotubes being conductive.
Given how far a TIP4P model is from being “realistic”, it seems sensible to either go “down” to three sites for ease or “up” to ab initio molecular dynamics for realism. You’re getting the worst of both worlds staying where you are, especially if you are relatively new to molecular dynamics.