Investigated a bit.
Basically I removed the liquid part of the system (as the poster said it can be reproduced with the solid part only), decreased the time step to 0.01fs (since sometimes this kind of problem is caused by large time steps), then remove the interactions one by one until getting good energy conservation.
The problem is in the dihedral forces. With the dihedral forces enabled I observed abrupt energy changes like this (timestep 0.01fs):
Step Temp Press Density TotEng
63045 494.6008 1086.4287 0.1646792 3839.9497
63046 494.61059 1086.4005 0.1646792 3839.9497
63047 494.62039 1086.3724 0.1646792 3839.9497
63048 494.63019 1086.3444 0.1646792 3839.9497
63049 494.64003 1086.3165 0.1646792 3839.9497
63050 494.6499 1086.2887 0.1646792 3839.9497
63051 494.65983 1086.261 0.1646792 3839.9497
63052 494.66985 1086.2335 0.1646792 3839.9497
63053 494.68001 1086.2062 0.1646792 3839.9498
63054 494.69037 1086.1791 0.1646792 3839.9498
63055 494.70105 1086.1523 0.1646792 3839.9499
63056 494.71226 1086.126 0.1646792 3839.95
63057 494.72434 1086.1003 0.1646792 3839.9502
63058 494.738 1086.0758 0.1646792 3839.9504
63059 494.75395 1086.0531 0.1646792 3839.9497
63060 494.7628 1086.026 0.1646792 3839.932
63061 494.74659 1085.9832 0.1646792 3839.9624
63062 494.72575 1085.9384 0.1646792 3839.9594
63063 494.72157 1085.9053 0.1646792 3839.9561
63064 494.72414 1085.877 0.1646792 3839.9552
63065 494.7294 1085.8505 0.1646792 3839.9548
63066 494.73592 1085.8251 0.1646792 3839.9547
63067 494.74313 1085.8002 0.1646792 3839.9546
63068 494.75075 1085.7756 0.1646792 3839.9545
63069 494.75861 1085.7513 0.1646792 3839.9545
63070 494.76663 1085.7273 0.1646792 3839.9545
63071 494.77476 1085.7033 0.1646792 3839.9545
63072 494.78295 1085.6796 0.1646792 3839.9545
63073 494.79119 1085.6559 0.1646792 3839.9545
When I replace the dihedral style with dihedral_style zero nocoeff the problem is gone.
I also tried disabling dihedral only while keeping every other force (lj, electrostatic, bond, etc.) untouched; the energy drift is <0.5kcal/mol after 100ps with a 1fs timestep (only the slab; with the liquid a smaller timestep is likely needed due to the hydrogen atoms).
Likely another case of The NVE integrator does not conserve total energy when a dihedral potential is applied to a toy polymer system . I don’t know if it can be avoided, though.
Also as a side note, you can convert a hexagonal lattice to an orthogonal lattice with a lattice transformation, which may speed up the simulation (WARNING: Triclinic box skew is large. LAMMPS will run inefficiently. (src/domain.cpp:221)).