Following water box in LAMMPS guide, H2O molecules look strange in dump file

LAMMPS LAMMPS Beginners
,
1

Dear all,

LAMMPS Model being used: 2024-03-02 18:23 LAMMPS-Win10-64bit-GUI-stable.exe

I am attempting to follow a water box guide from Polymer in water -
I have managed to run the file after altering the code slightly.
The code I changed was the following:

molecule h2omol H2O-SPCFw.mol
create_atoms 0 random 1050 87910 NULL mol &
    h2omol 454756 overlap 1.0 maxtry 50

Changed to:

molecule h2omol H2O-SPCFw.mol
create_atoms 0 random 1050 87910 NULL mol h2omol 454756 overlap 1.0 maxtry 50

When I try to view the dump file the hydrogen and oxygen molecules are not bonded in a proper manner. Below is a picture of what it looks like.
H20 dump

I will also provide all the code I am using to run the simulation here.

units real #masses are in grams per mole, distances in Ångstroms, time in femtoseconds, energies in Kcal/mole.
atom_style full #each atom is a dot with a mass and a charge that can be linked by bonds, angles, dihedrals and/or impropers
#bond_style, angle_style, and dihedral_style commands define the potentials for the bonds, angles, and dihedrals used in the simulation, here harmonic.
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
pair_style lj/cut/coul/long 12 #atoms interact through both a Lennard-Jones (LJ) potential and Coulombic interactions. The value of 12Å is the cutoff.
kspace_style pppm 1e-5 #defines the long-range solver for the (long) Coulombic interactions. The pppm style refers to particle-particle particle-mesh 
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 1.0 angle yes  #cancels the Lennard-Jones interactions between the closest atoms of the same molecule.

#create a 3D simulation box of dimensions 9x3x3nm^3, and make space for 9 atom types (2 for the water molecule, and 7 for the polymer molecule)
#7 bond types, 8 angle types, and 4 dihedral types.
region box block -45 45 -15 15 -15 15
create_box 9 box &
bond/types 7 &
angle/types 8 &
dihedral/types 4 &
extra/bond/per/atom 3 &
extra/angle/per/atom 6 &
extra/dihedral/per/atom 10 &
extra/special/per/atom 14


#file contains all the parameters (masses, interaction energies, bond equilibrium distances, etc). types 8 and 9 are for water and the atoms of types 1 to 7 are for the polymer molecule
include PARM.lammps

#create water molecules
molecule h2omol H2O-SPCFw.mol #molecule template
create_atoms 0 random 1050 87910 NULL mol h2omol 454756 overlap 1.0 maxtry 50 


group H2O type 8 9 #organize the atoms of types 8 and 9 of the water molecules in a group named H2O
minimize 1.0e-4 1.0e-6 100 1000 
reset_timestep 0    



fix mynpt all npt temp 300 300 100 iso 1 1 1000 

#print the atom positions in a .lammpstrj file every 1000 steps (i.e. 1 ps), print the temperature volume, 
#and density every 100 steps in 3 separate data files, and print the information in the terminal every 1000 steps

dump mydmp all atom 1000 dump.lammpstrj 
variable mytemp equal temp
variable myvol equal vol
fix myat1 all ave/time 10 10 100 v_mytemp file temperature.dat #print the temperature
fix myat2 all ave/time 10 10 100 v_myvol file volume.dat #print the volume
variable myoxy equal count(H2O)/3 

write_data H2O.data
2

How did you create this visualization?

3

Please see the forum guidelines post for details of how to correctly quote text files. Yours is partially misformatted.

4

Hello,

Isn’t that just a bad choice of VMD representation? What if you choose DynamicBonds or VDW ?
The fact that you had to change the SEED to make it work is a concern, I will take a look in the future.
Simon

5

I used CPK as the drawing method, VMD describes it as " a combination of both Bonds' and VDW’ in that it draws the atoms as spheres and the bonds as cylinders."
I may not understand completely but if I choose VDW it would just draw the atoms as spheres and leave out the bonds. Dynamic bonds guesses bonds based on distance.

6

It was created in VMD using CPK as the drawing method.

7

I will review the forum guidelines and edit my text files in the post.

8

The problem that you see is caused by the LAMMPS trajectory file having no information about which elements you have and no information about bonds.

For typical VMD use cases (i.e. processing trajectories created by NAMD using the CHARMM force field) you have a .psf file with the bond topology and element and type information and then a .dcd file with coordinate data.

The bond topology in LAMMPS is in the data file. Through the TopoTools plugin in VMD you can load data files and then use the VMD atom selection language to convert the numeric types that LAMMPS uses to convert them into element names. An automated method based on atom masses is part of TopoTools. You can then save the topology into a .psf file and load it together with the LAMMPS trajectory file. Mind you, when you write a trajectory with “wrapped” coordinates, you will still see those unexpected long bonds, since you have bonded atoms wrapped around periodic boundaries. The solution to this is to either write out unwrapped coordinates or also write out the image flags, or you can use the DynamicBonds representation in combination with VDW (scaled down) to visualize only the bonds within the box.

Another option would be to use OVITO for visualization. It is more in tune with what LAMMPS needs and provides than VMD.

9

Thank you for your response and time. I will look into all of this.

10

@Micah_Watson and future googlers… if you compile lammps with the code in PR 4245, then you have a write_psf command that will create a .psf to import into visualization software like vmd, chimerax (my personal favourite), and MolecularNodes for Blender (highest level of difficulty but best render quality). Please note you need to import topology information with read_psf first, otherwise you will get an error.