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

Rolf David rolf.d... at gmail.com
Sun Mar 22 07:23:29 UTC 2015


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.
>>>>>>>
>>>>>>
>>>>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20150322/7c6d3552/attachment.htm>


More information about the CP2K-user mailing list