[lammps-users] GCMC of CO2 in zeolite NaY

Dear lammps users,
I am using Lammps March 3, 2020 version. I am running the GCMC simulation of CO2 in zeolite NaY. But after running, I am getting the following error.
ERROR on proc 0: Shake determinant = 0.0 (src/RIGID/fix_shake.cpp:2315)
I input file as follows. I also attached the data file and CO2_zeo.txt file. Kindly help me to rectify this problem. Thank you well in advance.

#LAMMPS Input file for GCMC CO2 adsorbed on zeolite NaY

zeo 2.43

#variable available on command line

variable mu index -8.1
variable disp index 0.1
variable temp index 300.0
variable lbox index 24.8536
variable spacing index 1.0

Intialization

units real
dimension 3
boundary p p p
atom_style full
atom_modify map array sort 500 2.0

Atom Definition

#read_restart restart.500000.equil

read_data data.zeoco2
replicate 1 1 1

box, start molecules on simple cubic lattice

lattice sc {spacing} region box block 0 {lbox} 0 {lbox} 0 {lbox} units box
#create_box 2 box &

bond/types 1 &

angle/types 1 &

extra/bond/per/atom 2 &

extra/angle/per/atom 1 &

extra/special/per/atom 2

molecule co2mol CO2_zeo.txt
#create_atoms 0 box mol co2mol 464563 units box

bond_style harmonic
bond_coeff 1 300.290000000000 1.600000000000000 #Si-O
bond_coeff 2 222.200000000000 1.70000000000000 #Al-O
bond_coeff 3 30.0200000000000 3.12000000000000 #Si-O-Si
bond_coeff 4 30.0200000000000 3.18000000000000 #Al-O-Si
bond_coeff 5 1000.0000000000000 1.16100000000000 #C-O dummy K

angle_style harmonic
angle_coeff 1 7.01000000000000 149.500000000000 #Si-O-Si
angle_coeff 2 5.00000000000000 149.500000000000 #Al-O-Si
angle_coeff 3 75.0700000000000 109.500000000000 #O-Si-O
angle_coeff 4 65.0500000000000 109.500000000000 #O-Al-O
angle_coeff 5 100.0000000000000 109.490000000000 # O-C-O dummy K

#pair_style lj/long/coul/long long long 12.0
pair_style lj/cut/coul/long 12.0
kspace_style pppm 1.0e-4
#kspace_style ewald/disp 1.0e-4
pair_coeff 1 1 0.6004 3.9200
pair_coeff 1 2 0.6250 3.9200
pair_coeff 1 3 0.3032 3.5350
pair_coeff 1 4 0.3032 3.5350
pair_coeff 1 5 0.3101 3.1750
pair_coeff 1 6 0.1860 3.3525
pair_coeff 1 7 0.3142 3.4920
pair_coeff 2 2 0.6507 3.9200
pair_coeff 2 3 0.3156 3.5350
pair_coeff 2 4 0.3156 3.5350
pair_coeff 2 5 0.3229 3.1750
pair_coeff 2 6 0.1936 3.3525
pair_coeff 2 7 0.3271 3.4920
pair_coeff 3 3 0.1531 3.1500
pair_coeff 3 4 0.1531 3.1500
pair_coeff 3 5 0.1566 2.7900
pair_coeff 3 6 0.0939 2.9675
pair_coeff 3 7 0.1586 3.1070
pair_coeff 4 4 0.1531 3.1500
pair_coeff 4 5 0.1566 2.7900
pair_coeff 4 6 0.0939 2.9675
pair_coeff 4 7 0.1586 3.1070
pair_coeff 5 5 0.1602 2.4300
pair_coeff 5 6 0.0961 2.6075
pair_coeff 5 7 0.1623 2.7470
pair_coeff 6 6 0.0576 2.7850
pair_coeff 6 7 0.0973 2.9245
pair_coeff 7 7 0.1644 3.0640

group zeolite type 1:5
group co2 type 6:7
group freeze id 577:632

mass of atoms of co2

mass 6 12.0111
mass 7 15.9994

MD Settings

#velocity uranyl create 300.0 4928459 rot yes dist gaussian
#velocity nitrate create 300.0 4928458 rot yes dist gaussian
#velocity water create 300.0 4928457 rot yes dist gaussian

#velocity guest scale 200.0
velocity freeze set 0.0 0.0 0.0 units box
velocity co2 create ${temp} 54654

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.01

minimize 0.0 0.0 100 100 #1.0e-6 1.0e-6 1000 1000
reset_timestep 0
run 0

Rigid constraints with thermostst

fix mynvt all nvt temp {temp} {temp} 100
fix cshake co2 shake 0.0001 10 0 b 5 a 5 mol co2mol
#fix 3 water shake 0.0001 10 0 b 7 a 7
fix 4 freeze setforce 0.0 0.0 0.0
#unfix 3

important to make tempreture dofs dynamic

compute_modify thermo_temp dynamic/dof yes
compute_modify mynvt_temp dynamic/dof yes

run 1000
reset_timestep 0

gcmc

#unfix cshake
variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
fix mygcmc co2 gcmc 1 100 0 0 99 {temp} {mu} {disp} mol co2mol maxangle 180 full_energy & tfac_insert {tfac} group co2 shake cshake
#maxangle 180 full_energy
#fix cshake

atom counts

variable carbon atom “type==6”
variable oxygen atom “type==7”
group carbon dynamic all var carbon
group oxygen dynamic all var oxygen
variable nC equal count(carbon)
variable nO equal count(oxygen)

Output

variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)

thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc v_nC v_nO

thermo 1000

run 20000

With regards

Shamimul Ahsan

CO2_zeo.txt (591 Bytes)

data.zeoco2 (181 KB)

in.zeoco2_gcmc (5.23 KB)

You cannot use fix shake on linear molecules. please see: https://docs.lammps.org/fix_shake.html#restrictions

Hi,
I have used rigid instead of shake. It is now running but the number of molecules is not fluctuating since in GCMC simulation it is supposed to fluctuate.
My input file as follows. Kindly help me. But in system with out zeolite , number of atoms are fluctuating. Thanks well in advance. I have attached my outputfile also.

#LAMMPS Input file for GCMC CO2 adsorbed on zeolite NaY

zeo 2.43

#variable available on command line

variable mu index -10.0
variable disp index 0.5
variable temp index 300.0
variable lbox index 24.8536
variable spacing index 1.0

Intialization

units real
dimension 3
boundary p p p
atom_style full
#atom_modify map array sort 500 2.0

Atom Definition

#read_restart restart.500000.equil

read_data data.zeo_co2
#replicate 1 1 1

box, start molecules on simple cubic lattice

#lattice sc {spacing} #region box block 0 {lbox} 0 {lbox} 0 {lbox} units box
#create_box 2 box &

bond/types 1 &

angle/types 1 &

extra/bond/per/atom 2 &

extra/angle/per/atom 1 &

extra/special/per/atom 2

molecule co2mol CO2.txt
#create_atoms 0 box mol co2mol 464563 units box

bond_style harmonic
bond_coeff 1 300.290000000000 1.600000000000000 #Si-O
bond_coeff 2 222.200000000000 1.70000000000000 #Al-O
bond_coeff 3 30.0200000000000 3.12000000000000 #Si-O-Si
bond_coeff 4 30.0200000000000 3.18000000000000 #Al-O-Si
bond_coeff 5 1000.0000000000000 1.16100000000000 #C-O dummy K

angle_style harmonic
angle_coeff 1 7.01000000000000 149.500000000000 #Si-O-Si
angle_coeff 2 5.00000000000000 149.500000000000 #Al-O-Si
angle_coeff 3 75.0700000000000 109.500000000000 #O-Si-O
angle_coeff 4 65.0500000000000 109.500000000000 #O-Al-O
angle_coeff 5 200.0000000000000 180.00000000000 # O-C-O dummy K

#pair_style lj/long/coul/long long long 12.0
pair_style lj/cut/coul/long 12.0
kspace_style pppm 1.0e-4
#kspace_style ewald/disp 1.0e-4
pair_coeff 1 1 0.6004 3.9200
pair_coeff 1 2 0.6250 3.9200
pair_coeff 1 3 0.3032 3.5350
pair_coeff 1 4 0.3032 3.5350
pair_coeff 1 5 0.3101 3.1750
pair_coeff 1 6 0.1860 3.3525
pair_coeff 1 7 0.3142 3.4920
pair_coeff 2 2 0.6507 3.9200
pair_coeff 2 3 0.3156 3.5350
pair_coeff 2 4 0.3156 3.5350
pair_coeff 2 5 0.3229 3.1750
pair_coeff 2 6 0.1936 3.3525
pair_coeff 2 7 0.3271 3.4920
pair_coeff 3 3 0.1531 3.1500
pair_coeff 3 4 0.1531 3.1500
pair_coeff 3 5 0.1566 2.7900
pair_coeff 3 6 0.0939 2.9675
pair_coeff 3 7 0.1586 3.1070
pair_coeff 4 4 0.1531 3.1500
pair_coeff 4 5 0.1566 2.7900
pair_coeff 4 6 0.0939 2.9675
pair_coeff 4 7 0.1586 3.1070
pair_coeff 5 5 0.1602 2.4300
pair_coeff 5 6 0.0961 2.6075
pair_coeff 5 7 0.1623 2.7470
pair_coeff 6 6 0.0576 2.7850
pair_coeff 6 7 0.0973 2.9245
pair_coeff 7 7 0.1644 3.0640

group zeolite type 1:5
group co2 type 6:7
group freeze id 577:632

mass of atoms of co2

#mass 6 12.0111
#mass 7 15.9994

MD Settings

#velocity uranyl create 300.0 4928459 rot yes dist gaussian
#velocity nitrate create 300.0 4928458 rot yes dist gaussian
#velocity water create 300.0 4928457 rot yes dist gaussian

#velocity guest scale 200.0
velocity freeze set 0.0 0.0 0.0 units box
velocity co2 create ${temp} 54654

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
timestep 1.0

minimize 0.0 0.0 100 100 #1.0e-6 1.0e-6 1000 1000
reset_timestep 0
#run 0

Rigid constraints with thermostst

group notco2 subtract all co2
fix mynvt notco2 nvt temp {temp} {temp} 100
fix crigid co2 rigid/nvt/small molecule temp {temp} {temp} 100 mol co2mol

fix_modify mynvt dynamic/dof yes
fix_modify crigid dynamic/dof yes

#variable 10 co2mol CO2.txt
#fix 10 co2 rigid custom v_10
#fix cshake co2 shake 0.0001 10 0 b 5 a 5 mol co2mol
#fix 3 water shake 0.0001 10 0 b 7 a 7
fix 4 freeze setforce 0.0 0.0 0.0
#unfix 3

dynamically update fix rigid tempreture ndof

#fix_modify crigid dynamic/dof yes

important to make tempreture dofs dynamic

#compute_modify thermo_temp dynamic/dof yes
#compute_modify mynvt_temp dynamic/dof yes

#run 100
#reset_timestep 0

gcmc

#unfix cshake
variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
fix mygcmc co2 gcmc 100 100 100 0 99 {temp} {mu} {disp} mol co2mol & tfac_insert {tfac} group co2 #rigid crigid
#maxangle 180 full_energy
#fix cshake

atom counts

variable carbon atom “type==6”
variable oxygen atom “type==7”
group carbon dynamic all var carbon
group oxygen dynamic all var oxygen
variable nC equal count(carbon)
variable nO equal count(oxygen)

Output

variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)

compute_modify thermo_temp dynamic/dof yes

thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc v_nC v_nO

thermo 100

run 2000

log.lammps (15.2 KB)

in your fix gcmc command, the “rigid” option seems to be commented out. if you don’t have fix gcmc tell fix rigid that you are adding/removing molecules, there will be no update since fix rigid maintains the list of the rigid bodies.

Hi,
If I make rigid option active in fix gcmc command then i am getting the following error.
ERROR: Cannot use fix gcmc rigid with MC moves (src/MC/fix_gcmc.cpp:189).

With Regards
Shamimul Ahsan

By commenting the option out, you don’t fix anything. Fix gcmc will still no be able to do MC displacement moves of the rigid objects and will no be able to update fix rigid about inserted or deleted molecules.

The code is telling you about a real limitation.