Dear LAMMPS forum,
Appreciate your assistance.
I intended to calculate water crystal taking the source materialsproject hexagonal mp-703459 geometry with a,b= 7.50 Å,c= 7.06 Å; α,β=90.00 º,ɣ=120.00 º. RCO=0.995A, Angle HOH=107.8 deg.
Thus I built a water LAMMPS input crystal that contains 94770 atoms (31590 H2O molecules) using Biovia, VMD, and some more tools.
As a first approach I’ve tried to make simulations based on the official LAMMPS package gle examples (TIP4P/2005f, input files in.h2o-quantum and in.h2o-smart are attached). Another sample (for TIP4P/2005 potential) was taken at LAMMPS course (LAMMPS Online ). I also examined simulations using Second-Generation ReaxFF for water (J. Phys. Chem. B 2017, 121, 6021−6032).
There are some questions (and I would apologize if some of you find these trivial) that I haven’t found an answer to in the LAMMPS manual, the internet, or papers.
Both TIP4P/2005 and TIP4P/2005f contain O-M distance parameter, where M is compensating negative charge. The question is whether it is must to add M center into the atom’s data file or LAMMPS “knows” to find its position upon the corresponding parameter “qdist”. It is because all (TIP4P/2005 and TIP4P/2005f) examined example files contained H and O centers only with a negative charge located on O atom. It also could be a “headache” to build a big OHM data file.
Another question concerns the bond_coeff and angle_coeff parameters.
For TIP4P/2005 these are bond_coeff 1 0.0 0.9572 and angle_coeff 1 0.0 104.52, where 1 – is O, 0.9572 is Req(OH), and 104.52 is HOH eq. angle. Could you please explain what is 0.0?
For TIP4P/2005f these are bond_coeff 1 0.9419 607.19354 -1388.6516 1852.577
angle_coeff 1 43.93 107.4. Again, 0.9419 is Req(OH), 107.4 is HOH eq. angle. What bond_coeff 607.19354 -1388.6516 1852.577 and angle_coeff 43.93 parameters mean?
Another question concerns the results of the water crystal thermalization results that were calculated by me using all, TIP4P/2005, TIP4P/2005f, and ReaxFF. I did try various schemes, were the best (in my opinion) gives the scheme that includes: Energy minimization (minimize 1.0e-5 1.0e-6 1000 100000), npt (fix H2Onpt water npt temp 08.0 07.0 20.0 iso 0.05 0.01 300.0), and then nvt (fix H2Onvt water nvt temp 007.0 005.0 200.0).
The resulting slab looks like it retains a hexagonal structure except for the outer edges that look quite distorted (shown on the attached picture) that is probably due to uncompensated partial H and O charges. Does anyone met that effect, or it’s could be an artificial one? Any comments would be very appreciated.
No need to add the M point to the data file, LAMMPS re-calculate the position of the M point when needed.
It is the energy, put to 0 for rigid models.
Don’t you use periodic boundary conditions in a way that the crystal is effectively infinite? (if not, it is strange that you involve npt during equilibration)
This can be looked up in the documentation for the corresponding bond style and angle style.
LAMMPS will compute the forces of the system (i.e. atom positions, bonding topology, force field, box dimensions, boundary conditions, simulations settings) and then act on them as instructed. If this is not the system you aim to model, you will have to adjust your input accordingly. However, this is not really a LAMMPS problem, but rather a problem of the applied physics or the model.
Thank you very much for a quite complete explanations.
I will try to better understand the reason for using the bond_coeff and angle_coeff parameters. In some examples it is explained that for a rigid model (such as TIP4P/2005) these parameters are not important, and in some examples they are omitted at all.
Yes, I did use boundary p p p. I supposed that since materialsproject H2O geometry is slightly different from the TIP4P equilibrium one, npt could assist to equilibrate the system. Indeed, using nvt only results in more structural distortions.
Do you mean that insufficient box size can affect the edges of the system? The box size was 500x500x500A when the crystal one - more or less 100x100x100.
If this is a crystal, the box dimensions have to perfectly match the molecular structure and provide a perfect continuation. Even if the molecules are too close or too far from each other, the system must remain homogeneous for symmetry reasons, if set up properly.
Running MD at such low tempertures is rather pointless. Your results should not differ much from what you get with plain minimization or minmization with fix box/relax.
If I may ask one more question, I do not see the way to implement the TIP4P/2005f Dr and beta Morse interaction parameters and also the angle bending strength constant K(theta)
No idea what you mean by that. I suppose flexible water model are straightforward to implement, remove the shake constraint and add the right bond and angular constraints with energy different from 0.
For the functional forms, you can always look at the documentation of the corresponding bond and angle styles. Different research groups have different naming conventions in the potential function, but you can figure those out easily. For Lennard-Jones and Morse there also are different ways to present the functional form, but you can always convert the corresponding parameters from one to the other. It is just a simple exercise in calculus.
