I was wondering if someone could help me understand some discrepancies I am getting between lammps and accepted water cluster energies using tip4p. Specifically, I am comparing energy values against values from the following site:
citation:
D. J. Wales and M. P. Hodges, Chem. Phys. Lett., 286, 65 (1998)
I have provided the data file for the 2 molecular cluster and input file I am using. I have noticed discrepancies between the two water cluster and twenty one water cluster. I have not tested any of the others yet.
When I run the provided files, I got a total potential energy of (using 4.1868 conversion factor from kcal/mol to j/mol)
lammps total potential energy = -5.7754 (kcal/mol)
lammps total potential energy = -24.18 kj/mol
literature value = -26.08757 kJ/mol.
if you take the 21 molecular cluster, and create a data file and use my input file, I get
lammps (kcal/mol) :: -217.1854
lammps (kj/mol) :: -909.311
literature value (kj/mol) = -916.07
Unfortunately, I cannot find anywhere in the paper where wales mentions the cutoff length used. He mentions he used the original values in the tip4p paper (jorgensen, revised tips for simulations of liquid water and aqueous solutions).
It does appear to be an issue with the cutoff, as you suspected, although it has me slightly confused. If you bump the cutoff values to 150.00 the total energy predicted is -26.10 kj/mol, which is close the literature value of -26.08 kj/mol
Could you explain to me this effect of the cut off here? I was under the impression that a cutoff of 9.5 for the vdw and 9.0 for the electrostatics was sufficient for energy convergence in tip4p water. At least I have seen papers mentioning the same value of 8.5 for an electrostatic cutoff.
Yes, energy conservation is unrelated to this issue.
If your cut-off is 9.5 Ångströms for Van der Waals, it can be that there are some water molecules more than 9.5 apart, and then you “miss out” on their energy. The potential values at the cut-off are typically so low that they barely influence the dynamics or structures you get, but they are missing, which can cause small differences in total energy when comparing with other results that did use a larger cut-off. In the cases I noticed it in my systems, it was roughly a few percent or so, that is why my hunch was that it could be cut-off related.
Think for example of a 2D hexagonal Lennard-Jones lattice. The contributions of the first neighbours to the potential energy are 3 epsilon per atom, but there is also contributions by the second neighbours (about 0.1 epsilon), etc., all the way to infinity. The potential decays as R^{-6}, so these contributions tend to be small (except for tricky things like pressures), but the number of interactions scales with R^3 (or R^2 in 2D), so the contributions to the total energy do not decay as quickly, which is why it sometimes makes a differences compared to a way larger cut-off.
I understand the effects of the cutoff for the vdw energy. If you play with the cut off in the input file i sent, you see the vdw energy never changes. That because all atoms are within the cutoff. The distances between atom pairs (atoms 4 and 8 are the m atoms) below. So its an energy convergence of the electrostatics issue. And because I thought this parameter set was optimized for an electrostatic cutoff of around 8.5 or 9.0, I cant remember, I would have thought the energy would be converging at the cutoff value. Perhaps I am incorrect.
atom number, atom number, distance(angstroms)
found this distance 1 5 2.767255
found this distance 1 6 1.894373
found this distance 1 7 3.274261
found this distance 1 8 2.848844
found this distance 2 5 3.186686
found this distance 2 6 2.350208
found this distance 2 7 3.734235
found this distance 2 8 3.253453
found this distance 3 5 3.186093
found this distance 3 6 2.349259
found this distance 3 7 3.733571
found this distance 3 8 3.252919
found this distance 4 5 2.661562
found this distance 4 6 1.791233
found this distance 4 7 3.145984
found this distance 4 8 2.748051
found this distance 5 1 2.767255
found this distance 5 2 3.186686
found this distance 5 3 3.186093
found this distance 5 4 2.661562
found this distance 6 1 1.894373
found this distance 6 2 2.350208
found this distance 6 3 2.349259
found this distance 6 4 1.791233
found this distance 7 1 3.274261
found this distance 7 2 3.734235
found this distance 7 3 3.733571
found this distance 7 4 3.145984
found this distance 8 1 2.848844
found this distance 8 2 3.253453
found this distance 8 3 3.252919
found this distance 8 4 2.748051
Hmm, in that case I am not really sure what is happening. From the documentation it seems like electrostatics for all atoms within the cut-off are computed directly, and all distances you have are also well within the coulombic cut-off. Furthermore, there is no shifting of potentials, so then changing the cut-offs should not matter at all, except if the damping term applied to the coulombic forces with coul/long is somehow cut-off related.
Thanks so much stefan and Yi. So the answer to the discrepancy is that the inclusion of the long range electrostatic solver no longer makes the water cluster in vacuum. The long range electrostatic solver will pick up interactions between the water cluster in the main box, and its periodic images. These periodic image calculations are not included in the literature values, as I believe they are done in vacuum. So counter to what we deal with 99% of the time, the direct cut off method is indeed the correct way of handling the electrostatics here.
The tip4p water model only includes vdw interactions between oxygens correct? Theres no vdw interactions at all associated with the m atom? The vdw interaction energy here looked funky to me, so I hand calculated the distance between the two oxygen atoms
The Wales group always provides very reliable data as far as my experience goes. BTW, sorry for scooping up here at the end with my unrelated comment but just wanted to drop these lines to encourage future people to further dig into prof Wales site when looking for bullet-proof research.
If you system is large, you can also use kspace_style MSM for fully non-periodic BCs. Though for a small cluster, the direct cutoff method is probably faster.
The tip4p water model only includes vdw interactions between oxygens
correct? Theres no vdw interactions at all associated with the m atom? The
vdw interaction energy here looked funky to me, so I hand calculated the
distance between the two oxygen atoms
