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
++++=============================+++++