[lammps-users] (MD / LAMMPS Beginner Question) Water Molecules

I would like to run a simulation of a few hundred water molecules in a simulation box. First with the TIP4P model.

I’m reading the lammps section on the water models, and my question is, in order to get the water molecule models that i want, i need to create two types of atoms in my simulation (Hydrogen and Oxygen), and use the bond, angle, and other parameters between the molecular atoms given by lammps for whichever water model I want?

Do I also have the choice of doing this in the input script, or a data file?

forgive the first long drawn out question.

very respectfully,

LAMMPS won't make the configuration for you. I use a program called packmol. It's free. You can find it with google. There is an online howto for tip4p - http://lammps.sandia.gov/doc/Section_howto.html#4_8 . Be sure to use the right parameters for the electrostatics method you use. The document lists both.
Cheers,
Matt

Quoting Daniel Casimir <[email protected]>:

In particular, there is one other issue specific to the TIP4P models. Because of the way the algorithm for finding the location of the point charges is written, you must specify the water molecules in the order O–H1–H2, while the angles have to be specified with the atoms in the order H1–O--H2. If you get the first wrong, the system will give you an error message pretty quickly stating that an “H” atom couldn’t be found; if you don’t do the second, LAMMPS will not give an error message, but the mistake will cause the simulation to crash fairly quickly.

–AEI

Is packmol compatible with lammps?

Specifically, I’m asking if you use packmol to create molecules, then run your simulation in lammps?

Or, do you just use packmol as an alternative md simulator to do everything including the actual running of the simulation?

Packmol makes a simple xyz type coordinate file from a molecule template and an input file that specifies the number of molecules and the box size - and that is all it does. I then take that xyz file and turn it into a LAMMPS format configuration input. You have to add all the LAMMPS header information, box boundaries, atomic masses, and so forth. You also have to add the atom id, molecule id, charge, and so forth to the positions - as well as the bonding and angle information to the end of the file. I usually end up writing a simple fortran code to do all that. Not exactly pretty, but it works.
Matt

Quoting Daniel Casimir <[email protected]>:

Matt,
This is intriguing me. I need to solvate a large molecule in a water box. Packmol will do this? I just need to provide it with one of each type of molecule? Basically my concerns are:

  1. If my large molecule is in some irregular shape, Packmol will know to avoid the excluded volume when placing water molecules?
  2. Are the HOH coordinates given in order in the xyz output, so bond and angle information should be straight forward? In other words, I don’t have to search through positions to see which atoms are potentially bonded to which atoms, right?

Thanks!

Jonathan

Packmol makes a simple xyz type coordinate file from a molecule
template and an input file that specifies the number of molecules and
the box size - and that is all it does. I then take that xyz file and
turn it into a LAMMPS format configuration input. You have to add all
the LAMMPS header information, box boundaries, atomic masses, and so
forth. You also have to add the atom id, molecule id, charge, and so
forth to the positions - as well as the bonding and angle information
to the end of the file. I usually end up writing a simple fortran
code to do all that. Not exactly pretty, but it works.

if you want to do it in a more elegant way, i suggest to have a look
at VMD scripting and particularly the topotools plugin for VMD. this
plugin can derive angle and dihedral definition for many simple
force field cases directly from the bonding pattern derived from
a distance based guess. you can determine unique bond, angle,
dihedral types and have it write out a LAMMPS compatible data
file. assignment of potential parameters will (currently) have to be
done manually, but i have some script fragments that can help
there as well.

cheers,
    axel.

Yes to 1 and 2. There is some functionality for setting a precision to avoid overlap. I usually use a low precision (which makes the code run faster) and just do a quick minimization before an equilibration run. Packmol writes out the molecular atomic coordinates in the same order they occur in the input template. There are likely better codes for solvating large molecules, though. But I don't do bio stuff so I'm not really aware of any.
Matt

Quoting Jonathan Lee <[email protected]>: