Respected All,
I have written a script in Lammps, for GCMC. My system is adsorption of gas on Zeolite (taken from literature). After running the simulation, I am getting an error “Lost atoms”.
Here I am pasting the molecule input, LAMMPS input, data files, and the Output for your reference. It would be helpful to me if someone could suggest a way to clear the issue.
p.s: In the data file, I have included the coordinates of 6912 atoms that belong to the framework, while the remaining 4 atoms belong to the molecule.
LAMMPS version: September 2015 (colvars)
Thank you,
Yours sincerely,
Vardhan
Input:
Chabazite
units real
dimension 3
atom_style full
bond_style none
angle_style none
boundary p p p
pair_style lj/cut/coul/long 14.0 14.0 ## cut-off distance
pair_modify mix arithmetic shift no tail yes
kspace_style ewald 1.0e-6
read_data zeolite.data ## zeolite
region chabazite prism 0.0000 54.1169 0.0000 46.8666 0.0000 58.9931 -27.0584 0.0000 0.0000 side in units box
molecule h2smol h2s.txt
pair_coeff 1 1 0.0441 2.310 14.0 ## silicon
pair_coeff 2 2 0.1057 3.304 14.0 ## oxygen
pair_coeff 1 1 0.2424 3.600 14.0 ## S_h2s
pair_coeff 2 2 0.0994 2.500 14.0 ## H_h2s
pair_coeff 3 3 0.0000 0.000 14.0 ## S_com
mass 3 32.0650 ## S_h2s
mass 4 1.00794 ## H1_h2s
mass 5 1.e-20 ## S_com
group chabazite type 1 2
group h2s type 3 4 5
fix freeze chabazite setforce 0.0 0.0 0.0
velocity chabazite set 0.0 0.0 0.0 sum no
neighbor 2.0 bin
neigh_modify every 1 delay 10 check yes
velocity h2s create 298.0 9821144
timestep 1.0
fix h2s_rigidnvt h2s rigid/nvt/small molecule temp 298.0 298.0 1 mol h2smol
fix_modify h2s_rigidnvt dynamic/dof yes
fix h2s_gcmc h2s gcmc 1 100 100 0 9821144 298.0 1.0 20.0 mol h2smol pressure 0.986923 fugacity_coeff 0.992430 full_energy
variable sulfur atom “type==3”
variable hydrogen atom “type==4”
variable offsite atom “type==5”
group sulfur dynamic h2s var sulfur
group hydrogen dynamic h2s var hydrogen
group offsite dynamic h2s var offsite
variable nS equal count(sulfur)
variable nH equal count(hydrogen)
variable nX equal count(offsite)
variable tacc equal f_h2s_gcmc[2]/(f_h2s_gcmc[1]+0.1)
variable iacc equal f_h2s_gcmc[4]/(f_h2s_gcmc[3]+0.1)
variable dacc equal f_h2s_gcmc[6]/(f_h2s_gcmc[5]+0.1)
variable racc equal f_h2s_gcmc[8]/(f_h2s_gcmc[7]+0.1)
compute_modify thermo_temp dynamic/dof yes
thermo_modify lost ignore
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc v_nS v_nH v_nX
thermo 1
dump h2s_traj all xyz 2 dump.xyz
dump h2s_traj all image 2 dump.jpeg type type
dump_modify h2s_traj flush yes
run 20
molecule:
h2s molecule
4 atoms
3 bonds
3 angles
Coords
1 0.0000 0.0000 0.0000 ## S_h2s
2 0.9640 -0.9309 0.0000 ## H1_h2s
3 -0.9640 -0.9309 0.0000 ## H2_h2s
4 0.0000 -0.3000 0.0000 ## S_com
Types
1 1 ## S_h2s
2 2 ## H1_h2s
3 2 ## H2_h2s
4 3 ## S_com
Charges
1 0.00 ## S_h2s
2 0.21 ## H1_h2s
3 0.21 ## H2_h2s
4 -0.42 ## S_com
Bonds
1 1 1 2
2 1 1 3
3 2 1 4
Angles
1 1 2 1 3
2 2 2 1 4
3 2 4 1 3
Data file:
Chabazite
6916 atoms ## number of atoms
3 bonds
3 angles
0 dihedrals
0 impropers
5 atom types ## atom types, Si & O
2 bond types
2 angle types
0 dihedral types
0 improper types
0.0000 54.1169 xlo xhi
0.0000 46.8666 ylo yhi
0.0000 58.9931 zlo zhi
-27.0584 0.0000 0.0000 xy xz yz ## tilt factors for triclinic
Masses
1 28.0855 ## silicon
2 15.9990 ## oxygen
3 32.0650 ## S_h2s
4 1.00794 ## H1_h2s
5 1.e-20 ## S_com
Atoms
1 1 1 1.5 17.2009 42.9612 55.6093
2 1 1 1.5 21.8405 45.6394 55.6093
3 1 1 1.5 21.8401 40.2823 55.6093
.
.
.
.
.6911 1 2 -0.75 4.3759 0.2318 2.4581
6912 1 2 -0.75 -4.3338 10.3132 1.9160
6913 2 3 0.00 0.1000 0.1000 0.0000
6914 2 4 0.21 1.4400 0.1000 0.0000
6915 2 5 -0.42 0.3084 0.3158 0.0000
6916 2 4 0.21 0.0532 1.4392 0.0000
Pair Coeffs
1 0.0441 2.310 ## silicon
2 0.1057 3.304 ## oxygen
3 0.2424 3.600 ## S_h2s
4 0.0994 2.500 ## H_h2s
5 0.0000 0.000 ## S_com
Bonds
1 1 6913 6914
2 1 6913 6916
3 2 6913 6915
Angles
1 1 6914 6913 6916
2 2 6914 6913 6915
3 2 6915 6913 6916
Output
Ewald initialization …
G vector (1/distance) = 0.245045
estimated absolute RMS force accuracy = 0.000374866
estimated relative force accuracy = 1.1289e-06
KSpace vectors: actual max1d max3d = 6652 20 34460
kxmax kymax kzmax = 15 20 16
0 atoms in group FixGCMC:gcmc_exclusion_group:h2s_gcmc
0 atoms in group FixGCMC:rotation_gas_atoms:h2s_gcmc
WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (…/neighbor.cpp:441)
Neighbor list info …
1 neighbor list requests
update every 1 steps, delay 10 steps, check yes
master list distance cutoff = 16
ghost atom cutoff = 16
Memory usage per processor = 18.508 Mbytes
Step Temp Press PotEng KinEng Density Atoms iacc dacc tacc racc nS nH nX
0 0.015904701 3747097.1 906768.17 0.32773768 1.5367248 6916 0 0 0 0 1 2 1
ERROR: Lost atoms: original 7128 current 6912 (…/thermo.cpp:395)
THE END