[CP2K:6271] Energy conservation in periodic QM/MM

Rolf David rolf.d... at gmail.com
Mon Mar 23 08:53:56 UTC 2015


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

On Sunday, March 22, 2015 at 9:38:05 AM UTC+1, Teo wrote:
>
> 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
>
> On 22 Mar 2015, at 08:23, Rolf David <rolf... at gmail.com <javascript:>> 
> wrote:
>
> 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.
>
>
>
> On Thursday, March 19, 2015 at 12:48:29 PM UTC+1, Jonggu Jeon wrote:
>>
>> 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
>>
>>
>> <https://lh6.googleusercontent.com/-P0Moyh0OT5k/VQq0qY--uYI/AAAAAAAAA3w/6xnksES21jA/s1600/te.png>
>>
>>
>>
>> 2015년 3월 19일 목요일 오후 8시 17분 54초 UTC+9, 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
>>>
>>>  
>>>
>>>> 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 님의 말:
>>>>>
>>>>> 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
>>>>>
>>>>> On Wednesday, March 18, 2015 at 11:22:10 AM UTC, Jonggu Jeon wrote:
>>>>>>
>>>>>> 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 님의 말:
>>>>>>>
>>>>>>> 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
>>>>>>>  
>>>>>>>
>>>>>>> 2015-03-18 8:46 GMT+01:00 Jonggu Jeon <jeon... at gmail.com>:
>>>>>>>
>>>>>>>> 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
>>>>>>>>
>>>>>>>>
>>>>>>>> -- 
>>>>>>>> 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+... at googlegroups.com.
>>>>>>>> To post to this group, send email to cp... at googlegroups.com.
>>>>>>>> Visit this group at http://groups.google.com/group/cp2k.
>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>
>>>>>>>
>>>>>>>
> -- 
> 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+... at googlegroups.com <javascript:>.
> To post to this group, send email to cp... at googlegroups.com <javascript:>.
> Visit this group at http://groups.google.com/group/cp2k.
> For more options, visit https://groups.google.com/d/optout.
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20150323/821beebf/attachment.htm>


More information about the CP2K-user mailing list