diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/barostat_types.F src/barostat_types.F
*** src_colvar_restraint_constraint_print/barostat_types.F	2012-06-04 15:50:28.986510482 +0100
--- src/barostat_types.F	2012-06-04 15:51:28.518682699 +0100
***************
*** 21,27 ****
                                               nph_uniaxial_damped_ensemble,&
                                               nph_uniaxial_ensemble,&
                                               npt_f_ensemble,&
!                                              npt_i_ensemble
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
                                               section_vals_type,&
--- 21,28 ----
                                               nph_uniaxial_damped_ensemble,&
                                               nph_uniaxial_ensemble,&
                                               npt_f_ensemble,&
!                                              npt_i_ensemble,&
!                                              npt_y_ensemble
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
                                               section_vals_type,&
***************
*** 85,90 ****
--- 86,92 ----
      CALL section_vals_get(barostat_section, explicit=explicit, error=error)
      IF (simpar%ensemble == npt_i_ensemble .OR. &
          simpar%ensemble == npt_f_ensemble .OR. &
+         simpar%ensemble == npt_y_ensemble .OR. &
          simpar%ensemble == npe_f_ensemble .OR. &
          simpar%ensemble == npe_i_ensemble .OR. &
          simpar%ensemble == nph_uniaxial_ensemble .OR. &
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/barostat_utils.F src/barostat_utils.F
*** src_colvar_restraint_constraint_print/barostat_utils.F	2012-06-04 15:50:45.293919816 +0100
--- src/barostat_utils.F	2012-06-04 15:51:28.579160151 +0100
***************
*** 19,25 ****
                                               nph_uniaxial_damped_ensemble,&
                                               nph_uniaxial_ensemble,&
                                               npt_f_ensemble,&
!                                              npt_i_ensemble
    USE kinds,                           ONLY: default_string_length,&
                                               dp
    USE machine,                         ONLY: m_flush
--- 19,26 ----
                                               nph_uniaxial_damped_ensemble,&
                                               nph_uniaxial_ensemble,&
                                               npt_f_ensemble,&
!                                              npt_i_ensemble,&
!                                              npt_y_ensemble
    USE kinds,                           ONLY: default_string_length,&
                                               dp
    USE machine,                         ONLY: m_flush
***************
*** 80,85 ****
--- 81,89 ----
         baro_pot = -0.5_dp * v_shock * v_shock * ( 1._dp - cell%deth * iv0 )**2 - simpar%p0 * ( v0 - cell%deth )
         ! Valid only for orthorhombic cell
         baro_kin = 0.5_dp * npt ( 1, 1 )%v * npt ( 1, 1 )%v * npt ( 1, 1 )%mass
+     ELSEIF ( simpar%ensemble==npt_y_ensemble) THEN
+        baro_pot = simpar%p_ext * cell%deth
+        baro_kin = 0.5_dp * npt ( 1, 1 )%v**2 * npt ( 1, 1 )%mass
      END IF
  
    END SUBROUTINE get_baro_energies
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/constraint.F src/constraint.F
*** src_colvar_restraint_constraint_print/constraint.F	2012-06-04 15:48:43.382233527 +0100
--- src/constraint.F	2012-06-04 15:51:28.630623234 +0100
***************
*** 43,49 ****
    USE distribution_1d_types,           ONLY: distribution_1d_type
    USE f77_blas
    USE input_constants,                 ONLY: npt_f_ensemble,&
!                                              npt_i_ensemble
    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                               section_vals_type
    USE kinds,                           ONLY: default_string_length,&
--- 43,50 ----
    USE distribution_1d_types,           ONLY: distribution_1d_type
    USE f77_blas
    USE input_constants,                 ONLY: npt_f_ensemble,&
!                                              npt_i_ensemble,&
!                                              npt_y_ensemble
    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                               section_vals_type
    USE kinds,                           ONLY: default_string_length,&
***************
*** 341,346 ****
--- 342,349 ----
         CALL get_roll_matrix ('SHAKE', r_shake, v_shake, vector_r, vector_v )
      ELSE IF ( simpar % ensemble == npt_f_ensemble ) THEN
         CALL get_roll_matrix ('SHAKE', r_shake, v_shake, vector_r, vector_v, u )
+     ELSE IF ( simpar % ensemble == npt_y_ensemble ) THEN
+        CALL get_roll_matrix ('SHAKE', r_shake, v_shake, vector_r, vector_v )
      END IF
      DO k = 1, SIZE(pos,2)
         atomic_kind=>particle_set (k) % atomic_kind
***************
*** 479,484 ****
--- 482,489 ----
         CALL get_roll_matrix ('RATTLE', v_shake=r_rattle, vector_v=vector )
      ELSE IF ( simpar % ensemble == npt_f_ensemble ) THEN
         CALL get_roll_matrix ('RATTLE', v_shake=r_rattle, vector_v=vector, u=u )
+     ELSE IF ( simpar % ensemble == npt_y_ensemble ) THEN
+        CALL get_roll_matrix ('RATTLE', v_shake=r_rattle, vector_v=vector)
      END IF
      DO k = 1, SIZE(vel,2)
         atomic_kind=>particle_set (k) % atomic_kind
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/csvr_system_mapping.F src/csvr_system_mapping.F
*** src_colvar_restraint_constraint_print/csvr_system_mapping.F	2012-06-04 15:50:27.318966827 +0100
--- src/csvr_system_mapping.F	2012-06-04 15:51:28.713328955 +0100
***************
*** 20,26 ****
         do_thermo_only_master, isokin_ensemble, langevin_ensemble, &
         npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, &
         nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, nve_ensemble, &
!        nvt_ensemble, reftraj_ensemble
    USE kinds,                           ONLY: dp
    USE message_passing,                 ONLY: mp_sum
    USE molecule_kind_types,             ONLY: molecule_kind_type
--- 20,26 ----
         do_thermo_only_master, isokin_ensemble, langevin_ensemble, &
         npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, &
         nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, nve_ensemble, &
!       npt_y_ensemble, nvt_ensemble, reftraj_ensemble
    USE kinds,                           ONLY: dp
    USE message_passing,                 ONLY: mp_sum
    USE molecule_kind_types,             ONLY: molecule_kind_type
***************
*** 66,72 ****
      SELECT CASE (simpar%ensemble)
      CASE  DEFAULT
         CALL stop_program(routineN,moduleN,__LINE__,'Never reach this point!')
!     CASE ( npt_i_ensemble, npt_f_ensemble )
         map_info => csvr%map_info
         map_info%dis_type = do_thermo_only_master
  
--- 66,72 ----
      SELECT CASE (simpar%ensemble)
      CASE  DEFAULT
         CALL stop_program(routineN,moduleN,__LINE__,'Never reach this point!')
!     CASE ( npt_i_ensemble, npt_f_ensemble, npt_y_ensemble )
         map_info => csvr%map_info
         map_info%dis_type = do_thermo_only_master
  
***************
*** 132,138 ****
      CASE ( nve_ensemble, isokin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_ensemble,&
             nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble )
         CALL stop_program(routineN,moduleN,__LINE__,'Never reach this point!')
!     CASE ( nvt_ensemble, npt_i_ensemble, npt_f_ensemble)
  
         CALL setup_csvr_thermostat(csvr, thermostat_info, deg_of_freedom,&
              massive_atom_list, molecule_kind_set, local_molecules, molecule_set,&
--- 132,138 ----
      CASE ( nve_ensemble, isokin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_ensemble,&
             nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble )
         CALL stop_program(routineN,moduleN,__LINE__,'Never reach this point!')
!     CASE ( nvt_ensemble, npt_i_ensemble, npt_f_ensemble, npt_y_ensemble )
  
         CALL setup_csvr_thermostat(csvr, thermostat_info, deg_of_freedom,&
              massive_atom_list, molecule_kind_set, local_molecules, molecule_set,&
***************
*** 265,271 ****
      CASE ( isokin_ensemble, nph_uniaxial_ensemble,&
             nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble )
         CALL stop_program(routineN,moduleN,__LINE__,'Never reach this point!')
!     CASE ( nve_ensemble, npe_f_ensemble, npe_i_ensemble, nvt_ensemble, npt_i_ensemble, npt_f_ensemble )
  
         CALL setup_csvr_thermostat(csvr, thermostat_info, deg_of_freedom, massive_shell_list,&
              molecule_kind_set, local_molecules, molecule_set, para_env, nshell_local, &
--- 265,271 ----
      CASE ( isokin_ensemble, nph_uniaxial_ensemble,&
             nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble )
         CALL stop_program(routineN,moduleN,__LINE__,'Never reach this point!')
!     CASE ( nve_ensemble, npe_f_ensemble, npe_i_ensemble, nvt_ensemble, npt_i_ensemble, npt_f_ensemble, npt_y_ensemble )
  
         CALL setup_csvr_thermostat(csvr, thermostat_info, deg_of_freedom, massive_shell_list,&
              molecule_kind_set, local_molecules, molecule_set, para_env, nshell_local, &
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/extended_system_init.F src/extended_system_init.F
*** src_colvar_restraint_constraint_print/extended_system_init.F	2012-06-04 15:49:15.862631379 +0100
--- src/extended_system_init.F	2012-06-04 15:51:28.881009742 +0100
***************
*** 32,38 ****
                                               nph_uniaxial_damped_ensemble,&
                                               nph_uniaxial_ensemble,&
                                               npt_f_ensemble,&
!                                              npt_i_ensemble
    USE input_cp2k_binary_restarts,      ONLY: read_binary_thermostats_nose
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
--- 32,39 ----
                                               nph_uniaxial_damped_ensemble,&
                                               nph_uniaxial_ensemble,&
                                               npt_f_ensemble,&
!                                              npt_i_ensemble,&
!                                              npt_y_ensemble
    USE input_cp2k_binary_restarts,      ONLY: read_binary_thermostats_nose
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
***************
*** 118,123 ****
--- 119,129 ----
         CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure)
         temp = simpar % temp_baro_ext
  
+     CASE ( npt_y_ensemble )
+        ALLOCATE ( npt_info ( 1, 1 ), STAT = stat )
+        CPPrecondition(stat==0,cp_fatal_level,routineP,error,failure)
+        temp = simpar % temp_baro_ext
+ 
      CASE  DEFAULT
         ! Do nothing..
         NULLIFY(npt_info)
***************
*** 832,837 ****
--- 838,845 ----
         npt ( :, : ) % mass =  REAL ( nfree + 3,KIND=dp) * temp_ext * tau_cell ** 2 / 3.0_dp
      CASE ( npe_i_ensemble )
         npt ( :, : ) % mass = REAL ( nfree + 3,KIND=dp) * temp_ext * tau_cell ** 2
+     CASE ( npt_y_ensemble )
+        npt ( :, : ) % mass = REAL ( nfree + 3,KIND=dp) * temp_ext * tau_cell ** 2
      END SELECT
      ! initializing velocities
      DO i = 1, SIZE ( npt,1)
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/input_constants.F src/input_constants.F
*** src_colvar_restraint_constraint_print/input_constants.F	2012-06-04 15:48:50.316723174 +0100
--- src/input_constants.F	2012-06-04 15:59:53.374916605 +0100
***************
*** 108,114 ****
                                                langevin_ensemble=10,&
                                                npe_f_ensemble=11,&
                                                npe_i_ensemble=12,&
!                                               nvt_adiabatic_ensemble=13
  
    INTEGER, PARAMETER, PUBLIC               :: do_thermo_same_as_part=100,&
                                                do_thermo_nose=101,&
--- 108,115 ----
                                                langevin_ensemble=10,&
                                                npe_f_ensemble=11,&
                                                npe_i_ensemble=12,&
!                                               nvt_adiabatic_ensemble=13,&
!                                               npt_y_ensemble=14
  
    INTEGER, PARAMETER, PUBLIC               :: do_thermo_same_as_part=100,&
                                                do_thermo_nose=101,&
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/input_cp2k_check.F src/input_cp2k_check.F
*** src_colvar_restraint_constraint_print/input_cp2k_check.F	2012-06-04 15:50:44.012369232 +0100
--- src/input_cp2k_check.F	2012-06-04 15:51:29.289714026 +0100
***************
*** 22,28 ****
    USE input_constants,                 ONLY: &
         do_region_global, do_thermo_csvr, do_thermo_gle, do_thermo_nose, &
         do_thermo_same_as_part, do_vwn5, id_release_version, npt_f_ensemble, &
!        npt_i_ensemble, xc_funct_b3lyp, xc_funct_blyp, xc_funct_bp, &
         xc_funct_hcth120, xc_funct_no_shortcut, xc_funct_olyp, xc_funct_pade, &
         xc_funct_pbe, xc_funct_pbe0, xc_funct_tpss, xc_funct_xwpbe, xc_none
    USE input_cp2k,                      ONLY: create_cp2k_root_section
--- 22,28 ----
    USE input_constants,                 ONLY: &
         do_region_global, do_thermo_csvr, do_thermo_gle, do_thermo_nose, &
         do_thermo_same_as_part, do_vwn5, id_release_version, npt_f_ensemble, &
!        npt_i_ensemble, npt_y_ensemble, xc_funct_b3lyp, xc_funct_blyp, xc_funct_bp, &
         xc_funct_hcth120, xc_funct_no_shortcut, xc_funct_olyp, xc_funct_pade, &
         xc_funct_pbe, xc_funct_pbe0, xc_funct_tpss, xc_funct_xwpbe, xc_none
    USE input_cp2k,                      ONLY: create_cp2k_root_section
***************
*** 563,569 ****
            flag = check_restart(input_file, restart_file, "MOTION%MD", error)
            IF(flag) THEN
               CALL section_vals_val_get(input_file,"MOTION%MD%ENSEMBLE",i_val=ensemble,error=error)
!              IF(ensemble == npt_i_ensemble .OR. ensemble==npt_f_ensemble) THEN
                  CALL section_vals_val_get(r_section,"RESTART_BAROSTAT_THERMOSTAT",l_val=flag,error=error)
                  check = check_restart(input_file, restart_file, "MOTION%MD%BAROSTAT", error)
                  CALL restart_thermostat(flag, input_file, restart_file, "MOTION%MD%BAROSTAT%THERMOSTAT", &
--- 563,569 ----
            flag = check_restart(input_file, restart_file, "MOTION%MD", error)
            IF(flag) THEN
               CALL section_vals_val_get(input_file,"MOTION%MD%ENSEMBLE",i_val=ensemble,error=error)
!              IF(ensemble == npt_i_ensemble .OR. ensemble==npt_f_ensemble .or. ensemble==npt_y_ensemble) THEN
                  CALL section_vals_val_get(r_section,"RESTART_BAROSTAT_THERMOSTAT",l_val=flag,error=error)
                  check = check_restart(input_file, restart_file, "MOTION%MD%BAROSTAT", error)
                  CALL restart_thermostat(flag, input_file, restart_file, "MOTION%MD%BAROSTAT%THERMOSTAT", &
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/input_cp2k_motion.F src/input_cp2k_motion.F
*** src_colvar_restraint_constraint_print/input_cp2k_motion.F	2012-06-04 15:49:09.407702791 +0100
--- src/input_cp2k_motion.F	2012-06-04 16:06:29.410115474 +0100
***************
*** 1518,1530 ****
              description="The ensemble/integrator that you want to use for MD propagation",&
              usage="ensemble nve",&
              default_i_val=nve_ensemble,&
!             enum_c_vals=s2a("NVE", "NVT", "NPT_I", "NPT_F", "MSST", "MSST_DAMPED",&
              "HYDROSTATICSHOCK", "ISOKIN", "REFTRAJ", "LANGEVIN", "NPE_F", &
              "NPE_I", "NVT_ADIABATIC" ),&
              enum_desc=s2a("constant energy (microcanonical)",&
              "constant temperature and volume (canonical)",&
              "constant temperature and pressure using an isotropic cell",&
              "constant temperature and pressure using a flexible cell",&
              "simulate steady shock (uniaxial)",&
              "simulate steady shock (uniaxial) with extra viscosity",&
              "simulate steady shock with hydrostatic pressure",&
--- 1518,1531 ----
              description="The ensemble/integrator that you want to use for MD propagation",&
              usage="ensemble nve",&
              default_i_val=nve_ensemble,&
!             enum_c_vals=s2a("NVE", "NVT", "NPT_I", "NPT_F", "NPT_Y", "MSST", "MSST_DAMPED",&
              "HYDROSTATICSHOCK", "ISOKIN", "REFTRAJ", "LANGEVIN", "NPE_F", &
              "NPE_I", "NVT_ADIABATIC" ),&
              enum_desc=s2a("constant energy (microcanonical)",&
              "constant temperature and volume (canonical)",&
              "constant temperature and pressure using an isotropic cell",&
              "constant temperature and pressure using a flexible cell",&
+             "constant temperature and pressure using a flexible cell in y direction",&
              "simulate steady shock (uniaxial)",&
              "simulate steady shock (uniaxial) with extra viscosity",&
              "simulate steady shock with hydrostatic pressure",&
***************
*** 1535,1541 ****
              "constant pressure ensemble using an isotropic cell (no thermostat)",&
              "adiabatic dynamics in constant temperature and volume ensemble"),&
              citations=(/Evans1983,Minary2003,Kolafa2004/),&
!             enum_i_vals=(/nve_ensemble,nvt_ensemble,npt_i_ensemble,npt_f_ensemble,&
              nph_uniaxial_ensemble,nph_uniaxial_damped_ensemble,nph_ensemble,isokin_ensemble,&
              reftraj_ensemble,langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
              nvt_adiabatic_ensemble/),error=error)
--- 1536,1542 ----
              "constant pressure ensemble using an isotropic cell (no thermostat)",&
              "adiabatic dynamics in constant temperature and volume ensemble"),&
              citations=(/Evans1983,Minary2003,Kolafa2004/),&
!             enum_i_vals=(/nve_ensemble,nvt_ensemble,npt_i_ensemble,npt_f_ensemble,npt_y_ensemble,&
              nph_uniaxial_ensemble,nph_uniaxial_damped_ensemble,nph_ensemble,isokin_ensemble,&
              reftraj_ensemble,langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
              nvt_adiabatic_ensemble/),error=error)
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/integrator.F src/integrator.F
*** src_colvar_restraint_constraint_print/integrator.F	2012-06-04 15:49:15.616372017 +0100
--- src/integrator.F	2012-06-04 16:13:18.578890405 +0100
***************
*** 100,106 ****
  
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'integrator'
  
!   PUBLIC :: isokin, langevin, nve, nvt, npt_i, npt_f, nve_respa
    PUBLIC :: nph_uniaxial_damped, nph_uniaxial, nvt_adiabatic, reftraj
  
  
--- 100,106 ----
  
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'integrator'
  
!   PUBLIC :: isokin, langevin, nve, nvt, npt_i, npt_f, npt_y, nve_respa
    PUBLIC :: nph_uniaxial_damped, nph_uniaxial, nvt_adiabatic, reftraj
  
  
***************
*** 2667,2671 ****
--- 2667,3055 ----
  
    END SUBROUTINE nve_respa
  
+ ! *****************************************************************************
+ !> \brief npt_y integrator for particle positions & momenta
+ !>      isotropic box changes
+ !> \par History
+ !>      none
+ !> \author CJM
+ ! *****************************************************************************
+   SUBROUTINE npt_y ( md_env, globenv, error )
+ 
+     TYPE(md_environment_type), POINTER       :: md_env
+     TYPE(global_environment_type), POINTER   :: globenv
+     TYPE(cp_error_type), INTENT(inout)       :: error
+ 
+     CHARACTER(len=*), PARAMETER :: routineN = 'npt_y', &
+       routineP = moduleN//':'//routineN
+     LOGICAL, PARAMETER                       :: box_change = .TRUE.
+     REAL(KIND=dp), PARAMETER                 :: e2 = 1.0_dp/6.0_dp, &
+                                                 e4 = e2/20.0_dp, &
+                                                 e6 = e4/42.0_dp, &
+                                                 e8 = e6/72.0_dp
+ 
+     INTEGER                                  :: iroll, istat, ivar, &
+                                                 nparticle, nparticle_kind, &
+                                                 nshell
+     LOGICAL                                  :: failure, first, first_time, &
+                                                 shell_adiabatic, &
+                                                 shell_check_distance, &
+                                                 shell_present
+     REAL(KIND=dp)                            :: dt, infree, kin, roll_tol, &
+                                                 roll_tol_thrs
+     REAL(KIND=dp), DIMENSION(3)              :: vector_r, vector_v
+     REAL(KIND=dp), DIMENSION(3, 3)           :: pv_kin
+     REAL(KIND=dp), DIMENSION(:), POINTER     :: rand
+     REAL(KIND=dp), SAVE                      :: eps_0
+     TYPE(atomic_kind_list_type), POINTER     :: atomic_kinds
+     TYPE(atomic_kind_type), DIMENSION(:), &
+       POINTER                                :: atomic_kind_set
+     TYPE(cell_type), POINTER                 :: cell
+     TYPE(cp_para_env_type), POINTER          :: para_env
+     TYPE(cp_subsys_type), POINTER            :: subsys
+     TYPE(distribution_1d_type), POINTER      :: local_molecules, &
+                                                 local_particles
+     TYPE(force_env_type), POINTER            :: force_env
+     TYPE(global_constraint_type), POINTER    :: gci
+     TYPE(mol_kind_new_list_type), POINTER    :: molecule_kinds
+     TYPE(mol_new_list_type), POINTER         :: molecules
+     TYPE(molecule_kind_type), DIMENSION(:), &
+       POINTER                                :: molecule_kind_set
+     TYPE(molecule_type), DIMENSION(:), &
+       POINTER                                :: molecule_set
+     TYPE(npt_info_type), POINTER             :: npt(:,:)
+     TYPE(old_variables_type), POINTER        :: old
+     TYPE(particle_list_type), POINTER        :: core_particles, particles, &
+                                                 shell_particles
+     TYPE(particle_type), DIMENSION(:), &
+       POINTER                                :: core_particle_set, &
+                                                 particle_set, &
+                                                 shell_particle_set
+     TYPE(rng_stream_type), POINTER           :: rng_stream
+     TYPE(simpar_type), POINTER               :: simpar
+     TYPE(thermostat_type), POINTER           :: thermostat_baro, &
+                                                 thermostat_part, &
+                                                 thermostat_shell
+     TYPE(tmp_variables_type), POINTER        :: tmp
+     TYPE(virial_type), POINTER               :: virial
+ 
+     failure = .FALSE. 
+ 
+     NULLIFY(gci,thermostat_baro,thermostat_part,thermostat_shell,force_env)
+     NULLIFY(atomic_kinds,cell,para_env,subsys,local_molecules,local_particles)
+     NULLIFY(molecule_kinds,molecules,molecule_kind_set,npt)
+     NULLIFY(core_particles,particles,shell_particles, tmp, old)
+     NULLIFY(core_particle_set,particle_set,shell_particle_set)
+     NULLIFY(simpar,virial,rand)
+    
+     CALL get_md_env(md_env=md_env,simpar=simpar,force_env=force_env,&
+          thermostat_part=thermostat_part, thermostat_baro=thermostat_baro,&
+          thermostat_shell=thermostat_shell, npt=npt, first_time=first_time,&
+          para_env=para_env, error=error)
+     dt = simpar%dt
+     infree = 1.0_dp / REAL ( simpar%nfree,KIND=dp)
+ 
+     CALL force_env_get(force_env=force_env, subsys=subsys, virial=virial,&
+          cell=cell, error=error)
+ 
+     ! Do some checks on coordinates and box
+     CALL apply_qmmm_walls_reflective(force_env, error=error)
+ 
+     CALL cp_subsys_get(subsys=subsys,atomic_kinds=atomic_kinds,local_particles=local_particles,&
+          particles=particles,local_molecules_new=local_molecules,molecules_new=molecules, &
+          gci=gci,molecule_kinds_new=molecule_kinds ,error=error)
+ 
+     nparticle_kind = atomic_kinds%n_els
+     atomic_kind_set => atomic_kinds%els
+     molecule_kind_set => molecule_kinds%els
+ 
+     nparticle = particles%n_els
+     particle_set => particles%els
+     molecule_set => molecules%els
+ 
+     CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,&
+          shell_present=shell_present,shell_adiabatic=shell_adiabatic,&
+          shell_check_distance=shell_check_distance)
+ 
+     IF ( first_time ) THEN
+        CALL virial_evaluate ( atomic_kind_set, particle_set,  &
+             local_particles, virial, para_env%group, error=error)
+     END IF
+ 
+     ! Allocate work storage for positions and velocities
+     CALL allocate_old ( old, particle_set, npt, error=error )
+ 
+     IF (ASSOCIATED(force_env%meta_env)) THEN
+ print *,'meta_env is not implemented with npt_y.'
+ stop
+        ! Allocate random number for Langevin Thermostat acting on COLVARS
+        IF(force_env%meta_env%langevin) THEN
+           ALLOCATE (rand(force_env%meta_env%n_colvar),STAT=istat)
+           CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
+           rand(:) = 0.0_dp
+        ENDIF
+     ENDIF
+ 
+     IF(shell_present) THEN
+ print *,'shell is not implemented with npt_y.'
+ stop
+        CALL cp_subsys_get(subsys=subsys,&
+             shell_particles=shell_particles, core_particles=core_particles,&
+             error=error)
+        shell_particle_set =>  shell_particles%els
+        nshell = SIZE(shell_particles%els)
+        IF(shell_adiabatic) THEN
+           core_particle_set =>  core_particles%els
+        END IF
+     END IF
+ 
+     CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic, error=error)
+ 
+     ! Initialize eps_0 the first time through
+     IF ( first_time ) eps_0 = npt (1,1)%eps
+ if (first_time) print *,'npt_y first eps_0 ',eps_0
+ 
+     ! Apply thermostat to  barostat
+     ! rescale barostat kinetic energy according to thermostat checked csvr barostat, not lnhc
+     CALL apply_thermostat_baro( thermostat_baro, npt, para_env%group, error)
+ 
+     ! Apply Thermostat over the full set of particles
+     IF(simpar% ensemble /= npe_i_ensemble) THEN
+       IF(shell_adiabatic) THEN
+ print *,'shell is not implemented with npt_y.'
+ stop
+         CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,&
+              particle_set, local_molecules, para_env%group,shell_adiabatic=shell_adiabatic,&
+              shell_particle_set=shell_particle_set, core_particle_set=core_particle_set,&
+              error=error)
+ 
+       ELSE
+         CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,&
+              particle_set, local_molecules, para_env%group ,error=error)
+       END IF
+     END IF
+ 
+     ! Apply Thermostat over the core-shell motion 
+     CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set,&
+          local_particles, para_env%group, shell_particle_set=shell_particle_set,&
+          core_particle_set=core_particle_set, error=error)
+ 
+     IF (simpar%constraint) THEN
+        ! Possibly update the target values
+        CALL shake_update_targets(gci, local_molecules, molecule_set, &
+             molecule_kind_set, dt, force_env%root_section, error)
+     END IF
+ 
+     ! setting up for ROLL: saving old variables
+     IF (simpar%constraint) THEN
+        roll_tol_thrs = simpar%roll_tol
+        iroll = 1
+        CALL set ( old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F' )
+        CALL getold ( gci, local_molecules, molecule_set,  &
+             molecule_kind_set, particle_set, cell, error)
+     ELSE
+        roll_tol_thrs = EPSILON(0.0_dp)
+     ENDIF
+     roll_tol = -roll_tol_thrs
+ 
+     !    *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2)
+     IF (ASSOCIATED(force_env%meta_env)) THEN
+ print *,'meta_env is not implemented with npt_y.'
+ stop
+        IF(force_env%meta_env%langevin) THEN
+           DO ivar = 1 , force_env%meta_env%n_colvar
+              rng_stream => force_env%meta_env%rng(ivar)%stream
+              rand(ivar)=next_random_number(rng_stream,error=error)
+           ENDDO
+           CALL metadyn_velocities_colvar(force_env,rand,error=error)
+        ENDIF
+     ENDIF
+ 
+     SR: DO WHILE (ABS(roll_tol)>=roll_tol_thrs) ! SHAKE-ROLL LOOP
+ 
+        IF (simpar%constraint) THEN
+           CALL set ( old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B' )
+        END IF
+ 
+        CALL update_pv ( gci, simpar, atomic_kind_set, particle_set, &
+             local_molecules, molecule_set, molecule_kind_set, &
+             local_particles, kin, pv_kin, virial, para_env%group )
+        CALL update_veps ( cell, npt, simpar, pv_kin, kin, virial, infree )
+ 
+        !POS scaling -- special case of npt_f, similarly to npt_i or nph_uniaxial
+        tmp%arg_r(1) = ( 0.5_dp * npt ( 1, 1 )%v * dt ) * &
+             ( 0.5_dp * npt ( 1, 1 )%v * dt )
+ 
+        tmp%poly_r(2) = 1.0_dp + e2*tmp%arg_r(1) + e4*tmp%arg_r(1)*tmp%arg_r(1) + &
+             e6*tmp%arg_r(1)**3 + e8*tmp%arg_r(1)**4
+        tmp%poly_r(1) = 1.0_dp
+        tmp%poly_r(3) = 1.0_dp
+ 
+        tmp%scale_r(1) = 1.0_dp
+        tmp%scale_r(2) = EXP( 0.5_dp * dt * npt ( 1, 1 )%v )
+        tmp%scale_r(3) = 1.0_dp
+ 
+        !VEL scaling -- special case of npt_f, similarly to npt_i or nph_uniaxial
+        tmp%arg_v(1) = ( 0.25_dp * npt ( 1, 1 )%v * dt *  &
+             (1.0_dp + infree ) )*( 0.25_dp * npt ( 1, 1 )%v * &
+             dt * ( 1.0_dp + infree ) )
+        tmp%arg_v(2) = ( 0.25_dp * npt ( 1, 1 )%v * dt * infree) * &
+             ( 0.25_dp * npt ( 1, 1 )%v * dt * infree )
+ 
+        tmp%poly_v(2) = 1.0_dp + e2*tmp%arg_v(1) + e4*tmp%arg_v(1)*tmp%arg_v(1) + &
+             e6*tmp%arg_v(1)**3 + e8*tmp%arg_v(1)**4
+        tmp%poly_v(1) = 1.0_dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + &
+             e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4
+        tmp%poly_v(3) = 1.0_dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + &
+             e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4
+ 
+        !this would be from nph_uniaxial
+        tmp%scale_v(2) = EXP( -0.25_dp * dt * npt ( 1, 1 )%v * &
+             (1.0_dp + infree ) )
+        tmp%scale_v(1) = EXP( -0.25_dp * dt * npt ( 1, 1 )%v * infree )
+        tmp%scale_v(3) = EXP( -0.25_dp * dt * npt ( 1, 1 )%v * infree )
+ 
+        ! first half of velocity verlet
+        CALL vv_first( tmp, atomic_kind_set, local_particles, particle_set, &
+             core_particle_set, shell_particle_set, nparticle_kind,&
+             shell_adiabatic, dt, error=error )
+ 
+        IF(simpar%variable_dt) CALL  variable_timestep(md_env,tmp, dt, simpar,para_env,&
+             atomic_kind_set, local_particles, particle_set, core_particle_set,&
+             shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt,error=error)
+ 
+        roll_tol = 0.0_dp
+        vector_r ( : ) = 0._dp !scale_r * poly_r would be 1._dp
+        vector_r ( 2 ) = tmp%scale_r(2) * tmp%poly_r(2)
+        vector_v ( : ) = tmp%scale_v(:) * tmp%poly_v(:)
+ 
+        IF ( simpar%constraint ) then
+           CALL shake_roll_control( gci, local_molecules, &
+             molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar,  &
+             roll_tol, iroll, vector_r, vector_v, para_env%group, cell=cell, &
+             local_particles=local_particles, error=error )
+        endif
+     END DO SR
+ 
+     ! Update eps:
+     npt ( :, : )%eps =  npt ( :, : )%eps + dt * npt ( :, : )%v
+ 
+     ! Update h_mat
+     cell%hmat ( 2, 2 ) = cell%hmat ( 2, 2 ) * EXP ( npt ( 1, 1 )%eps - eps_0 )
+ 
+     eps_0 =  npt ( 1, 1 )%eps
+ 
+     ! Update the inverse
+     CALL init_cell ( cell )
+ 
+     ! Broadcast the new particle positions and deallocate the pos components of temporary
+     CALL update_dealloc_tmp ( tmp, particle_set, shell_particle_set, &
+          core_particle_set, para_env, shell_adiabatic, pos=.TRUE.,  error=error )
+ 
+     IF(shell_adiabatic .AND. shell_check_distance) THEN
+         CALL optimize_shell_core(force_env,  particle_set,&
+              shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE., error=error)
+     END IF
+ 
+     ! Update forces
+     CALL force_env_set_cell(force_env, cell, error)
+     CALL force_env_calc_energy_force(force_env,error=error)
+ 
+     ! Metadynamics
+     CALL metadyn_integrator(force_env, tmp%vel, rand=rand, error=error)
+ 
+     ! Velocity Verlet (second part)
+     CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
+          core_particle_set, shell_particle_set, nparticle_kind,&
+          shell_adiabatic, dt, error=error )
+ 
+     IF (simpar%constraint) THEN
+        roll_tol_thrs = simpar%roll_tol
+        first = .TRUE.
+        iroll = 1
+        CALL set ( old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F' )
+     ELSE
+        roll_tol_thrs = EPSILON(0.0_dp)
+     ENDIF
+     roll_tol = -roll_tol_thrs
+ 
+     RR: DO WHILE (ABS(roll_tol)>=roll_tol_thrs) ! RATTLE-ROLL LOOP
+        roll_tol = 0.0_dp
+        IF ( simpar%constraint ) CALL rattle_roll_setup ( old, gci, atomic_kind_set, &
+             particle_set, local_particles, molecule_kind_set, molecule_set,  &
+             local_molecules, tmp%vel, dt,  cell, npt, simpar, virial, vector_v,  &
+             roll_tol, iroll, infree, first, para_env,error=error)
+ 
+        CALL update_pv ( gci, simpar, atomic_kind_set, tmp%vel, particle_set, &
+             local_molecules, molecule_set, molecule_kind_set, &
+             local_particles, kin, pv_kin, virial, para_env%group )
+        CALL update_veps ( cell, npt, simpar, pv_kin, kin, virial, infree )
+     END DO RR
+ 
+     ! Apply Thermostat over the full set of particles
+     IF(simpar% ensemble /= npe_i_ensemble) THEN
+       IF (shell_adiabatic) THEN
+         CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, &
+              particle_set, local_molecules, para_env%group, shell_adiabatic=shell_adiabatic,&
+              vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel, error=error)
+       ELSE
+         CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, &
+              particle_set, local_molecules, para_env%group, vel=tmp%vel, error=error)
+       END IF
+     END IF
+ 
+ 
+     ! Apply Thermostat over the core-shell motion 
+     IF(ASSOCIATED(thermostat_shell)) THEN
+       CALL apply_thermostat_shells(thermostat_shell,atomic_kind_set, particle_set,&
+            local_particles, para_env%group, vel=tmp%vel, shell_vel=tmp%shell_vel, &
+            core_vel=tmp%core_vel, error=error)
+     END IF
+ 
+     ! Apply Thermostat to Barostat
+     CALL apply_thermostat_baro( thermostat_baro, npt, para_env%group, error)
+ 
+     ! Annealing of particle velocities is only possible when no thermostat is active
+     IF (simpar% ensemble == npe_i_ensemble .AND. simpar%annealing) THEN
+        tmp%vel(:,:)=tmp%vel(:,:)*simpar%f_annealing
+        IF (shell_adiabatic) THEN
+           CALL shell_scale_comv(atomic_kind_set,local_particles,particle_set,&
+                tmp%vel,tmp%shell_vel,tmp%core_vel)
+        END IF
+     END IF
+     ! Annealing of CELL velocities is only possible when no thermostat is active
+     IF (simpar% ensemble == npe_i_ensemble .AND. simpar%annealing_cell) THEN
+        npt(1,1)%v = npt(1,1)%v * simpar%f_annealing_cell
+     END IF
+ 
+ 
+     ! Broadcast the new particle velocities and deallocate temporary
+     CALL update_dealloc_tmp ( tmp, particle_set, shell_particle_set, &
+          core_particle_set, para_env, shell_adiabatic, vel=.TRUE.,  error=error )
+ 
+     ! Update constraint virial
+     IF ( simpar%constraint ) CALL pv_constraint ( gci, local_molecules, &
+          molecule_set, molecule_kind_set, particle_set, virial, para_env%group )
+ 
+     CALL virial_evaluate ( atomic_kind_set, particle_set,  &
+          local_particles, virial, para_env%group , error=error)
+ 
+     ! Deallocate old variables
+     CALL deallocate_old (old, error=error )
+ 
+     IF (ASSOCIATED(force_env%meta_env)) THEN
+        IF(force_env%meta_env%langevin) THEN
+           DEALLOCATE(rand,STAT=istat)    
+           CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
+        ENDIF
+     ENDIF
+ 
+     IF (first_time) THEN
+        first_time = .FALSE.
+        CALL set_md_env(md_env, first_time=first_time,error=error)
+     END IF
+ 
+   END SUBROUTINE npt_y
+ 
  END MODULE integrator
  
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/integrator_utils.F src/integrator_utils.F
*** src_colvar_restraint_constraint_print/integrator_utils.F	2012-06-04 15:50:58.408028158 +0100
--- src/integrator_utils.F	2012-06-04 15:51:29.903193206 +0100
***************
*** 26,32 ****
         do_clv_x, do_clv_xy, do_clv_xyz, do_clv_xz, do_clv_y, do_clv_yz, &
         do_clv_z, isokin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
!        npt_i_ensemble, nve_ensemble, nvt_ensemble
    USE kinds,                           ONLY: dp
    USE mathlib,                         ONLY: matvec_3x3,&
                                               transpose_3d
--- 26,32 ----
         do_clv_x, do_clv_xy, do_clv_xyz, do_clv_xz, do_clv_y, do_clv_yz, &
         do_clv_z, isokin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
!        npt_i_ensemble, nve_ensemble, nvt_ensemble, npt_y_ensemble
    USE kinds,                           ONLY: dp
    USE mathlib,                         ONLY: matvec_3x3,&
                                               transpose_3d
***************
*** 750,755 ****
--- 750,762 ----
         npt_loc ( 3, 3 ) % mass = npt ( 1, 1 ) % mass
      CASE ( npt_f_ensemble )
         npt_loc = npt
+     CASE ( npt_y_ensemble)
+ !       print *,'Constraints with npt_y_ensemble are not implemented'
+ !       stop
+        npt_loc ( :, : ) % v = 0.0_dp
+        npt_loc ( :, : ) % mass = 0.0_dp
+        npt_loc ( 2, 2 ) % v = npt ( 1, 1 ) % v
+        npt_loc ( 2, 2 ) % mass = npt ( 1, 1 ) % mass
      END SELECT
  
      ! resetting
***************
*** 960,965 ****
--- 967,978 ----
                 pv_kin ( 1, 1 ) +  virial % pv_constraint ( 1, 1 ) -  &
                 simpar % p0 * box % deth + infree * kin
         END IF
+     ELSEIF ( simpar % ensemble == npt_y_ensemble ) THEN
+        ! get force on barostat, as a special case of npt_f
+        npt ( 1, 1 ) % f = virial % pv_virial ( 2, 2 ) +  &
+             pv_kin ( 2, 2 ) +  virial % pv_constraint ( 2, 2 ) -  &
+             simpar % p_ext * box % deth +  &
+             infree * kin
      ENDIF
  
      ! update barostat velocities
***************
*** 1668,1673 ****
--- 1681,1689 ----
              core_particle_set, shell_particle_set, nparticle_kind,&
              shell_adiabatic, dt, error=error)
  
+     CASE (npt_y_ensemble)
+         print *,'Rescaled time step with npt_y_ensemble is not implemented'
+         stop
      END SELECT
  
    END SUBROUTINE rescaled_vv_first
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/md_conserved_quantities.F src/md_conserved_quantities.F
*** src_colvar_restraint_constraint_print/md_conserved_quantities.F	2012-06-04 15:50:35.478583359 +0100
--- src/md_conserved_quantities.F	2012-06-04 16:17:28.355311019 +0100
***************
*** 26,32 ****
         isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
         npt_i_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, &
!        reftraj_ensemble
    USE kinds,                           ONLY: dp
    USE mathconstants,                   ONLY: zero
    USE md_ener_types,                   ONLY: md_ener_type,&
--- 26,32 ----
         isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
         npt_i_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, &
!        reftraj_ensemble, npt_y_ensemble
    USE kinds,                           ONLY: dp
    USE mathconstants,                   ONLY: zero
    USE md_ener_types,                   ONLY: md_ener_type,&
***************
*** 150,156 ****
         CALL get_econs_nve (md_env, md_ener,  para_env,error=error)
      CASE ( nvt_ensemble)
         CALL get_econs_nvt (md_env, md_ener,  para_env,error=error)
!     CASE (  npt_i_ensemble,  npt_f_ensemble )
         CALL get_econs_npt (md_env, md_ener,  para_env, error=error)
         md_ener%temp_baro = md_ener%temp_baro*kelvin
      CASE (  nph_uniaxial_ensemble )
--- 150,156 ----
         CALL get_econs_nve (md_env, md_ener,  para_env,error=error)
      CASE ( nvt_ensemble)
         CALL get_econs_nvt (md_env, md_ener,  para_env,error=error)
!     CASE (  npt_i_ensemble,  npt_f_ensemble,  npt_y_ensemble )
         CALL get_econs_npt (md_env, md_ener,  para_env, error=error)
         md_ener%temp_baro = md_ener%temp_baro*kelvin
      CASE (  nph_uniaxial_ensemble )
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/md_energies.F src/md_energies.F
*** src_colvar_restraint_constraint_print/md_energies.F	2012-06-04 15:51:02.033330966 +0100
--- src/md_energies.F	2012-06-04 15:51:30.599774674 +0100
***************
*** 39,44 ****
--- 39,45 ----
                                               nph_uniaxial_ensemble,&
                                               npt_f_ensemble,&
                                               npt_i_ensemble,&
+                                              npt_y_ensemble,&
                                               reftraj_ensemble
    USE input_cp2k_motion,               ONLY: create_md_section
    USE input_enumeration_types,         ONLY: enum_i2c,&
***************
*** 752,757 ****
--- 753,759 ----
            IF ( simpar%ensemble == nph_uniaxial_ensemble .OR.&
                 simpar%ensemble == nph_uniaxial_damped_ensemble.OR.&
                 simpar%ensemble == npt_i_ensemble.OR.&
+                simpar%ensemble == npt_y_ensemble.OR.&
                 simpar%ensemble == npt_f_ensemble.OR.&
                 simpar%ensemble == npe_i_ensemble.OR.&
                 simpar%ensemble == npe_f_ensemble)  &
***************
*** 831,836 ****
--- 833,839 ----
            IF ( simpar%ensemble == nph_uniaxial_ensemble .OR.&
                 simpar%ensemble == nph_uniaxial_damped_ensemble.OR.&
                 simpar%ensemble == npt_i_ensemble.OR.&
+                simpar%ensemble == npt_y_ensemble.OR.&
                 simpar%ensemble == npt_f_ensemble.OR.&
                 simpar%ensemble == npe_i_ensemble.OR.&
                 simpar%ensemble == npe_f_ensemble) THEN
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/md_run.F src/md_run.F
*** src_colvar_restraint_constraint_print/md_run.F	2012-06-04 15:48:45.676522739 +0100
--- src/md_run.F	2012-06-04 15:51:30.691770704 +0100
***************
*** 34,45 ****
                                               force_env_type
    USE free_energy_methods,             ONLY: free_energy_evaluate
    USE free_energy_types,               ONLY: fe_env_create,&
                                               free_energy_type
    USE global_types,                    ONLY: global_environment_type
    USE input_constants,                 ONLY: &
         ehrenfest, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
!        npt_i_ensemble, reftraj_ensemble
    USE input_cp2k_check,                ONLY: remove_restart_info
    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                               section_vals_remove_values,&
--- 34,46 ----
                                               force_env_type
    USE free_energy_methods,             ONLY: free_energy_evaluate
    USE free_energy_types,               ONLY: fe_env_create,&
+                                              fe_env_release, &
                                               free_energy_type
    USE global_types,                    ONLY: global_environment_type
    USE input_constants,                 ONLY: &
         ehrenfest, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
!        npt_i_ensemble, npt_y_ensemble, reftraj_ensemble
    USE input_cp2k_check,                ONLY: remove_restart_info
    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                               section_vals_remove_values,&
***************
*** 287,292 ****
--- 288,295 ----
      CALL set_md_env(md_env=md_env, simpar=simpar, fe_env=fe_env, cell=cell,&
                      force_env=force_env, error=error)
  
+ !call fe_env_release(fe_env, error=error)
+ 
      ! Possibly initialize Wiener processes
      IF (simpar%ensemble == langevin_ensemble) CALL create_wiener_process(md_env,error)
      time_iter_start=m_walltime()
***************
*** 304,309 ****
--- 307,313 ----
      ! multiple force_evals
      IF ( (simpar%ensemble==npt_i_ensemble).OR.&
           (simpar%ensemble==npt_f_ensemble).OR.&
+          (simpar%ensemble==npt_y_ensemble).OR.&
           (simpar%ensemble==npe_f_ensemble).OR.&
           (simpar%ensemble==npe_i_ensemble).OR.&
           (simpar%ensemble==nph_uniaxial_ensemble).OR.&
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/md_vel_utils.F src/md_vel_utils.F
*** src_colvar_restraint_constraint_print/md_vel_utils.F	2012-06-04 15:51:00.685814130 +0100
--- src/md_vel_utils.F	2012-06-04 15:51:30.755776140 +0100
***************
*** 36,42 ****
    USE global_types,                    ONLY: global_environment_type
    USE input_constants,                 ONLY: &
         npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, &
!        nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
         reftraj_ensemble, use_perd_none, use_perd_x, use_perd_xy, &
         use_perd_xyz, use_perd_xz, use_perd_y, use_perd_yz, use_perd_z
    USE input_cp2k_binary_restarts,      ONLY: read_binary_velocities
--- 36,42 ----
    USE global_types,                    ONLY: global_environment_type
    USE input_constants,                 ONLY: &
         npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, &
!        nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, npt_y_ensemble, &
         reftraj_ensemble, use_perd_none, use_perd_x, use_perd_xy, &
         use_perd_xyz, use_perd_xz, use_perd_y, use_perd_yz, use_perd_z
    USE input_cp2k_binary_restarts,      ONLY: read_binary_velocities
***************
*** 1262,1267 ****
--- 1262,1270 ----
        IF ( simpar%ensemble==npt_i_ensemble .OR. simpar%ensemble==npe_i_ensemble) THEN
           npt ( 1, 1 )%v = npt ( 1, 1 )%v*scale
           md_ener%baro_kin = 0.5_dp * npt ( 1, 1 )%v**2 * npt ( 1, 1 )%mass
+        ELSE IF (simpar%ensemble==npt_y_ensemble) THEN
+          npt ( 1, 1 )%v = npt ( 1, 1 )%v*scale
+          md_ener%baro_kin = 0.5_dp * npt ( 1, 1 )%v**2 * npt ( 1, 1 )%mass
         ELSE IF (simpar%ensemble==npt_f_ensemble .OR. simpar%ensemble==npe_f_ensemble) THEN
           md_ener%baro_kin = 0.0_dp
           DO i = 1, 3
***************
*** 1330,1336 ****
      ! Control cell motion
      SELECT CASE (simpar%ensemble)
      CASE( nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
!             npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble)
         IF(simpar%temp_baro_tol > 0.0_dp) THEN
            CALL scale_velocity_baro(md_env, md_ener, simpar%temp_baro_ext, simpar%temp_baro_tol, iw, error)
         END IF
--- 1333,1339 ----
      ! Control cell motion
      SELECT CASE (simpar%ensemble)
      CASE( nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
!             npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble, npt_y_ensemble)
         IF(simpar%temp_baro_tol > 0.0_dp) THEN
            CALL scale_velocity_baro(md_env, md_ener, simpar%temp_baro_ext, simpar%temp_baro_tol, iw, error)
         END IF
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/simpar_methods.F src/simpar_methods.F
*** src_colvar_restraint_constraint_print/simpar_methods.F	2012-06-04 15:48:43.939548941 +0100
--- src/simpar_methods.F	2012-06-04 15:51:30.819432574 +0100
***************
*** 23,29 ****
    USE input_constants,                 ONLY: &
         isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
!        npt_i_ensemble, nvt_ensemble, reftraj_ensemble
    USE input_cp2k_motion,               ONLY: create_md_section
    USE input_enumeration_types,         ONLY: enum_i2c,&
                                               enumeration_type
--- 23,29 ----
    USE input_constants,                 ONLY: &
         isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
!        npt_i_ensemble, npt_y_ensemble, nvt_ensemble, reftraj_ensemble
    USE input_cp2k_motion,               ONLY: create_md_section
    USE input_enumeration_types,         ONLY: enum_i2c,&
                                               enumeration_type
***************
*** 133,138 ****
--- 133,139 ----
              simpar%f_annealing_cell
         IF (  simpar%ensemble == npt_i_ensemble .OR. &
               simpar%ensemble == npt_f_ensemble .OR. &
+              simpar%ensemble == npt_y_ensemble .OR. &
               simpar%ensemble == npe_i_ensemble .OR. &
               simpar%ensemble == npe_f_ensemble ) THEN
            tmp_r1 = cp_unit_from_cp2k(simpar%p_ext,"bar",error=error)
***************
*** 169,174 ****
--- 170,176 ----
              simpar%ensemble == nph_uniaxial_damped_ensemble.OR.&
              simpar%ensemble == npt_i_ensemble.OR.&
              simpar%ensemble == npt_f_ensemble .OR. &
+             simpar%ensemble == npt_y_ensemble .OR. &
              simpar%ensemble == npe_i_ensemble .OR. &
              simpar%ensemble == npe_f_ensemble) THEN
  
***************
*** 229,238 ****
         IF(simpar%annealing ) THEN
           IF ( simpar%ensemble == nvt_ensemble .OR. &
                simpar%ensemble == npt_i_ensemble .OR. &
                simpar%ensemble == npt_f_ensemble )  THEN
                CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,&
                "Annealing of the ions has been required "//&
!               "even if the thermostat is active (nvt or npt_i or npt_f) "//&
                "These two methods to control the temperature act one against the other. "//&
  CPSourceFileRef,&
                 only_ionode=.TRUE.)
--- 231,241 ----
         IF(simpar%annealing ) THEN
           IF ( simpar%ensemble == nvt_ensemble .OR. &
                simpar%ensemble == npt_i_ensemble .OR. &
+               simpar%ensemble == npt_y_ensemble .OR. &
                simpar%ensemble == npt_f_ensemble )  THEN
                CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,&
                "Annealing of the ions has been required "//&
!               "even if the thermostat is active (nvt or npt_i or npt_f or npt_y) "//&
                "These two methods to control the temperature act one against the other. "//&
  CPSourceFileRef,&
                 only_ionode=.TRUE.)
***************
*** 332,338 ****
  
         SELECT CASE (simpar%ensemble)
         CASE( nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
!                npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble)
            tmp_section => section_vals_get_subs_vals(md_section,"BAROSTAT",error=error)
            CALL section_vals_val_get(tmp_section,"PRESSURE",r_val=simpar%p_ext,error=error)
            CALL section_vals_val_get(tmp_section,"TIMECON",r_val=simpar%tau_cell,error=error)
--- 335,341 ----
  
         SELECT CASE (simpar%ensemble)
         CASE( nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
!                npt_f_ensemble, npt_i_ensemble, npt_y_ensemble, npe_f_ensemble, npe_i_ensemble)
            tmp_section => section_vals_get_subs_vals(md_section,"BAROSTAT",error=error)
            CALL section_vals_val_get(tmp_section,"PRESSURE",r_val=simpar%p_ext,error=error)
            CALL section_vals_val_get(tmp_section,"TIMECON",r_val=simpar%tau_cell,error=error)
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/thermal_region_utils.F src/thermal_region_utils.F
*** src_colvar_restraint_constraint_print/thermal_region_utils.F	2012-06-04 15:50:28.050834349 +0100
--- src/thermal_region_utils.F	2012-06-04 15:51:30.859732935 +0100
***************
*** 19,24 ****
--- 19,25 ----
                                               force_env_type
    USE input_constants,                 ONLY: npt_f_ensemble,&
                                               npt_i_ensemble,&
+                                              npt_y_ensemble,&
                                               nvt_ensemble
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
***************
*** 79,84 ****
--- 80,86 ----
      IF (explicit)THEN
         apply_thermostat =(simpar%ensemble == nvt_ensemble)  .OR.&
                           (simpar%ensemble == npt_f_ensemble).OR.&
+                          (simpar%ensemble == npt_y_ensemble).OR.&
                           (simpar%ensemble == npt_i_ensemble)
         IF(apply_thermostat) THEN
           CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,&
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/thermostat_methods.F src/thermostat_methods.F
*** src_colvar_restraint_constraint_print/thermostat_methods.F	2012-06-04 15:49:15.529675783 +0100
--- src/thermostat_methods.F	2012-06-04 16:24:50.148507288 +0100
***************
*** 44,50 ****
         do_region_global, do_thermo_csvr, do_thermo_gle, do_thermo_nose, &
         do_thermo_same_as_part, npe_f_ensemble, npe_i_ensemble, &
         npt_f_ensemble, npt_i_ensemble, nve_ensemble, nvt_adiabatic_ensemble, &
!        nvt_ensemble
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
                                               section_vals_remove_values,&
--- 44,50 ----
         do_region_global, do_thermo_csvr, do_thermo_gle, do_thermo_nose, &
         do_thermo_same_as_part, npe_f_ensemble, npe_i_ensemble, &
         npt_f_ensemble, npt_i_ensemble, nve_ensemble, nvt_adiabatic_ensemble, &
!        nvt_ensemble, npt_y_ensemble
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
                                               section_vals_remove_values,&
***************
*** 154,159 ****
--- 154,160 ----
      apply_thermo_adiabatic = (simpar%ensemble==nvt_adiabatic_ensemble)
  
      apply_thermo_baro = (simpar%ensemble == npt_f_ensemble).OR.&
+                         (simpar%ensemble == npt_y_ensemble).OR.&
                          (simpar%ensemble == npt_i_ensemble).AND.&
                          (.NOT. apply_thermo_adiabatic )
  
***************
*** 164,169 ****
--- 165,171 ----
                           (simpar%ensemble == nvt_ensemble)  .OR.&
                           (simpar%ensemble == npt_f_ensemble).OR.&
                           (simpar%ensemble == npt_i_ensemble).OR.&
+                          (simpar%ensemble == npt_y_ensemble).OR.&
                           (simpar%ensemble == npe_i_ensemble).OR.&
                           (simpar%ensemble == npe_f_ensemble).AND.&
                           (.NOT. apply_thermo_adiabatic )
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/thermostat_utils.F src/thermostat_utils.F
*** src_colvar_restraint_constraint_print/thermostat_utils.F	2012-06-04 15:50:53.405514273 +0100
--- src/thermostat_utils.F	2012-06-04 16:27:54.957661681 +0100
***************
*** 30,36 ****
         langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
         npt_i_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, &
!        reftraj_ensemble
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
                                               section_vals_type,&
--- 30,36 ----
         langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
         npt_i_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, &
!        reftraj_ensemble, npt_y_ensemble
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
                                               section_vals_type,&
***************
*** 574,580 ****
           reftraj_ensemble, langevin_ensemble)
         ! Do Nothing
      CASE(nve_ensemble, nvt_ensemble, nvt_adiabatic_ensemble, npt_i_ensemble, &
!         npt_f_ensemble, npe_i_ensemble, npe_f_ensemble)
         IF (ensemble==nve_ensemble) check = do_shell
         IF (check) THEN
            SELECT CASE(region)
--- 574,580 ----
           reftraj_ensemble, langevin_ensemble)
         ! Do Nothing
      CASE(nve_ensemble, nvt_ensemble, nvt_adiabatic_ensemble, npt_i_ensemble, &
!         npt_f_ensemble, npe_i_ensemble, npe_f_ensemble, npt_y_ensemble)
         IF (ensemble==nve_ensemble) check = do_shell
         IF (check) THEN
            SELECT CASE(region)
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_colvar_restraint_constraint_print/velocity_verlet_control.F src/velocity_verlet_control.F
*** src_colvar_restraint_constraint_print/velocity_verlet_control.F	2012-06-04 15:50:56.843703288 +0100
--- src/velocity_verlet_control.F	2012-06-04 16:30:13.519095964 +0100
***************
*** 17,26 ****
         isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
         npt_i_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, &
!        reftraj_ensemble
    USE integrator,                      ONLY: &
         isokin, langevin, nph_uniaxial, nph_uniaxial_damped, npt_f, npt_i, &
!        nve, nve_respa, nvt, nvt_adiabatic, reftraj
    USE md_environment_types,            ONLY: get_md_env,&
                                               md_environment_type
    USE simpar_types,                    ONLY: simpar_type
--- 17,26 ----
         isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, &
         npt_i_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, &
!        reftraj_ensemble, npt_y_ensemble
    USE integrator,                      ONLY: &
         isokin, langevin, nph_uniaxial, nph_uniaxial_damped, npt_f, npt_i, &
!        nve, nve_respa, nvt, nvt_adiabatic, reftraj, npt_y
    USE md_environment_types,            ONLY: get_md_env,&
                                               md_environment_type
    USE simpar_types,                    ONLY: simpar_type
***************
*** 85,90 ****
--- 85,92 ----
         CALL isokin (md_env, error)
      CASE (npt_i_ensemble)
         CALL npt_i (md_env, globenv, error)
+     CASE (npt_y_ensemble)
+        CALL npt_y (md_env, globenv, error)
      CASE (npt_f_ensemble)
         CALL npt_f (md_env, globenv, error)
      CASE (nph_uniaxial_ensemble)
