Thermostat overwrites selective dynamics

Hello everyone,

I want to fix a specific atoms while allowing other atoms to move. However, when I use a thermostat, the restriction on the fixed atoms is lifted. On the other hand, if I do not use a thermostat, atoms behave rather unphysically, i.e., metal slabs as well as molecules tend to dissociate, even without any external perturbations.

Is there any easy remedy for one of these problems that I might be missing?

Many thanks for the help!

P.S.: I have attached a test file (for some reason I cannot attach a file anymore)

CalculationMode =td 
UnitsOutput = eV_Angstrom

%Coordinates
 "O" |  1.482769e-01*angstrom | -3.685994e-01*angstrom |  7.162566e-01*angstrom | yes
 "O" |  1.482769e-01*angstrom | -3.685994e-01*angstrom | -7.160650e-01*angstrom | yes
 "H" | -7.641828e-01*angstrom | -6.817663e-02*angstrom |  9.003930e-01*angstrom | no
 "H" |  7.641828e-01*angstrom |  3.685994e-01*angstrom | -9.003930e-01*angstrom | no
%

Dimensions = 3
ExperimentalFeatures = yes
XCFunctional = lda_x_rel + lda_c_pw
PseudopotentialSet = hgh_lda
UnitsOutput = ev_angstrom
ExtraStates = 6
ExtraStatesToConverge = 4
MoveIons = Yes

Thermostat = nose_hoover
%TDFunctions
 "temperature" | tdf_cw | 300
%

BoxShape = sphere
Radius = 4*angstrom
Spacing = 0.1*angstrom

Tf  = 20 * femtosecond
dt = 0.005 * femtosecond

TDPropagator = crank_nicolson 
TDMaxSteps = Tf/dt
TDTimeStep = dt

freq = 0.0

%TDExternalFields
  electric_field | 0 | 0 | 1 | freq | "envelope_function"
%

amp = 1
tau0 = 200*femtosecond
t0 = 105*femtosecond
tau1 = 5*femtosecond

%TDFunctions
  "envelope_function" | tdf_trapezoidal | amp | tau0 | t0 | tau1
%

%TDOutput
 geometry
%