I performed a regular NPT MD simulation for 864 water molecules at 300K. The water model is SPC/E.
When I looked at E_coul in log file, it shows a positive value, which is strange to me.
Another question is the potential energy in the simulation does not show to be equal to the summation of E_bond, E_angle, E_dihed, E_impro, E_vdwl, and E_coul.
The following is the log file from the simulation. Could you help me to understand this?
The Coul eng you are thinking of is
the sum of Ecoul and Elong. The latter
is also included in the total PE. You should
also make sure that the special_bonds settings
are what they should be for water.
Thank you for your prompt response! After I take the Elong into account, everything looks good now.
For SPC/E model, I didn’t see anyone using any special bond setting.
It seems all the properties of water calculated from the simulation are good so far.
Dear Leo,
Generally, special bond setting is not used for water. as 1-2 and 1-3 interactions are the only ones which are present in water (only 3 atoms per molecule) which are in 0 0 in almost all cases. The 1-4 interaction is not applicable to water.
Thank you very much for your comments. I noticed for most forcefields such as amber or charmm, the 1-2 and 1-3 interaction weighting factor are set to zero. 1-4 interaction depends on forcefield. I am about to use cvff force field. However, I was not able to find the weighting factor so far. Is this available in LAMMPS documents? If I can’t find it, is it safe to assume it is zero?
Dear Leo,
Run the NPT simulation with PPPM with high tolerance (1e-6) and see what density of water you get with your model and then match with those other’s have predicted with same force field. Radial distribution function is also good to compare with published literature. I am sure you can find a lot of papers with SPC/E model.
With my current script for water runs good, I can get properties such as g® and diffusion similar to those in literature. Eventually, I will use some small organic molecules in water. I will keep looking for the weight factor of 1-4 interactions for cvff forcefield.