########################################################
#                   water slab system
########################################################
# water system
# This is an ordinary system of (SPCE) water which has 
# been pre-equilibrated at pressure = 1 atm, and temp = 300K.
# 
#
#    TWO WATER TYPES:
# Ordinarily, water has 2 atom types (oxygen, hydrogen).
# However there are 4 atom types in this file, corresponding to 2 types 
# of water.  The middle region of the box consists of water made from
# atoms of type 3 and 4, and is immobile.  The thinner strip above and
# below consists of water made from atoms 1 and 2 and is allowed to move.
# Atom types were assigned based on the position of the water center of mass.
# (Using different atom types makes it easier later on to customize the forces 
#  acting on the waters in the two regions...although we don't do that yet..)
#
#
#   mobile      |212  212  212  212  212 |
#   layer       |   212   212  212    212|       
#               | 212   212  212   212   |       
#            /  |434  434  434  434  434 |      Z direction
#           |   |   434  434  434  434   |      
#           |   | 434  434  434  434  434|      /|\
#           |   |   434  434  434  434   |       |
#  immobile |   |434  434  434  434  434 |       |
#  layer   <    |   434  434  434  434   |       |
#           |   |434  434  434  434  434 |       |
#           |   |  434  434  434  434 434|
#           |   |   434  434  434  434   |
#           |   |434  434   434  434  434|
#            \  |   434  434  434  434   |
#   mobile      | 212  212  212  212  212|
#   layer       |   212  212  212  212   |
#               |212  212  212  212  212 |
#
#
# -- We will immobilize atoms of type 3,4


units          real
atom_style     full
pair_style     lj/charmm/coul/long 8.0 10.0 10.0
bond_style     harmonic
angle_style    charmm
kspace_style   pppm 0.0001
boundary       p p p

read_data system.data

# Both water types are identical:

pair_coeff 1 1  0.1553 3.166 0.1553 3.166    # (standard SPCE water)
pair_coeff 2 2  0.0    1.0   0.0    1.0      # (standard SPCE water)

pair_coeff 3 3  0.1553 3.166 0.1553 3.166    # (standard SPCE water)
pair_coeff 4 4  0.0    1.0   0.0    1.0      # (standard SPCE water)

# Mixing rules are used to define interactions between unlike pairs.

pair_modify mix arithmetic


timestep     1.0

group        mobile   type 1 2
group        immobile type 3 4

#dump         1 all custom 100 traj_nvt.lammpstrj id mol type x y z ix iy iz

dump	1 all xyz 10000 x_t*.xyz
#dump	2 all dcd 10 trajectory.dcd

thermo       100


#   Formally, SHAKE is required for SPCE water.  We must be careful not 
#   to apply "fix shake" and "fix rigid" to the same atoms.
#   If this added complication bothers you, I will mention that
#   turning it off has no significant qualitative effect on the test.
#   Because the system was pre-equilibrated with shake on, I leave it on.

fix Fshake mobile shake 0.0001 10 0 b 1 a 1

####
#### COMMENTING OUT FIX RIGID
####
#####  Use fix "rigid" to make a rigid body out of the "immobile" group:
####
#velocity immobile zero angular
#velocity immobile zero linear
##### (I'm not sure if the two lines above are necessary, but they don't hurt.)
#fix Ffreezestuff immobile rigid single force * off off off torque * off off off
####


# IMPORTANT for NPT: Use "neigh_modify" to turn off calculation 
#                    of the forces between immobilized atoms.
neigh_modify exclude group immobile immobile


# ------------------- NPT -----------------------
#   The next two lines print optional diagnostics:
compute tempMobile  mobile temp
compute pressMobile all    pressure tempMobile

thermo_style custom step c_tempMobile c_pressMobile temp press vol

fix Fmovestuff mobile npt temp 300.0 300.0 100.0  z 1.0 1.0 100.0 dilate partial
#fix Fmovestuff mobile nvt temp 300.0 300.0 100.0

fix_modify Fmovestuff temp tempMobile


# --------- FOR NVT, USE THIS INSTEAD: ----------
#   The next two lines print optional diagnostics:
# compute tempMobile mobile temp
# thermo_style custom step c_tempMobile
# fix Fmovestuff mobile nvt temp 300.0 300.0 100.0


restart         20000 restart_immobile

run             100000

write_restart  system_after_fix_rigid.rst