Energy conservation in periodic QM/MM

[Jonggu Jeon]

Dear CP2K users,

I'd like to ask some questions on the capabilities of the QM/MM part of 
CP2K (I'm using version 2.5.1, svn:13632).

I've been attempting a QM/MM MD simulation with CP2K. My system consists of 
one NMA molecule (TZV2P/BLYP DFT) and 48 methanol molecules (OPLS MM force 
field).
The relevant parts of my input file are attached at the end of the message. 
Basically, I'm following all the advices from this group for a fully 
periodic QM/MM setup and set the cell sizes for the QM (in QMMM/CELL) and 
MM (in SUBSYS/CELL) parts equal.

The following is the summary of my findings.

1. With FORCE_EVAL/QMMM/CENTER EVERY_STEP, the total energy is not 
conserved, as pointed out by Dr. Laino in this group before.
2. With FORCE_EVAL/QMMM/CENTER SETUP_ONLY, the total energy is very well 
conserved within 0.1 kcal/mol for a few ps until the QM molecule eventually 
crosses the cell boundary.  At this point system energy shoots up by 
several kcal/mol and the output begins to produce the following messages:

 *** 16:11:03 WARNING in qmmm_util:apply_qmmm_walls_reflective :: One or  
***
 *** few QM atoms are within the SKIN of the quantum box. Check your run  
***
 *** and you may possibly consider: the activation of the QMMM WALLS      
***
 *** around the QM box, switching ON the centering of the QM box or       
***
 *** increase the size of the QM cell. CP2K CONTINUE but results could be 
***
 *** meaningless. qmmm_util.F line 206                                    
***

My questions are as follows:

1. Is it possible to achieve the energy conservation in a fully periodic 
DFT/MM MD without worrying about the QM molecule crossing the cell 
boundary? I believe this is possible in full DFT MD but my experience on 
CP2K QM/MM so far indicates otherwise. I'd appreciate if someone can 
provide an answer and keywords to achieve it.

2. The QMMM/CENTER option affects energy conservation a lot. When CENTER 
EVERY_STEP is used, would it improve energy conservation if I use finer 
grids by increasing DFT/MGRID/CUTOFF value?

The skeletal form of my input file is attached below.
Thank you.

Jonggu Jeon



 &MOTION
   &MD
     ENSEMBLE  NVE
     STEPS  40000
     TIMESTEP     0.5
     STEP_START_VAL  0
     TIME_START_VAL    0
   &END MD
 &END MOTION
 &FORCE_EVAL
   METHOD  QMMM
   &DFT
     BASIS_SET_FILE_NAME /opt/cp2k/2.5.1/tests/QS/GTH_BASIS_SETS
     POTENTIAL_FILE_NAME /opt/cp2k/2.5.1/tests/QS/POTENTIAL
     WFN_RESTART_FILE_NAME wfn.rst
     CHARGE  0
     &SCF
       MAX_SCF  30
       EPS_SCF     9.9999999999999995E-07
       SCF_GUESS  RESTART
       &OT  T
         MINIMIZER  DIIS
         PRECONDITIONER  FULL_ALL
         ENERGY_GAP     1.0000000000000000E-03
       &END OT
       &OUTER_SCF  T
         EPS_SCF     9.9999999999999995E-07
         MAX_SCF  1
       &END OUTER_SCF
     &END SCF
     &QS
       EPS_DEFAULT     1.0000000000000000E-10
       MAP_CONSISTENT  T
       EXTRAPOLATION  ASPC
       EXTRAPOLATION_ORDER  3
       METHOD  GPW
     &END QS
     &MGRID
       CUTOFF     2.8000000000000000E+02
       COMMENSURATE  T
     &END MGRID
     &XC
       &XC_GRID
         XC_SMOOTH_RHO  NN50
         XC_DERIV  SPLINE2_SMOOTH
       &END XC_GRID
       &XC_FUNCTIONAL  NO_SHORTCUT
         &BECKE88  T
         &END BECKE88
         &LYP  T
         &END LYP
       &END XC_FUNCTIONAL
       &VDW_POTENTIAL
         POTENTIAL_TYPE  PAIR_POTENTIAL
         &PAIR_POTENTIAL
           TYPE  DFTD3
           PARAMETER_FILE_NAME /opt/cp2k/2.5.1/tests/QS/dftd3.dat
           REFERENCE_FUNCTIONAL BLYP
           CALCULATE_C9_TERM  F
           REFERENCE_C9_TERM  F
           LONG_RANGE_CORRECTION  F
         &END PAIR_POTENTIAL
       &END VDW_POTENTIAL
     &END XC
   &END DFT
   &MM
     &FORCEFIELD
       PARMTYPE  AMBER
       PARM_FILE_NAME ndmd48.ao.prmtop
       &SPLINE
         EMAX_SPLINE     1.0000000000000000E+00
       &END SPLINE
     &END FORCEFIELD
     &POISSON
       &EWALD
         EWALD_TYPE  SPME
         ALPHA     3.4999999999999998E-01
         GMAX  18
       &END EWALD
     &END POISSON
   &END MM
   &QMMM
     E_COUPL  GAUSS
     USE_GEEP_LIB  7
     NOCOMPATIBILITY  T
     CENTER  SETUP_ONLY
     INITIAL_TRANSLATION_VECTOR    -1.0888775605157457E+01   
-2.9066675351256084E+00    7.3922596274008097E+00
     &CELL
       ABC     1.5084868700000001E+01    1.5084868700000001E+01    
1.5084868700000001E+01
       PERIODIC  XYZ
     &END CELL
     &PERIODIC
       GMAX     5.0000000000000000E-01
       &MULTIPOLE
         RCUT     1.2000000000000000E+01
         EWALD_PRECISION     4.9999999999999998E-07
         ANALYTICAL_GTERM  T
         NGRIDS  50 50 50
       &END MULTIPOLE
     &END PERIODIC
   &END QMMM
   &SUBSYS
     &CELL
       A     1.5084868700000001E+01    0.0000000000000000E+00    
0.0000000000000000E+00
       B     0.0000000000000000E+00    1.5084868700000001E+01    
0.0000000000000000E+00
       C     0.0000000000000000E+00    0.0000000000000000E+00    
1.5084868700000001E+01
       MULTIPLE_UNIT_CELL  1 1 1
     &END CELL
     &TOPOLOGY
       NUMBER_OF_ATOMS  156
       CONN_FILE_NAME ndmd48.ao.prmtop
       CONN_FILE_FORMAT  AMBER
       MULTIPLE_UNIT_CELL  1 1 1
     &END TOPOLOGY
   &END SUBSYS
 &END FORCE_EVAL

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20150318/b85f4460/attachment.htm>