[CP2K:8626] Questions about AdBF-QM/MM in cp2k

[Noam Bernstein]

On Thu, Feb 2, 2017 at 12:59 AM, Ivan Gladich <igla... at gmail.com> wrote:
> Dear Cp2K users,
>
> I am a newbie of CP2K and I am trying to use Adaptive Buffered force QM-MM
> (AdBF-QM/MM) in cp2k
>
> https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4351341/
>
>  to study a radical specie in water
> I have a couple of questions.
>
> 1) Using
>
> ...
>
>   &QMMM
>
>
>     &QM_KIND Br
>
>       MM_INDEX  1
>
>     &END QM_KIND
>
>
>
>     &QM_KIND O1
>
>       MM_INDEX  2
>
>     &END QM_KIND
>
>
>     &FORCE_MIXING
>
> MAX_N_QM 500
>
>          R_CORE 0.0 0.1 #dimension of the core region for the reduced QM/MM
> calculation. Could be 0 if you have the force field
>
>          R_QM   3.0 3.5  #Dimension of the dynamical core region
>
>          R_BUF  3.0 3.5  #dimensione of the buffere region
>
>          QM_KIND_ELEMENT_MAPPING Br Br #atoms that can be adaptively
> selectefor the QM region
>
>          QM_KIND_ELEMENT_MAPPING O O  #atoms that can be adaptively selected
> for the QM region
>
>          QM_KIND_ELEMENT_MAPPING H H  #atoms that can be adaptively
> selectedfor the QM region
>
>          EXTENDED_DELTA_CHARGE  0 # this is the charge of the QM_dynamical +
> buffer region.
>
>          MOMENTUM_CONSERVATION_REGION  QM #only to dynamical core, you can
> also choose BUFFER for QM+BUFFER
>
>          MOMENTUM_CONSERVATION_TYPE   EQUAL_A #equal acceleration on each
> atom, or EQUAL_F, equal force
>
>    &END FORCE_MIXING
>
> ....
>
> And I got in my output
>
> QMMM FORCE MIXING final count (not including links):  N_QM core_list      2
> N_QM core     0 N_QM extended    12 N_QM buffered    63
>
> QMMM FORCE MIXING final count (not including links):  N_QM core_list      2
> N_QM core     0 N_QM extended    12 N_QM buffered    63
>
> QMMM FORCE MIXING final count (not including links):  N_QM core_list      2
> N_QM core     0 N_QM extended    12 N_QM buffered    63
>
> QMMM FORCE MIXING final count (not including links):  N_QM core_list      2
> N_QM core     0 N_QM extended    12 N_QM buffered    63
>
> QMMM FORCE MIXING final count (not including links):  N_QM core_list      2
> N_QM core     0 N_QM extended    12 N_QM buffered    63
>
>
> ...
>
>
> Is N_QM core = 0 the number of atoms in my dynamical QM region?

No, I apologize for the somewhat confusing output.  "N_QM core" is the
number added to the core list by the action of R_CORE.  I.e. The core
region consists of a fixed list (the atoms specified by the
conventional QM/MM directives) + an additional region, which is every
atom within R_CORE of the fixed list.  The core region (combining
those two, the list and the R_CORE based additions) is QM in both
reduced and extended calculations whose forces will be mixed.   The QM
dynamics region includes all those core atoms, plus whatever is
included by R_QM (based on distance from core atoms).

>
> Moreover, reading the paper N_QM extended "is constructed adding a buffer
> region around the QM dynamical core"
>
> But N_QM extended =12 and N_QM buffered is 63, so it cannot be (N_QM should
> be the same of the buffer plus the dynamical core)...

I think here the notation of the code and the paper are unfortunately
not quite consistent.  N_QM extended refers to what the paper calls
the dynamical QM region (I'm pretty sure it's actually just the part
that's in addition to the core).

It doesn't help that there are dynamically selected atoms in the core
region (as well as in the QM dynamics region and buffer region), but
that shouldn't be confused with the region that follows QM force
dynamics.

>
>
> 2) How does the R_BUF and R_QM works?
>
>
>  R_BUF is a radius around each atoms belonging to the QM dynamical core
> region
>
>
> OR
>
>
> R_BUF is centered around the QM_CORE (in my case determined by Br and O
> position) and, thus, I should impose always R_QM =< R_BUF ?

More like the former, but not exactly, at least not if you're using
the terminology of either the paper or the code.  As the code
documentation says, it's applied around every atom that should follow
QM forces, i.e. the core list + the core distance (R_CORE) based + the
QM dynamical region distance (R_QM) based.

>
>
> 3) Is it possible to write, at each time step of my run, a label or the
> indexes of the atoms in the dynamical and extended QM/MM region?
>
>
> I tried using
>
>
>       &FORCE_MIXING_LABELS
>
>          FORMAT XYZ
>
>          &EACH
>
>             MD 1
>
>          &END EACH
>
>          COMMON_ITERATION_LEVELS 0
>
>       &END FORCE_MIXING_LABELS
>
>
> but it seems that it labels only the atoms in the QM CORE, i.e.,
> *fmlabels-1_23.xyz
>
>
>     6482
>
>  i =       23, time =       11.500, E =       -21.9207476721
>
>  BR        10.0000000000       10.0000000000       10.0000000000
>
>   O        10.0000000000       10.0000000000       10.0000000000
>
>   O         0.0000000000        0.0000000000        0.0000000000
>
>   H         0.0000000000        0.0000000000        0.0000000000
>
>   H         0.0000000000        0.0000000000        0.0000000000
>

Good question.  I'll take a look at the source and remind myself what
kind of functionality is available.

     Noam