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

[Dorothea Golze]

Hi,

if you have a fully periodic setup where the QM and MM boxes are equal, you do not need the decoupling/recoupling stuff, i.e. switch off explicitly the MULTIPOLE section like that
&MULTIPOLE OFF

&END

Also set in the QMMM section

 &WALLS
     TYPE NONE
  &END WALLS

That should address your first question...

Cheers,

Dorothea

 

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To post to this group, send email to cp...@googlegroups.com.
Visit this group at http://groups.google.com/group/cp2k.
For more options, visit https://groups.google.com/d/optout.

[Jonggu Jeon]

Dear Dorothea,

Thank you for the invaluable suggestion. 
While your answer makes perfect sense to me, I just found out that CP2K 2.5.1 does not have option to turn off QMMM/PERIODIC/MULTIPOLE using the section parameter OFF.
Therefore, CP2K 2.5.1 produces the same (diverging) energy profile whether I set &MULTIPOLE OFF or not.

I will install version 2.6 and try it again.

Best regards,

Jonggu


2015년 3월 18일 수요일 오후 6시 22분 40초 UTC+9, Dorothea Golze 님의 말:

[Matt W]

Hi,

there is also the CENTER_GRID that has been introduced quite recently. I don't know how much it has been tested, maybe someone else can comment? But for my simple test water system it seems to improve the energy conservation considerably. The idea is, I think, to only center the QM molecule in a "quantized" manner on equivalent grid points of the (fine DFT?) grid.

To get really good energy conservation you might need to increase your cutoff/rel_cutoff too.

Matt

[Jonggu Jeon]

Dear all,

Here is the follow-up to my previous message.

I have repeated the QM/MM MD using MULTIPOLE OFF using the latest CP2K trunk binary (svn:15171). The total energy from this new run is *identical* to the previous one using MULTIPOLE section. The only difference is that the decoupling/recoupling energies and the warning messages about QM atoms close to the boundary are not printed in the output file any more.

I have also tried the "CENTER_GRID T" keyword together with CENTER EVERY_STEP as Matt suggested above and its short time energy profile up to 100 fs is nearly the same as the one without CENTER_GRID.  I guess the total energy is still adversely affected by the centering even with the CENTER_GRID.

From these observations, I'm afraid it's likely that the current CP2K cannot handle gracefully the QM atoms crossing boundary.
I'd appreciate any comment or suggestion.

Jonggu


2015년 3월 18일 수요일 오후 9시 13분 29초 UTC+9, Matt W 님의 말:

[Matt W]

Try increasing the cutoff to 600-800 Ry? BLYP functional is known to be noisy.

I have repeated the QM/MM MD using MULTIPOLE OFF using the latest CP2K trunk binary (svn:15171). The total energy from this new run is *identical* to the previous one using MULTIPOLE section. The only difference is that the decoupling/recoupling energies and the warning messages about QM atoms close to the boundary are not printed in the output file any more.

Dorothea pointed out that, _if_ your QM box is the same size as the QMMM box (as you have in the input above), you don't need the de/recoupling - it _should_ give identical results. It just slows the calculation a little.

HTH

Matt

[Jonggu Jeon]

Thank you for the suggestion.

In fact, I'm currently running a new job with MGRID/CUTOFF=480 and CENTER EVERY_STEP. As you can see in the following figure, the early result indicates a very good energy conservation  at the same level as the previous best example (CUTOFF=280 and CENTER SETUP_ONLY).
What troubles me is the fact that computation time increases rapidly with larger CUTOFF and the same level of energy conservation is already possible with smaller cutoff if no centering is used.

However, given the circumstances, I agree that using centering and larger cutoff is the best way to go.

Regards,

Jonggu

2015년 3월 19일 목요일 오후 8시 17분 54초 UTC+9, Matt W 님의 말:

[Rolf David]

I've done test on my systems and reach out some conclusions.

GAPW / B3LYP / 280/40 cutoff / CENTER_EVERY_STEP gave me a slow drift and some fluctuations (for 2ps)

GAPW / B3LYP / 300/50 cutoff / CENTER_EVERY_STEP and with CENTER_GRID TRUE, almost no drift/fluctuations (for 2ps)

I don't know if it's the cutoff of the center_grid true which made an ompatch, but the computational cost is negligible (don't see a difference). I'll see next month what of it when i reach 15/20ps.

[Teodoro Laino]

Hi Rolf,

the point is that the QM part is - strictly speaking - not translational invariant for arbitrary translation vectors. This means that if you re-center at every step for an arbitrary translation vector you will see a drift, given by the fact that an atom - depending on its relative position w.r.t. grid points - may experience numerically different forces (with a deviation, hopefully lower than the norm). Moreover, the situation is even more dramatic because the recenter vector will not be random, but rather subject to a constant direction bias, which will cause a pretty visible drift. 

Your combination: 300 and CENTER_GRID are both going in the direction of: 1) minimising the error in the atomic forces by increasing cutoff 2)Re-center using a translation vectors which is commensurate to grid spacing. Both should help minimise drift.

CENTER_GRID has no computational overhead. 

300 Ryd cutoff instead should marginally increase the cost.

Teo

[Rolf David]

Hi Teo,

That's the conclusion I reached. Thanks for confirming it.

And GPAW is less dependant of the cut-off than GPW because of the hard-density being calculated on local atomic grids. (i'm still trying to "tame the beast" so I might be wrong)

Rolf