Solvate a water box using lammps

LAMMPS LAMMPS General Discussion
1

Dear LAMMPS Users

I am trying to create a box of solvated water (NaCl) with a submerged nanostructure. I have created the water box and submerged the nanostructure, however I still can’t dissolve the water.

1.-My idea is to randomly remove water molecules according to the amount needed to solvate the volume of the box using the following syntax:

group water type 1 2
group tube type 3 4
delete_atoms random count 150 yes all water 482793 mol yes

However I get this error:

ERROR: Illegal delete_atoms command (…/delete_atoms.cpp:72)

2.-If I can remove the necessary water molecules, then I am going to replace the removed water molecules with sodium and chlorine ions using the following syntax:

variable counter equal count(water)
variable wtmol equal v_counter/3
variable wtion equal v_wtmol*0.15/55.5556 #Number of Na+ Ions
variable wtion1 equal round(v_wtion)

create_atoms 5 random ${wtion1} 12345 water overlap 2.0 maxtry 50

I also get problem in this syntax:

ERROR: Create_atoms region ID does not exist (…/create_atoms.cpp:104)

I need some ideas to achieve my goal please.

Sincerely

Eduardo Estevez

2

Both of these stem from the fact that you don’t know what a region is. Both commands require a region, but you give them a group, “water”.

3

This sounds very confusing. Water is the solvent so you cannot “solvate” it. From the context it seems that what you want to do is dissolve the NaCl in water.

Besides the syntax issues that were already pointed out, neither of your proposed strategies will work. The problem is that you have a dense system and you cannot easily find the empty spots where you can insert the ions after you delete the water molecules.

The way how this can be done is the following:

  • pick a set of water oxygen atoms and change their type and charge to that of the desired ions
  • now find all hydrogens attached to those former oxygens (now ions) and delete those and also the bonds and angles they are involved in.

This could be done with VMD scripting, for example. That is what I have been done in the past. There is a plugin (topotools) that can help with some of the steps. Most kinds of system manipulations that require changes to the topology are easier done with standalone tools.

However, this is also possible from within LAMMPS, but it is a bit complex. The following assumes that water oxygen atoms are type 1, water hydrogen atoms are type 2, sodium cations are type 3, chloride anions are type 4, and that all force field parameters and masses are properly set up. This will convert 20 randomly picked water molecules into 10 sodium and 10 chloride atoms.

group owat type 1
group hwat type 2
# change exclusions temporarily so that bonded atoms are in neighbor list
special_bonds lj/coul 1.0 1.0 1.0
# convert water oxygens to ions
set type 1 type/subset 3 10 6357834
set type 1 type/subset 4 10 6357834
# define groups for ions and set charges
group sod type 3
group chl type 4
set group sod charge 1.0
set group chl charge -1.0
# delete hydrogens (still) bound to ions and their bonds and angles
delete_atoms overlap 1.1 hwat sod compress no bond yes mol no
# this is needed to reset the internal neighbor list data to work around a current limitation in LAMMPS
run 0 post no
delete_atoms overlap 1.1 hwat chl compress no bond yes mol no
# restore exclusion settings to default and renumber atoms
special_bonds   lj/coul 0.0 0.0 0.0
reset_atom_ids

After this, there should be 10 water molecules changed to sodium, 10 water molecules to chloride and thus in total 40 fewer atoms (2x20 hydrogens), 40 fewer bonds (of the deleted hydrogens) and 20 fewer angles (involving the deleted hydrogens).

4

Thank you Dr Axel

I tried to use the code that you helped me, however I get an unexpected error:

ERROR: Invalid value in set command (…/set.cpp:134)
Last command: set type 1 type/subset 3 10 6357834

I attach the inpuscript. I understand that the “set” command converts oxygen to sodium ions, but LAMMPS indicates as an error:

set type 1 type/subset 3 10 6357834

Sincerely

Eduardo Estevez

#======================= Input script ==========================================

units real
boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic

pair_style lj/cut/coul/long 12

variable basename string CVFF

read_data ${basename}.lammpsdata

pair_modify tail no
kspace_style pppm 1e-6

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes one 11330 page 200000
timestep 1.0

group owat type 1
group hwat type 2

change exclusions temporarily so that bonded atoms are in neighbor list

special_bonds lj/coul 1.0 1.0 1.0

convert water oxygens to ions

set type 1 type/subset 3 10 6357834
set type 1 type/subset 4 10 6357834

define groups for ions and set charges

group sod type 3
group chl type 4
set group sod charge 1.0
set group chl charge -1.0

delete hydrogens (still) bound to ions and their bonds and angles

delete_atoms overlap 1.1 hwat sod compress no bond yes mol no

this is needed to reset the internal neighbor list data to work around a current limitation in LAMMPS

run 0 post no
delete_atoms overlap 1.1 hwat chl compress no bond yes mol no

restore exclusion settings to default and renumber atoms

special_bonds lj/coul 0.0 0.0 0.0
reset_atom_ids

write_data data.water_NaCl
run 0

5

are the atom types for the ions defined in the data file?
if not, you need to add extra/atom/types 2 to your read_data command.

6

Dr. Axel, if I actually put 4 types of atoms in the input data and put the extra/atom/types 2 command, but I get the error:

ERROR: Expected floating point parameter instead of ‘lj/cut/coul/long’ in input script or data file (…/pair_lj_cut_coul_long.cpp:628)

I enclose the following link with the inputscript and the data, since I do not have the permissions to attach directly to the forum.

https://drive.google.com/drive/folders/1S9U7w3IGtHyfIcE9Fu1M8AV1wN35MRCh?usp=sharing

Sincerely

Eduardo Estevez

7

Please think before making changes and carefully study the documentation. You should do only one of those changes, not both: that will create new problems.

Please do your own debugging! I don’t have the time to fix everybody’s inputs. You should only need to ask for help with your inputs when the behavior disagrees with the documentation. In those cases, however, you should also not “dump” your entire input deck on people, but rather try to construct a very minimal input that only produces the error and point out where it conflicts with the documentation.

Normally LAMMPS will print the offending line, so you can check it.
If this is while reading the data file, then this may be due to incorrectly formatting the Pair Coeffs section (obviously).

8

Thanks for your patience Dr. Axel, I had a little syntax error in the potential pairs (I shouldn’t put the name of the potential style in each pair of values when it is a single type).

Finally it worked correctly, the option I chose was to leave only the hydrogen and oxygen information of the water in the input data, and then inside the inputscript place the command:

read_data file.data extra/atom/types 2

## add the value of the masses of the sodium and chlorine ions

mass 3 22.98977 # SOD
mass 4 35.45 #CL

The only part I couldn’t get was resetting the number of atoms, I don’t know if it’s the version of my LAMMPS (Mar 3, 2020), I got the following error:

ERROR: Unknown command: reset_atom_ids (…/input.cpp:232)
Last command: reset_atom_ids

Sincerely

Eduardo Estevez

9

Yes, you should use a more recent version of LAMMPS.