output force on a frozen molecule

Hi all,

I’m trying to store and output the z direction force acting a frozen water molecule (the water was only froze on z direction ), The frozen water molecule was specified as group ‘freezewater’. All of the water in my system was specified as group ‘water’. The id numbers of H and O in that water were 6718, 6719 and 6720. This is the main body in my input script:

fix 2 water nvt temp 300.0 300.0 100.0
fix storeforcewater freezewater store/force
variable force1z equal f_storeforcewater[6718][3]
variable force2z equal f_storeforcewater[6719][3]
variable force3z equal f_storeforcewater[6720][3]
variable forcewater equal v_force1z+v_force2z+v_force3z

fix outputforce freezewater ave/time 5 1 5 v_forcewater file forcewater.txt

#neigh_modify exclude group freezewater moveable

fix setforcewater freezewater setforce NULL NULL 0
run 500000
unfix setforcewater
unfix outputforce
unfix storeforcewater
unfix 2

By using this input script, I always got the error: ERROR: Out of range atoms - cannot compute PPPM.

Only if I turned off the pair interaction between group ‘water’ and group ‘watera’ by activating the command line ’ neigh_modify exclude group freezewater moveable ', I can get rid of that error. However, I think the subsequently calculated force wasn’t the one that I really needed because of the excluded pair interaction.

How can I correctly calculated the z-direction force acting on the water molecule ?

Many thanks in advance,