diff --git a/sample_inputs/input.in b/sample_inputs/input.in
index a7b4e8a9..315efbb1 100644
--- a/sample_inputs/input.in
+++ b/sample_inputs/input.in
@@ -31,9 +31,9 @@ nwritef=0,		! how often we print forces?
 /
 
 ! Thermostat options
-&nhcopt
+&thermostat
 temp=298.15,		! temperature [K] for Maxwell-Boltzmann sampling and thermostat
-inose=0,		! Thermostating: Nose-Hoover 1, microcanonical 0,GLE 2, Langevin 3
+therm='none',		! "none"=NVE ensemble; "nhc"=Nose-Hoover Chains thermostat, "gle"=Generalized Langevin thermostat
 nchain=4,		! number of nose-hoover chains
 tau0=0.0015,		! relaxation time of NHC thermostat
 /
diff --git a/sample_inputs/input.in.cmd b/sample_inputs/input.in.cmd
index 7a343e3b..19bfaef4 100644
--- a/sample_inputs/input.in.cmd
+++ b/sample_inputs/input.in.cmd
@@ -18,11 +18,11 @@ nwritev=0,		! how often we print velocities? (by default ABIN does not print vel
 /
 
 
-&nhcopt
+&thermostat
 temp=298.15,		! temperature [K] for Maxwell-Boltzmann sampling and thermostat
-inose=1,		! Thermostat options: 
-                        ! 0 = NVE( microcanonical)
-                        ! 1 = Nose-Hoover
-                        ! 2 = Generalized Langevin Equation (GLE)
+therm='nhc',		! Thermostat options: 
+                        ! 'none'= NVE( microcanonical)
+                        ! 'nhc' = Nose-Hoover Chains
+                        ! 'gle' = Generalized Langevin Equation (GLE)
 tau0=0.001,		! relaxation time of NHC thermostat in picoseconds
 /
diff --git a/sample_inputs/input.in.pigle b/sample_inputs/input.in.pigle
index 1d89e957..acb50877 100644
--- a/sample_inputs/input.in.pigle
+++ b/sample_inputs/input.in.pigle
@@ -17,8 +17,8 @@ nrest=5,		! how often we print restart files?
 nwritev=0,		! how often we print velocities?
 /
 
-&nhcopt
+&thermostat
 temp=298.15,		! temperature (Kelvins) for initial Maxwell-Boltzmann sampling of velocities
                         ! WARNING: The actual temperature of the simulation is determined via the GLE parameters!
-inose=2,		! GLE thermostat for QT or PI+GLE, needs additional parameters in "GLE-A" and "GLE-C" files
+therm='gle',		! GLE thermostat for QT or PI+GLE, needs additional parameters in "GLE-A" and "GLE-C" files
 /
diff --git a/sample_inputs/input.in.pimd b/sample_inputs/input.in.pimd
index 33eb364d..6a4a7974 100644
--- a/sample_inputs/input.in.pimd
+++ b/sample_inputs/input.in.pimd
@@ -16,8 +16,8 @@ nwritex=5,		! how often should we print coordinates?
 nrest=5,		! how often we print restart files?
 /
 
-&nhcopt
+&thermostat
 temp=298.15,		! temperature [K] for Maxwell-Boltzmann sampling and thermostat
-inose=1,		! NHC thermostat 
+therm='nhc',		! NHC thermostat 
 tau0=0.001,             ! thermostat relaxation time (in picoseconds)
 /
diff --git a/sample_inputs/input.in.qt b/sample_inputs/input.in.qt
index 7e2b3d02..2526c9b0 100644
--- a/sample_inputs/input.in.qt
+++ b/sample_inputs/input.in.qt
@@ -21,8 +21,8 @@ nrest=5,		! how often to print restart files?
 nwritev=0,		! how often to print velocities?
 /
 
-&nhcopt
+&thermostat
 temp=298.15,		! temperature (Kelvins) for initial Maxwell-Boltzmann sampling of velocities
                         ! WARNING: The actual temperature of the simulation is determined via the GLE parameters!
-inose=2,		! GLE thermostat for QT or PI+GLE, needs additional parameters in "GLE-A" and "GLE-C" files
+therm='gle',		! GLE thermostam for QT or PI+GLE, needs additional parameters in "GLE-A" and "GLE-C" files
 /
diff --git a/sample_inputs/input.in.sh b/sample_inputs/input.in.sh
index 548ec0b0..ac8776a7 100644
--- a/sample_inputs/input.in.sh
+++ b/sample_inputs/input.in.sh
@@ -14,8 +14,8 @@ nrest=1,		! how often to print restart files?
 nwritev=0,		! how often to print velocities?
 /
 
-&nhcopt
-inose=0,		! NVE ensemble (thermostat turned OFF)
+&thermostat
+therm='none',		! NVE ensemble (thermostat turned OFF)
 !temp=0.00,		! Usually, you would read initial velocities from previous ground state sampling (-v option)
 /
 
diff --git a/sample_inputs/input.in.shake b/sample_inputs/input.in.shake
index f1831fbe..33104462 100644
--- a/sample_inputs/input.in.shake
+++ b/sample_inputs/input.in.shake
@@ -19,12 +19,12 @@ ishake1=1,1,4,4,8,8,11,11,  ! list of bond constraints
 ishake2=2,3,5,6,9,10,12,13,
 /
 
-&nhcopt
+&thermostat
 nmolt=53,               ! number of separate systems(partitions) to couple NHC thermostat
 natmolt=53*3            ! number of atoms in each partition
 nshakemol=2,2,0,2,2,3*0,  ! number of constraints in each partition
 imasst=0,               ! DO NOT use massive thermostat when using SHAKE
-inose=1,		! Only global Nose-Hoover is available for SHAKE
+therm='nhc',		! Only global Nose-Hoover is available for SHAKE
 temp=298.15,		! initial temperature for Maxwell-Boltzmann sampling [au]
 tau0=0.001
-/
+
diff --git a/src/init.F90 b/src/init.F90
index 19378f9e..e4598fd3 100644
--- a/src/init.F90
+++ b/src/init.F90
@@ -104,6 +104,7 @@ subroutine init(dt)
       character(len=20) :: xyz_units
       character(len=60) :: chdivider
       character(len=60) :: mdtype
+      character(len=60) :: therm
       character(len=1024) :: tc_server_name, tcpb_input_file, tcpb_host
       integer :: tcpb_port
       logical :: file_exists
@@ -115,12 +116,12 @@ subroutine init(dt)
       ! in the form of the standard Fortran namelist syntax.
       ! The input parameters are grouped in several different namelists:
       !
-      ! - general: Most basic MD settings + misc
-      ! - nhcopt:  parameters for thermostats
-      ! - remd:    parameters for Replica Exchange MD
-      ! - sh:      parameters for Surface Hopping, moved to mod_sh module
-      ! - system:  system-specific parameters for model potentials, masses, SHAKE constraints...
-      ! - qmmm:    parameters for internal QMMM (not really tested).
+      ! - general:      Most basic MD settings + misc
+      ! - thermostat:   parameters for thermostats (previously nhcopt)
+      ! - remd:         parameters for Replica Exchange MD
+      ! - sh:           parameters for Surface Hopping, moved to mod_sh module
+      ! - system:       system-specific parameters for model potentials, masses, SHAKE constraints...
+      ! - qmmm:         parameters for internal QMMM (not really tested).
       !
       ! All namelists need to be in a single input file, and the code
       ! in this subroutine must ensure that the namelists can be in any order.
@@ -133,7 +134,12 @@ subroutine init(dt)
 
       namelist /remd/ nswap, nreplica, deltaT, Tmax, temp_list
 
-      namelist /nhcopt/ inose, temp, temp0, nchain, tau0, tau0_langevin, imasst, nrespnose, nyosh, &
+      ! new namelist section for thermostat, previously nhcopt
+      namelist /thermostat/ inose, therm, temp, temp0, nchain, tau0, tau0_langevin, imasst, nrespnose, nyosh, &
+         scaleveloc, readNHC, nmolt, natmolt, nshakemol, rem_comrot, rem_comvel, gle_test
+
+      ! old namelist section for thermostat, kept for backwards compatibility
+      namelist /nhcopt/ inose, therm, temp, temp0, nchain, tau0, tau0_langevin, imasst, nrespnose, nyosh, &
          scaleveloc, readNHC, nmolt, natmolt, nshakemol, rem_comrot, rem_comvel, gle_test
 
       namelist /system/ masses, massnames, ndist, dist1, dist2, &
@@ -152,6 +158,7 @@ subroutine init(dt)
       tcpb_input_file = ''
       tcpb_port = -1
       mdtype = ''
+      therm = ''
       dt = -1
       nproc = 1
       iplumed = 0
@@ -359,8 +366,14 @@ subroutine init(dt)
       natmolt(1) = natom ! default for global NHC thermostat
       nshakemol = 0
 
-      read (uinput, nhcopt)
+      ! read &thermostat section
+      read (uinput, thermostat, iostat=iost)
       rewind (uinput)
+      ! if &thermostat is not found, try the older &nhcopt
+      if (iost /= 0) then
+         read (uinput, nhcopt)
+         rewind (uinput)
+      end if
 
       if (ipimd == 2) then
          call read_sh_input(uinput)
@@ -396,6 +409,23 @@ subroutine init(dt)
          call initialize_tcpb(natqm, atnames, tcpb_port, tcpb_host, tcpb_input_file)
       end if
 
+      ! first check if 'therm' keyword was used and then transfer to inose
+      if (therm /= '') then
+         therm = tolower(therm)
+         select case (therm)
+         case ('none')
+            inose = 0
+         case ('nhc')
+            inose = 1
+         case ('gle')
+            inose = 2
+         case ('langevine')
+            inose = 3
+         case default
+            call fatal_error(__FILE__, __LINE__, 'invalid thermostat (therm) in '//trim(chinput))
+         end select
+      end if
+
       ! Check whether input parameters make sense
       call check_inputsanity()
 
diff --git a/tests/DOUBLEWELL/input.in b/tests/DOUBLEWELL/input.in
index d6976c16..d320200a 100644
--- a/tests/DOUBLEWELL/input.in
+++ b/tests/DOUBLEWELL/input.in
@@ -15,8 +15,8 @@ nwritef=1,
 idebug=1,
 /
 
-&nhcopt
-inose=0,
+&thermostat 
+therm='none', 
 temp=298.15,
 /
 
diff --git a/tests/LANGEVIN/input.in b/tests/LANGEVIN/input.in
index b05e0d78..fbff49b2 100644
--- a/tests/LANGEVIN/input.in
+++ b/tests/LANGEVIN/input.in
@@ -14,8 +14,8 @@ nwritev=1,
 nwritef=1,
 /
 
-&nhcopt
-inose=3,
+&thermostat 
+therm='langevine', 
 temp=300.0
 tau0_langevin=0.001
 /
diff --git a/tests/NHC-GLOBAL/input.in b/tests/NHC-GLOBAL/input.in
index 5dc17d44..0ab0eb38 100644
--- a/tests/NHC-GLOBAL/input.in
+++ b/tests/NHC-GLOBAL/input.in
@@ -15,8 +15,8 @@ nwritef=1,
 idebug=1,
 /
 
-&nhcopt
-inose=1,
+&thermostat 
+therm='nhc',
 imasst=0,
 nyosh=3,
 temp=298.15,
diff --git a/tests/PIGLE/input.in b/tests/PIGLE/input.in
index 59bf37bd..f85e1b61 100644
--- a/tests/PIGLE/input.in
+++ b/tests/PIGLE/input.in
@@ -17,8 +17,8 @@ nwritef=1,
 nrest=1,
 /
 
-&nhcopt
-inose=2,
+&thermostat 
+therm='gle',
 temp=200,
 rem_comrot=.false.
 gle_test=.true.