Hi,<div><br></div><div>I am trying a NVT MD simulation for a QM/MM system with QM water molecules surrounding a spatially fixed MM molecule, which are in turn surrounded by MM water molecules. I am using a time step of 5 a.u. (for particular reasons) and the QM box is much smaller than compared to full MM box so within the timescale of MD simulation the QM molecules do not cross the QM box boundary. I am attaching the relevant sections of input file:</div><div>----</div><div><div> &MOTION</div><div>   &MD</div><div>     ENSEMBLE  NVT</div><div>     STEPS  10000</div><div>     TIMESTEP     1.21E-01</div><div>     TEMPERATURE     3.00E+02<br></div><div>     &THERMOSTAT</div><div>       TYPE  NOSE</div><div>       REGION  DEFINED</div><div>       &DEFINE_REGION</div><div>         QM_SUBSYS  ATOMIC</div><div>       &END DEFINE_REGION</div><div>       &DEFINE_REGION</div><div>         MM_SUBSYS  ATOMIC</div><div>       &END DEFINE_REGION</div><div>       &NOSE</div><div>         TIMECON     9.99E+02</div><div>       &END NOSE</div><div>     &END THERMOSTAT</div></div><div><br></div><div><div> &FORCE_EVAL</div><div>   METHOD  QMMM</div><div>   &DFT</div><div>     BASIS_SET_FILE_NAME ./GTH_BASIS_SETS</div><div>     POTENTIAL_FILE_NAME ./GTH_POTENTIALS</div><div>     CHARGE  0</div><div>     &SCF</div><div>       MAX_SCF  60</div><div>       EPS_SCF     2.0E-07</div><div>       SCF_GUESS  RESTART</div><div>       &OT  T</div><div>         MINIMIZER  DIIS</div><div>         PRECONDITIONER  FULL_SINGLE_INVERSE</div><div>       &END OT</div><div>       &OUTER_SCF  T</div><div>         EPS_SCF     2.0E-07</div><div>         MAX_SCF  40</div><div>       &END OUTER_SCF</div><div>       &PRINT</div><div>         &RESTART  SILENT</div><div>           BACKUP_COPIES  1</div><div>         &END RESTART</div><div>       &END PRINT</div><div>     &END SCF</div><div>     &QS</div><div>       EXTRAPOLATION  ASPC</div><div>       EXTRAPOLATION_ORDER  4</div><div>     &END QS</div><div>     &MGRID</div><div>       CUTOFF     3.2E+02</div><div>       COMMENSURATE  T</div><div>     &END MGRID</div><div>     &XC</div><div>       DENSITY_CUTOFF     1.0E-10</div><div>       GRADIENT_CUTOFF     1.0E-10</div><div>       TAU_CUTOFF     1.0E-10</div><div>       &XC_GRID</div><div>         XC_SMOOTH_RHO  NN50</div><div>         XC_DERIV  SPLINE2_SMOOTH</div><div>       &END XC_GRID</div><div>       &XC_FUNCTIONAL  NO_SHORTCUT</div><div>         &PBE  T</div><div>           PARAMETRIZATION  REVPBE</div><div>         &END PBE</div><div>       &END XC_FUNCTIONAL</div><div>       &VDW_POTENTIAL</div><div>         POTENTIAL_TYPE  PAIR_POTENTIAL</div><div>         &PAIR_POTENTIAL</div><div>           R_CUTOFF     1.9E+01</div><div>           TYPE  DFTD3</div><div>           PARAMETER_FILE_NAME ./dftd3.dat</div><div>           REFERENCE_FUNCTIONAL revPBE</div><div>           CALCULATE_C9_TERM  T</div><div>           REFERENCE_C9_TERM  T</div><div>           LONG_RANGE_CORRECTION  F</div><div>         &END PAIR_POTENTIAL</div><div>       &END VDW_POTENTIAL</div><div>     &END XC</div><div>   &END DFT</div><div>   &MM</div><div>     &FORCEFIELD</div><div>       ...</div><div>     &END FORCEFIELD</div><div>     &POISSON</div><div>       &EWALD</div><div>         EWALD_TYPE  SPME</div><div>         RCUT     1.2E+01</div><div>         ALPHA     2.917E-01</div><div>         GMAX  75</div><div>         O_SPLINE  4</div><div>       &END EWALD</div><div>     &END POISSON</div><div>   &END MM</div><div>   &QMMM</div><div>     E_COUPL  GAUSS</div><div>     USE_GEEP_LIB  6</div><div>     NOCOMPATIBILITY  T</div><div>     CENTER  NEVER</div><div>     INITIAL_TRANSLATION_VECTOR     0.0E+00    0.0E+00    0.0E+00</div><div>     &QM_KIND O</div><div>       MM_INDEX  ...</div><div>     &END QM_KIND</div><div>     &QM_KIND H</div><div>       MM_INDEX  ...</div><div>     &END QM_KIND</div><div>     &MM_KIND H</div><div>       RADIUS     4.4E-01</div><div>     &END MM_KIND</div><div>     &MM_KIND O</div><div>       RADIUS     1.2E+00</div><div>     &END MM_KIND</div><div>     .....</div><div>     &CELL</div><div>       ABC     4.0E+01    3.0E+01    3.0E+01</div><div>       PERIODIC  XYZ</div><div>     &END CELL</div><div>     &PERIODIC</div><div>       GMAX     5.0E-01</div><div>       &MULTIPOLE  ON</div><div>         RCUT     1.2E+01</div><div>         EWALD_PRECISION     1.0E-08</div><div>         ANALYTICAL_GTERM  T</div><div>       &END MULTIPOLE</div><div>     &END PERIODIC</div><div>   &END QMMM</div><div>   &SUBSYS</div><div>     &CELL</div><div>       A     5.0E+01    0.0E+00    0.0E+00</div><div>       B     0.0E+00    4.0E+01    0.0E+00      </div><div>       C     0.0E+00    0.0E+00    4.0E+01</div><div>       MULTIPLE_UNIT_CELL  1 1 1</div><div>     &END CELL</div><div>     &COORD</div></div><div><br></div><div>-------</div><div><br></div><div>With this setup for the production run I obtained the following trends for the potential and constant quantity as the graphs attached here. I have tested the QM setup for a fully QM calculation on bulk water, which works fine. My main aim here is to reduce the energy drift for the QM/MM system, so is there any other aspect of the calculation should I take care to reduce the drift? Any suggested corrections for the QM/MM file would also be helpful.</div><div><br></div><div>Best Regards,</div><div>Mayank Dodia</div><div><br></div><div><br></div><div><br></div>