Thermostats in QMMM with constraints

[shihao]

Dear all,

I'm trying to run QM/MM simulations with constraint.  The constraint
is placed on the distance difference of three specific atoms (these
three atoms are in two different molecules, which are both in QM
region).  When I use NOSE thermostat and user-defined regions, I got
the error message "*** User Defined thermostats with global
constraints not implemented!  ***". Since MASSIVE thermostat is also
impossible, it seems that the only possible choice is GLOBAL.
However, GLOBAL thermostat cannot effectively control the temperature
in QM region.  Is there any way to go around it?

Thanks a lot.
Shihao

P.S. Here's the input file in case you need to take a look:

@SET CURR_I  01
@SET REPLICA  001
@SET SEED     1000

&GLOBAL
  PROGRAM_NAME                 CP2K
  PROJECT_NAME                 run_${REPLICA}_${CURR_I}
  RUN_TYPE                     MD
  SEED                         ${SEED}
  PREFERRED_FFT_LIBRARY        FFTW
  PRINT_LEVEL                  LOW
  SAVE_MEM
&END GLOBAL

&MOTION
  &MD
    ENSEMBLE                   NVT
    STEPS                      1000
    TIMESTEP                   0.50
    TEMPERATURE                298.15
    &THERMOSTAT
      TYPE                     NOSE
      REGION                   DEFINED
      &NOSE
        LENGTH                 3
        YOSHIDA                3
        TIMECON                100
        MTS                    2
      &END NOSE
      @INCLUDE run_${REPLICA}_cp2k.tr1
    &END THERMOSTAT
    &PRINT
      &ENERGY
        &EACH
          MD                   20
        &END EACH
      &END ENERGY
    &END PRINT
  &END MD

  &PRINT
    &RESTART
      &EACH
        MD                     20
      &END EACH
    &END RESTART
    &RESTART_HISTORY           OFF
    &END RESTART_HISTORY

    &TRAJECTORY                SILENT
      FORMAT                   DCD
      &EACH
        MD                     5
      &END EACH
    &END TRAJECTORY
    &VELOCITIES                OFF
    &END VELOCITIES
    &FORCES                    OFF
    &END FORCES
  &END PRINT

  &CONSTRAINT
    &HBONDS
      EXCLUDE_QM               T
    &END HBONDS
    SHAKE_TOLERANCE            1.0E-9
    &LAGRANGE_MULTIPLIERS ON
      COMMON_ITERATION_LEVELS 3
    &END LAGRANGE_MULTIPLIERS
    &COLLECTIVE
     COLVAR 1
     INTERMOLECULAR T
     TARGET [angstrom]  1.0
    &END COLLECTIVE
  &END CONSTRAINT
&END MOTION

&FORCE_EVAL
  METHOD QMMM

  &DFT
    BASIS_SET_FILE_NAME        ./BASIS_SET
    POTENTIAL_FILE_NAME        ./POTENTIAL
    CHARGE                     1
    MULTIPLICITY               1
    UNRESTRICTED_KOHN_SHAM     F

    &SCF
      SCF_GUESS                ATOMIC
      EPS_SCF                  1.0E-6
      MAX_SCF                  50
      &OUTER_SCF
        MAX_SCF                10
      &END OUTER_SCF
      &OT
        PRECONDITIONER         FULL_SINGLE_INVERSE
        MINIMIZER              DIIS
        N_DIIS                 7
      &END OT
      &PRINT
        &RESTART
          &EACH
            MD                 20
          &END EACH
        &END RESTART
        &RESTART_HISTORY       OFF
        &END RESTART_HISTORY
      &END PRINT
    &END SCF

    &QS
      EPS_DEFAULT              1.0E-12
      EXTRAPOLATION            ASPC
    &END QS
    &MGRID
      COMMENSURATE
      CUTOFF                   280
    &END MGRID
    &POISSON
      POISSON_SOLVER           MT
      PERIODIC                 NONE
    &END POISSON

    &XC
      &XC_FUNCTIONAL           PBE
      &END XC_FUNCTIONAL
      &XC_GRID
        XC_SMOOTH_RHO          NN10
        XC_DERIV               SPLINE2_SMOOTH
      &END XC_GRID
    &END XC

    &PRINT
      &E_DENSITY_CUBE
        &EACH
          MD                   20
        &END EACH
      &END E_DENSITY_CUBE
    &END PRINT
  &END DFT

  &MM
    &FORCEFIELD
      PARMTYPE                 CHM
      PARM_FILE_NAME           ./par_all27_prot_lipid_amm.prm
      &SPLINE
        RCUT_NB                12.0
      &END SPLINE
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE             SPME
        ALPHA                  0.35
        GMAX                   80 80 80
        O_SPLINE               6
      &END EWALD
    &END POISSON
  &END MM

  &QMMM
    USE_GEEP_LIB               7
    E_COUPL                    GAUSS

    @INCLUDE run_${REPLICA}_cp2k.inp

    @INCLUDE mm_kinds

    &WALLS
      TYPE                     REFLECTIVE
      WALL_SKIN                1.5
    &END WALLS

    &PRINT
      &PROGRAM_RUN_INFO        SILENT
      &END PROGRAM_RUN_INFO
      &PERIODIC_INFO           SILENT
      &END PERIODIC_INFO
      &QMMM_LINK_INFO          SILENT
      &END QMMM_LINK_INFO
    &END PRINT
  &END QMMM

  &SUBSYS
    &CELL
      ABC                      92.3237 77.4554 91.24
      PERIODIC                 XYZ
    &END CELL
    &TOPOLOGY
      CONNECTIVITY             UPSF
      CONN_FILE_NAME           ./build.xplor_psf
      COORDINATE               PDB
      COORD_FILE_NAME          ./run_${REPLICA}_cp2k.pdb
      PARA_RES                 T
    &END TOPOLOGY

    ########################################  Basis sets and
pseudopotentials
    &KIND H
      BASIS_SET DZVP-GTH-PBE
      POTENTIAL GTH-PBE-q1
    &END KIND
    &KIND C
      BASIS_SET DZVP-GTH-PBE
      POTENTIAL GTH-PBE-q4
    &END KIND
    &KIND N
      BASIS_SET DZVP-GTH-PBE
      POTENTIAL GTH-PBE-q5
    &END KIND
    &KIND O
      BASIS_SET DZVP-GTH-PBE
      POTENTIAL GTH-PBE-q6
    &END KIND

    &COLVAR
      &COMBINE_COLVAR
        &COLVAR
          &DISTANCE
            ATOMS 2554 27602
          &END
        &END
        &COLVAR
          &DISTANCE
            ATOMS 27600 27602
          &END
        &END
        FUNCTION CV1-CV2
        VARIABLES CV1 CV2
        ERROR_LIMIT 1.0E-9
      &END
    &END
  &END SUBSYS
&END FORCE_EVAL