Why does special_bonds affect the intermolecular interaction?

[Jianghui]

I simulated a system with silicon dioxide and water, where a cluster of cylindrical water molecules was initially placed on the surface of silicon dioxide. The silicon dioxide was modeled using Emami et al.'s force field (https://pubs.acs.org/doi/10.1021/cm500365c), while water was modeled using the SPC/fw model. When I used the command ‘special_bonds 0 0 1’, the water molecules quickly dispersed and then gradually recondensed. However, when changed to ‘special_bonds 0 0 0’, the water molecules remained clustered. This suggests that the ‘special_bonds’ command altered the interactions between silicon dioxide and water or among the water molecules. This seems anomalous to me because ‘special_bonds’ should only affect interactions between directly or indirectly bonded atoms. Where might my understanding have gone wrong? I have attached my code below and upload the data file.
Q4_1000water.data (4.5 MB)

boundary p p p

units real
atom_style full

pair_style lj/cut/coul/long 12
bond_style harmonic
angle_style harmonic
pair_modify mix arithmetic

read_data Q4_1000water.data

group water type 5 6
group sio2 type 1 2 3 4

kspace_style pppm 1e-6

special_bonds lj/coul 0 0 1

velocity all create 298.15 85417 dist gaussian

fix 1 all nvt temp 298.15 298.15 $(100*dt)

timestep 1
thermo 1000
dump 1 all custom 1000 Q4_1000water.dump id type x y z vx vy vz
run 3000000

Thanks a lot for any help.

Jianghui

[akohlmey]

Your initial structure has rather high potential energy, so you should start your run with a minimization.

When you change the special_bonds setting, you significantly change the force field for your SiO2 slab. Even though you don’t have explicit dihedral interactions, changing special_bonds will exclude the 1-4 neighbors from the non-bonded computation. You can see the difference from the significantly higher pressure and the fact that the slab will buckle.

[Jianghui]

Thank you

After relaxing the system with ‘special_bonds 0 0 0’ before running it, I observed similar results to those with ‘special_bonds 0 0 1.’

Upon testing various scenarios, I’ve concluded that the issue doesn’t lie in the alteration of the ‘special_bonds’ setting but rather in the distinct temperature control effects among different systems. For instance, using the command:

fix 1 all nvt temp 298.15 298.15 $(100*dt)

Resulted in the silicon dioxide reaching temperatures below 200 K initially, while the water’s temperature exceeded 1,000 K.

However, by employing the following commands:

fix 1 water nvt temp 298.15 298.15 $(100*dt)

fix 2 sio2 nvt temp 298.15 298.15 $(100*dt)

The temperature control was effective, and the water molecules ceased their dispersion, aligning with my expectations. My confusion has been resolved.

Jianghui