Hello Developers/Users,
I am trying to simulate water desalination from C1N1 layer using lammps. The system comprises of water , water+ions seperated by C1N1 layer and the graphene layers on the end of water layer to apply pressure. For C1N1 and graphene i am trying to use REAXFF forcefield and for the water its simple /lj/cut/tip4p/long. The lammps data file is in sequence :
Masses
1 12.0107 # C
2 35.453 # Cl
3 1.00794 # H
4 14.0067 # N
5 22.9898 # Na
6 15.9994 # O
Accordingly the potential is defined as :
pair_style hybrid reaxff lmp_control lj/cut/tip4p/long O H 1 1 0.1546 13.0
followed by pair_coeff for all possible interactions
pair_coeff * * ffield.reax C H N O
pair_coeff {NAtype} {Ntype} lj/cut/tip4p/long 0.150000 3.0
pair_coeff {CLtype} {Ntype} lj/cut/tip4p/long 0.040000 5.0
pair_coeff {Ctype} {CLtype} lj/cut/tip4p/long 0.031702 4.2821 # C-CL
pair_coeff {Ctype} {Htype} lj/cut/tip4p/long 0 0 # C-H
pair_coeff {Ctype} {Otype} lj/cut/tip4p/long 0.12613 3.2793 # C-O
pair_coeff {Ctype} {NAtype} lj/cut/tip4p/long 0.12027 2.8293 # C-NA
pair_coeff {CLtype} {CLtype} lj/cut/tip4p/long 0.0117 5.1645 # CL-CL
pair_coeff {CLtype} {Htype} lj/cut/tip4p/long 0 0 # CL-H
pair_coeff {CLtype} {Otype} lj/cut/tip4p/long 0.046549 4.1617 # CL-O
pair_coeff {CLtype} {NAtype} lj/cut/tip4p/long 0.044388 3.7117 # CL-NA
pair_coeff {Htype} {Htype} lj/cut/tip4p/long 0 0 # H-H
pair_coeff {Htype} {Otype} lj/cut/tip4p/long 0 0 # H-O
pair_coeff {Htype} {NAtype} lj/cut/tip4p/long 0 0 # H-NA
pair_coeff {Otype} {Otype} lj/cut/tip4p/long 0.1852 3.1589 # O-O
pair_coeff {Otype} {NAtype} lj/cut/tip4p/long 0.1766 2.7089 # O-NA
pair_coeff {NAtype} {NAtype} lj/cut/tip4p/long 0.1684 2.2589 # NA-NA
kspace_style pppm/tip4p 1.0e-5
However i still get an error message as :
ERROR: Expected hybrid sub-style instead of reax/c in pair_coeff command (src/pair_hybrid.cpp:530)
Last command: pair_coeff * * reax/c ffield.reax C H N O
Since I cant upload the input file i am also sharing the complete lammps input below
atom_style full
units real
dimension 3
boundary p p p
variable Ctype equal 1
variable CLtype equal 2
variable Htype equal 3
variable Ntype equal 4
variable NAtype equal 5
variable Otype equal 6
variable seed equal 1
variable pressure equal 0.02
variable 1atm equal -0.000041
variable waterbond equal 1
variable waterangle equal 1
read_data water-c1n1.data
variable natoms equal “count(all)”
Define interaction parameters
pair_style hybrid reaxff lmp_control lj/cut/tip4p/long {Otype} {Htype} {waterbond} {waterangle} 0.1546 13.0
pair_modify tail yes
bond_style harmonic
angle_style harmonic
dihedral_style none
------------- Regions and Groups ----------------
group carbon type {Ctype} # all carbon atoms
group nitrogen type {Ntype}
group ox type {Otype}
group hy type {Htype}
group na type {NAtype}
group cl type {CLtype}
Defining the positions of all four carbon planes
variable zpiston1 equal -5
variable zpiston2 equal 71
Defining the piston and the membrane carbons
variable zmin equal {zpiston1}-1.5
variable zmax equal {zpiston1}+1.5
region piston1zone block INF INF INF INF {zmin} {zmax} units box
variable zmin equal {zpiston2}-1.5
variable zmax equal {zpiston2}+1.5
region piston2zone block INF INF INF INF {zmin} {zmax} units box
group piston1 region piston1zone # whole piston1 (carbon)
group piston2 region piston2zone # whole piston1 (carbon)
group bothpistons union piston1 piston2
defining water groups
group water union ox hy
group saltwater union water na cl
group notsaltwater subtract all saltwater
group carbon_nitride union carbon nitrogen
Defining the group of atoms that go into the Nose-Hoover thermostat
group thermostat_target union saltwater carbon_nitride
Setting water charges
set group ox charge -1.0484
set group hy charge 0.5242
set group na charge 1.
set group cl charge -1.
------------- Coefficients ----------------
pair_coeff * * ffield.reax C H N O
pair_coeff {NAtype} {Ntype} lj/cut/tip4p/long 0.150000 3.0
pair_coeff {CLtype} {Ntype} lj/cut/tip4p/long 0.040000 5.0
pair_coeff {Ctype} {CLtype} lj/cut/tip4p/long 0.031702 4.2821 # C-CL
pair_coeff {Ctype} {Htype} lj/cut/tip4p/long 0 0 # C-H
pair_coeff {Ctype} {Otype} lj/cut/tip4p/long 0.12613 3.2793 # C-O
pair_coeff {Ctype} {NAtype} lj/cut/tip4p/long 0.12027 2.8293 # C-NA
pair_coeff {CLtype} {CLtype} lj/cut/tip4p/long 0.0117 5.1645 # CL-CL
pair_coeff {CLtype} {Htype} lj/cut/tip4p/long 0 0 # CL-H
pair_coeff {CLtype} {Otype} lj/cut/tip4p/long 0.046549 4.1617 # CL-O
pair_coeff {CLtype} {NAtype} lj/cut/tip4p/long 0.044388 3.7117 # CL-NA
pair_coeff {Htype} {Htype} lj/cut/tip4p/long 0 0 # H-H
pair_coeff {Htype} {Otype} lj/cut/tip4p/long 0 0 # H-O
pair_coeff {Htype} {NAtype} lj/cut/tip4p/long 0 0 # H-NA
pair_coeff {Otype} {Otype} lj/cut/tip4p/long 0.1852 3.1589 # O-O
pair_coeff {Otype} {NAtype} lj/cut/tip4p/long 0.1766 2.7089 # O-NA
pair_coeff {NAtype} {NAtype} lj/cut/tip4p/long 0.1684 2.2589 # NA-NA
kspace_style pppm/tip4p 1.0e-5
Bond and angles coeffs
Syntax: angle_coeff N spring_constant theta_0
bond_coeff * 0.0 0.0 # Zero by default
angle_coeff * 0.0 0.0 # Zero by default
Water
bond_coeff {waterbond} 0.0 0.9572 # H2O bond (TIP4P-2005) angle_coeff {waterangle} 0.0 104.52 # H2O angle (TIP4P-2005)
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
------------- Setup ----------------
For the equilibration phase, keep the piston rigid.
fix pistonfreeze bothpistons setforce 0.0 0.0 0.0
Make sure the thermostat is looking at the relevant atoms by defining a new temperature calculation that only looks at the moveable atoms
compute selective_thermostat thermostat_target temp
And require LAMMPS to use this compute when doing anything related to temperature
thermo_modify temp selective_thermostat
------------- Minimization ----------------
thermo_style one
thermo 10
run 0
fix freezewater saltwater setforce 0.0 0.0 0.0
minimize 1.0e-6 1.0e-8 1000 1000
unfix freezewater
fix relax all box/relax x 0.0 y 0.0
minimize 1.0e-6 1.0e-8 100 1000
unfix relax
unfix pistonfreeze
------------- Equilibration ----------------
fix pistonkeep bothpistons setforce 0.0 0.0 NULL
fix piston1thrust piston1 aveforce NULL NULL 0.01 # squeezing the water between two pistons at about 150 MPa each
fix piston2thrust piston2 aveforce NULL NULL -0.01
fix 1 water shake 1.0e-4 100 0 b {waterbond} a {waterangle}
fix NVTequilib thermostat_target nvt temp 300 300 0.05
timestep 0.0005
thermo 100
dump equilibdump not_hy atom 100 equilib.lammpstrj
run 50000
write_restart monolayer.seed1.R3.p0.02.t300.dt0.5.2piston.hydrogen.tip4p.*.restart
print “Equilibration completed.”
unfix NVTequilib
undump equilibdump
------------- Dynamics --------
unfix piston1thrust
unfix piston2thrust
fix piston1push piston1 aveforce NULL NULL {pressure}
fix piston2push piston2 aveforce NULL NULL {1atm}
fix NVT thermostat_target nvt temp 300 300 0.05
thermo 100
thermo_style custom step temp etotal vol press
restart 2500000 dynamics.restart
dump fulldump all atom 1000 dynamics.lammpstrj
run 10000000
Any suggestions will be very helpful. I have checked on the forum but i still cant debug the error.
Thanks and Best Regards
Deepak Ojha