Your input deck is cluttered with problematic choices, errors, and inconsistencies.
This is all easily understandable from the documentation. It looks like you just copied existing examples without actually studying what they do and understanding what you need.
It should be quite visible from the LAMMPS output and the created data file, that things are not going the way you expect.
Some examples:
pair_coeff 3 4 0.0000 0.0000 # O(water)-H(water)
pair_coeff 1 4 0.0000 0.0000 # Ti-H(water)
pair_coeff 2 4 0.0000 0.0000 # O(anatase)-H(water)
What would prevent the water molecules to go on top of the anatase? While the Lennard-Jones interactions are turned off, the Coulomb interactions are not and thus are strongly attractive between oppositely charged atoms. This will crash LAMMPS with “lost atoms” and crazy large negative energies and pressure.
# Bond Coefficients for SPC/E Water
bond_coeff 1 469 1.4 # O-H bond (bond type 1 in SPC/E)
bond_coeff 2 1000.0 1.0 # O-H bond (bond type 2 in SPC/E, if applicable)
angle_coeff 1 63 120
angle_coeff 2 10000.0 109.47
There should be only one bond type needed because only the water has a bond and the same for angle types.
molecule water spce.mol toff 1 boff 1 aoff 1
This offsets the atom types by one (but your anatase has two atom types), so your SPC/E oxygen becomes atom type 2 and the hydrogen atom type 3. This clearly is not what you want.
On the other hand there is no reason to use an offset for bond any angle type. This way your bond type becomes 2 instead of 1 and the same for the angle type.
region liquid block 0.0 3.785 0.0 3.785 0.0 40.0
This region extends well beyond the box dimensions in z direction (which stops a z < 10.0). If you want the extra space, you have to have a larger box, either by editing the data file or by using change_box.
group H2O type 3 4_
This will only contain the water hydrogens for the reasons explained above.
fix walls slab setforce 0.0 0.0 0.0
When turning off forces like this, you should also exclude their computation from the neighbor list with neigh_modify exclude group slab slab
. Otherwise you have bogus contributions to forces and pressure.
fix myshk H2O shake 1.0e-3 100 10000 b 1 a 1
As is shown in the log file, there are no constrained bonds or angles, because of the issues discussed above.
velocity H2O create 100.0 465135
This will assign a velocity also to the anatase oxygen atoms and fix nvt will add time integration.
This is all because of earlier mistakes.
timestep 1
When losing atoms or have missing bond atoms, the first step is always to use a more conservative time step. This is discussed almost everywhere where the above error is discussed.
The temperature printed by the thermo output is misleading, since it shows the whole system. It would be smarter to use the temperature computed by fix nvt (after correcting its group to only cover water atoms) using thermo_modify temp mynvt_temp
There are more settings that can be tweaked. As explained in many places. The way to start a simulation and to “build” an input is to do this in small steps starting from the most minimal version (i.e. only your anatase system) and then adding “complications” step by step and validating at every step that the added command (or commands) do exactly what is expected.
Many of the issues I am mentioning can be easily spotted this way.
Please note that sorting these things out is effectively your job, so don’t expect that you’ll get a similarly detailed discussion of what to improve in the future. Usually we do this with people just once to show then how many things can go wrong even when they believe they have covered everything that could be checked. The point of this forum is to get you to a state of being able to sort this out for yourself and not you dumping an input and getting it fixed by us.