Rerun using GCMC trajectory

LAMMPS LAMMPS General Discussion
,
1

Hello,

I’m attempting to calculate the potential energy between adsorbed water molecules and a zeolite framework obtained from GCMC simulation (I used GOMC to generate the coordinates, not LAMMPS) using the rerun functionality in LAMMPS so I can omit the framework-framework and water-water interactions. The ultimate goal is to calculate the heat of adsorption at a given temperature and pressure.

I’ve gotten the calculation up and running and added in a dynamic variable to keep track of the water oxygen atoms. I also included the thermo variable “atoms” in my thermo_style for good measure. To my dismay both of these return a constant number of atoms that reflects the atom count in the first frame of the trajectory. But the potential energy does seem to be different for each frame. Is the calculation carrying on correctly but the atom count just isn’t updating?

To initialize the system I generated a data file from the first frame of the trajectory. Additionally, I’m aware that due to the long range electrostatic interactions I need to run additional separate calculations to remove the water-water and zeolite-zeolite contributions.

Below is my input and the corresponding log file showing the constant atom count. I’m using version LAMMPS (27 Oct 2021).

Input
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
units real
atom_style full
#bond_style harmonic
#angle_style harmonic
pair_style lj/cut/coul/long 12.0 12.0
pair_modify tail yes
kspace_style ewald 0.000001
boundary p p p

read_data bea_a-silica-tip5pe-water-25C-7240Pa.data
mass 1 15.999 # O TIP5P-E Water
mass 2 1.0078 # H TIP5P-E Water
mass 3 1e-30 # M TIP5P-E Water
mass 47 1.0 # Placeholders
mass 8 15.999 # O Zeo
mass 913 1.0 # Placeholders
mass 14 28.09 # Si Zeo
#pair_coeff 1 1 0.16 3.12 # O TIP5P-E Water
pair_coeff 1 1 0.0 0.0 # O-O TIP5P-E Water
pair_coeff 1 2 0.0 0.0 # O-H TIP5P-E Water
pair_coeff 1 3 0.0 0.0 # O-M TIP5P-E Water
pair_coeff 2 2 0.0 0.0 # H TIP5P-E Water
pair_coeff 2 3 0.0 0.0 # H-M TIP5P-E Water
pair_coeff 3 3 0.0 0.0 # M TIP5P-E Water
pair_coeff 1 8 0.1283745 3.0655 # O Water - O Zeo
pair_coeff 2 8 0.0 0.0 # H Water - O Zeo
pair_coeff 3 8 0.0 0.0 # M Water - O Zeo
pair_coeff 1 14 0.1011929 3.045 # O Water - Si Zeo
pair_coeff 2 14 0.0 0.0 # H Water - Si Zeo
pair_coeff 3 14 0.0 0.0 # M Water - Si Zeo
pair_coeff 4 4 0.0 0.0 # Necessary placeholders
pair_coeff 5 5 0.0 0.0 # Necessary placeholders
pair_coeff 6 6 0.0 0.0 # Necessary placeholders
pair_coeff 7 7 0.0 0.0 # Necessary placeholders
pair_coeff 8 8 0.0 0.0 # O Zeo
pair_coeff 8 14 0.0 0.0 # O Zeo - Si Zeo
pair_coeff 9 9 0.0 0.0 # Necessary placeholders
pair_coeff 10 10 0.0 0.0 # Necessary placeholders
pair_coeff 11 11 0.0 0.0 # Necessary placeholders
pair_coeff 12 12 0.0 0.0 # Necessary placeholders
pair_coeff 13 13 0.0 0.0 # Necessary placeholders
pair_coeff 14 14 0.0 0.0 # Si Zeo

group framework type 8 14
group water type 1 2 3
neigh_modify every 1 delay 0 check yes page 10000 one 1000
thermo 1
set type 1 charge 0.0 # O TIP5P-E Water
set type 2 charge 0.241 # H TIP5P-E Water
set type 3 charge -0.241 # M TIP5P-E Water
set type 47 charge 0.0 # Placeholders
set type 8 charge -0.7065 # O Zeo
set type 913 charge 0.0 # Placeholders
set type 14 charge 1.413 # Si Zeo
variable water_O atom “type==1”
group water_O dynamic all every 1 var water_O
variable nOw equal count(water_O)
thermo_style custom step pe atoms v_nOw
variable PE_mol equal pe
fix avg2 all ave/time 1 1 1 v_PE_mol file thermo-props-rerun-all-charges.dat
rerun bea_a-silica-tip5pe-water-25C-7240Pa.lammpstrj first 0 dump x y z box yes

Log
#~~~~~~~~~~~~~~~~~~~~~~~~~~
LAMMPS (27 Oct 2021)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
Reading data file …
orthogonal box = (0.0000000 0.0000000 0.0000000) to (37.983000 37.983000 52.812000)
1 by 1 by 1 MPI processor grid
reading atoms …
3491 atoms
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.020 seconds
read_data CPU = 0.033 seconds
3456 atoms in group framework
35 atoms in group water
Setting atom values …
7 settings made for charge
Setting atom values …
14 settings made for charge
Setting atom values …
14 settings made for charge
Setting atom values …
0 settings made for charge
Setting atom values …
2304 settings made for charge
Setting atom values …
0 settings made for charge
Setting atom values …
1152 settings made for charge
dynamic group water_O defined
Ewald initialization …
using 12-bit tables for long-range coulomb (src/kspace.cpp:340)
G vector (1/distance) = 0.28520638
estimated absolute RMS force accuracy = 0.00040031343
estimated relative force accuracy = 1.205532e-06
KSpace vectors: actual max1d max3d = 5313 17 21437
kxmax kymax kzmax = 12 12 17
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 1000, page size: 10000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 6 6 8
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 39.32 | 39.32 | 39.32 Mbytes
Step PotEng Atoms v_nOw
0 -522630.78 3491 7
1 -522630.33 3491 7
2 -522633.58 3491 7
3 -522631.16 3491 7
4 -522629.88 3491 7
5 -522640.83 3491 7
6 -522631.59 3491 7
7 -522633.07 3491 7
8 -522633.07 3491 7
9 -522627.94 3491 7
10 -522627.84 3491 7
11 -522627.84 3491 7
12 -522634.65 3491 7
13 -522630.27 3491 7
14 -522626.4 3491 7
15 -522632.24 3491 7
16 -522632.24 3491 7
17 -522632.24 3491 7
18 -522632.76 3491 7
19 -522632.69 3491 7
20 -522632.32 3491 7
Loop time of 4.34358 on 1 procs for 21 steps with 3491 atoms

Total wall time: 0:00:04

Any help would be appreciated!

2

I don’t see how this can work. If you add or remove molecules, you also need to add or remove the corresponding bonds and angles. However, there is no information about those stored in trajectory files, so this information is not available and thus you cannot recover the exact system that was present during the simulation.

The behavior of how to deal with added or removed atoms depends on the settings for dump file reader. These can be changed by using the “dump” keyword to “rerun” and you can see the available options in the documentation of the “read_dump” command. But as mentioned above, you cannot restore bond/angle information, since that is not present in the dump file. So even if you change the settings to create added atoms, you will not be able to post-process the trajectory correctly.

3

Thanks for the quick response!

I should have specified, I’m using the rigid TIP5P-E water model and a rigid framework model so I’m not really concerned about the bond/angle information. I was doing this with the hope that I could simply calculate the potential energy between water and zeolite atoms based on their non-bonded interactions.

Blockquote
The behavior of how to deal with added or removed atoms depends on the settings for dump file reader. These can be changed by using the “dump” keyword to “rerun” and you can see the available options in the documentation of the “read_dump” command. But as mentioned above, you cannot restore bond/angle information, since that is not present in the dump file. So even if you change the settings to create added atoms, you will not be able to post-process the trajectory correctly.
Blockquote

I did take note of these features on the “rerun” and “read_dump” pages. If I’m understanding your response correctly, these options should be able to handle an atom count that changes each frame? Or do they only take affect on the first frame? Sorry if I’m being a bit naive.

Thanks!

Edit: I added the trim, replace, and add options to my rerun command and it appears to be working correctly now.

4

Did you dump and restore charges? I didn’t see them in your quoted input. If not, you may have LAMMPS delete and create atoms, but the Coulomb interactions would not be the same. It is probably a good idea to check if the rerun can reproduce the potential energy values reported by the original GCMC run. Your input does have “set” commands to set charges, but those will be applied immediately and only to the atoms that are present at the time.

5

I’ll check the potential energies produced by both methods, that’s a good point.

In terms of the “dump and restore charges” comment, how exactly would that work? Would I need to have the charges already be present within the trajectory I’m passing through “rerun” or is there a different approach I can implement with my current input?

6

Yes. Properties that can be read back are listed in the read_dump documentation.

I cannot think of any.

Since in your case the charge is tied to the atom type, you could post-process your trajectory file to insert a column with the charges.

7

I’m still working on the potential energy validation between the two methods but for the sake of cataloguing the post-processing for others who may find it useful, I used the following awk and sed commands and so far everything is behaving as expected.

  • This assumes you only have x y z id and type as columns in your trajectory and want to add charges. Obviously you would need to adjust accordingly if you have additional quantities.

  • My overall approach here involved adding an extra column for the charges to the end of each relevant line “tmp” and using sed to edit in the real charge values.

  • I did also include two additional columns “tmp tmp” on each line so those lines would be different from everything else in the file, this way I could isolate them using awk conditions.

  • I then rearranged the lines containing atom coordinates and charges so they would be in the order of “x y z q id type” and then threw out the redundant lines that resulted from the previous awk commands.

  • Finally I changed the headers above the coordinates to read “ITEM: ATOMS x y z q id type”

awk ‘{print} {if (NF==5) {print $0, “tmp tmp tmp”}}’ file1 > file2

sed -i ‘’ “s/14 tmp/14 1.413/g;s/8 tmp/8 -0.7065/g;s/1 tmp/1 0.0/g;s/2 tmp/2 0.241/g;s/3 tmp/3 -0.241/g” file

awk ‘{print} {if (NF==8) {print $1, $2, $3, $6, $4, $5}}’ file1 | awk ‘{if (NF!=5) {print}}’ | awk ‘{if (NF!=8) {print}}’ > file2

sed -i ‘’ “s/ITEM: ATOMS x y z id type/ITEM: ATOMS x y z q id type/” file

(Here 1, 2, 3, 8, and 14 are my atom types and the numbers that look like charges in the first sed command are the charges)

I’m sure there’s a more elegant way to accomplish this. Regardless, there you go.