Adding hydroxide ion (OH-) and ammonium ion (NH4+) to the simulation box

LAMMPS LAMMPS General Discussion
1

Hi all,

Using VMD I filled between two graphene sheets with a water box, but VMD does not have NH4+ and OH- ions to ionize my water box. I want to use the topo tools plugin of vmd to produce a LAMMPS data file for this system.
I think we can include Param and Top files for NH4+ to the charmm toppar but it is not a straightforward task. I used OpenBabel to generate the PDB file of the desired ions and then used packmol to make the water box with the ion inside that. I can open that pdb file in vmd and everything looks fine but the pdb file format doesn’t look fine. Part of the pdb which is the end of H2O atoms and the beginning of OH is shown below.

ATOM 53994 O HOH B4999 2.995 23.380 -17.954 0.00 0.00
ATOM 53995 H HOH B5000 12.807 -23.484 64.827 0.00 0.00
ATOM 53996 H HOH B5000 11.194 -23.428 64.585 0.00 0.00
ATOM 53997 O HOH B5000 11.982 -22.932 64.948 0.00 0.00
ATOM 53998 H HOH B5001 -25.043 24.943 55.807 0.00 0.00
ATOM 53999 H HOH B5001 -25.452 26.018 56.966 0.00 0.00
ATOM 54000 O HOH B5001 -25.790 25.485 56.192 0.00 0.00

ATOM 54001 O HOH C 1 16.958 12.226 47.560 1.00 0.00 O
ATOM 54002 H HOH C9999 16.028 12.102 47.329 1.00 0.00 H
ATOM 54003 O HOH C 1 -23.213 -1.937 48.734 1.00 0.00 O
ATOM 54004 H HOH C9999 -22.898 -2.053 47.828 1.00 0.00 H
ATOM 54005 O HOH C 1 7.618 -18.981 -36.418 1.00 0.00 O
ATOM 54006 H HOH C9999 7.110 -18.458 -35.783 1.00 0.00 H

As it is shown the residue type for both is similar, HOH, and the residue sequence number for OH- are either 1 or 9999. My question is that is there any way to fix this problem like manually changing the residue type and sequence in pdb or the LAMMPS data file? Or no I should include the topology of the NH4+ and OH- negative to the charmm rtf file and regenerate the psf, and pdb file using that? I will really appreciate it if someone can help me with this.

I used Topo Tools to produce the LAMMPS data file but because my pdb has these issues it also has a similar problem like my pdb

Thank you
Sadra

2

Before making any suggestions one question:
were you able to get a properly working simulation with pure water?

3

I did not try with the pure water. You are right I should try it first, but I was too much obsessed with the ions that I forgot about that. I will try water alone without graphene and ions.

4

Then you should get that working first. It is pointless to discuss how to get ions into your system if you have not successfully done the simpler case first.

5

You are right, I will try it
thank you

6

Dear Dr. Kohlmeyer,

I considered just the water box and the ions.

1- I revised the PDB files for each component

2- Then I used Packmole to generate my initial setup

3- used topotools to generate lammps input file

4: set the LJ parameters as below:

I obtained the LJ parameters for OH- from toppar_water_ions.mod.str and for ammonium ion I consider the parameter of methylammonium (KK 03/10/92)) par_all22_prot.inp based on other people experience in here, here, here and some other discussion in CHARMM Forum.

I- Could you please let me know if my approach is correct?

II- I ran the script for a short time and a test system (500 water molecules and 10 molecules of each one of ions). The simulation was completed without any error but two of the water molecules are breaking apart. I am not sure if they are still in place and just the bonds are extended due to the PBC? I will appreciate any comment about that.

++++==============================++++
units real
dimension 3
#newton off
#processors * * *
boundary p p f
atom_style full

##Set Potential Parameters_______

pair_style lj/cut/coul/long 10
bond_style harmonic
angle_style harmonic
kspace_style pppm 1.0e-6
kspace_modify slab 3.0

read_data lammps-ion-wat-500.dat
include param.lammps

group wat type 2 5
group hyox type 3 6
group amon type 1 4

neighbor 3.0 bin
#neigh_modify delay 0 every 1 check yes
neigh_modify every 1 delay 5 check yes

thermo_style multi
thermo 1000
minimize 1.0e-4 1.0e-6 1000 1000
write_data data-min-wat-ions.lammps
dump mydmp all atom 1000 dump.lammpstrj
dump dcd_1 all dcd 1000 wat-ion.dcd

###run
run_style verlet
timestep 0.01 ##0.001
velocity all create 300.0 4928454 rot yes dist gaussian

fix 0 all nvt temp 300.0 300.0 100
#fix 1 wat shake 1.0e-4 20 10 b 2 a 2
#fix 2 hyox shake 1.0e-4 20 10 b 3
#fix 3 amon shake 1.0e-4 20 10 b 1 a 1 ## cannot do it in lammps^
fix 4 all shake 1.0e-4 20 10 b 2 3 a 1 2

++++=============================+++++

7

Your approach is correct, if it produces a correct system topology, geometry and correct simulation results. It is impossible to determine this from just your description or looking at some scripts. I would have to redo all your work and verify each step individually. This would be something that your should be doing in coordination with your adviser/supervisor/tutor. I don’t have the time to do this.

Did you do a pure water simulation? Did this have the same “symptoms”? With bond style harmonic it is “impossible” for bonds to break (on their own). Are you certain that what you are observing is not an artifact of your visualization? Please note that dump files do not contain topology information, so a visualization may not show the “real” bonds but only apply an empirical, distance based heuristic.
Overall, verifying and validating simulation results (on small test systems), is standard procedure when doing MD simulations and something that you should be doing on your own. As mentioned before, if you need help with this you need to ask your adviser/supervisor/tutor/colleagues for help.

8

Thank you for your reply
Yes I tested with water. I am not getting any error but some of the water molecules are breaking apart whihc as you said It could be due to visulaization tools.
Thank you for your help

9

And also keep in mind that you are using a rigid water model and keep it constrained with fix shake (if used correctly, that is), so the bonds should not even change in length. Since you have periodic boundary conditions, you may have a situation where the two atoms of a bond are on opposite sides of the simulation box, since you output “wrapped” coordinates.

1 Like