#----------------------------- simulation model and settings       
log npt_polar_2.0_delete_bonds.log
units real # time in fs; temperature in K; energy in Kcal/mole; pressure in atomspheres 
dimension 3
boundary p p p
atom_style full
box tilt   large

variable dt equal 0.25
variable temp equal 300 
                  
timestep ${dt}

bond_style      harmonic  
angle_style     hybrid cosine/periodic harmonic fourier
dihedral_style  harmonic
improper_style  fourier

pair_style      lj/cut/thole/long 2.0 12.5  
pair_modify     tail yes mix arithmetic 
special_bonds   lj/coul 0.0 0.0 1.0
kspace_style    ewald 1.0e-6
neighbor        2.0 bin
neigh_modify    every 1 delay 0 check yes
read_data hydrated_apical_pristine_400_polar.data extra/special/per/atom 9
comm_modify vel yes

group gHK molecule 444

group ATOMS_mof type 1 2 3 4 5 6
group ATOMS_w type 7 8 9

group ATOMS_mof_pair type 1 2 3 4 5 6 10 11 12 13
group ATOMS_w_pair type 7 8 9 14

group gCORES type 1:9
group gDRUDES type 10:14

group DRUDES_w type 14
group DRUDES_mof type 10:13
fix DRUDE all drude C C C N C N C N N D D D D D

dump DUMP all custom 100 apical_npt_polar.lammpstrj id type x y z q
dump DUMP2 ATOMS_mof_pair custom 100 npt_mof.lammpstrj id type x y z q
dump DUMP3 ATOMS_w_pair custom 100 npt_w.lammpstrj id type x y z q

velocity        gCORES create 300. 87287 mom yes rot yes 
velocity        gDRUDES create 1.0 87287 mom yes rot yes

delete_bonds ATOMS_w multi
fix MOMENTUM all momentum 100 linear 1 1 1

compute TDRUDE all temp/drude
variable TDRUDE equal c_TDRUDE[1]

thermo_style custom step temp press vol etotal density pe ke c_TDRUDE[1] c_TDRUDE[2]
# LANGEVIN
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
######################## equilibration run
fix NVT gDRUDES nve
fix NVT_MOF ATOMS_mof nve
fix RIGID ATOMS_w rigid/nve/small molecule

thermo 100
run 40000
unfix NVT
unfix NVT_MOF
unfix RIGID