Fully periodic QM/MM simulations

[Jadzia]

Dear All, 

I'm relatively (but not completely) new to CP2K, and am currently trying to run a QM/MM simulation (at the moment just a test system with water) which is fully periodic, using DFTB3 with dispersion correction, and electrostatic embedding (Coulomb).

As far as I could find out, with DFTB one can set the QM box size to the system size, and in this case one does not need to decouple the QM images in the QMMM section. Does this mean that I can leave away the &PERIODIC section completely in the QM/MM region? Because if I have &PERIODIC sections in the MM and in the DFT sections, then the periodicity should alraedy be set up correctly, right? I have already tried to switch off only the &MULTIPOLE feature, but then I get an error about to a missing keyword related to EWALD:

Value requested, but no value set getting value from keyword GMAX of section EWALD. 

Therefore I guess I need to renmove the &PERIODIC keyword in the QMMM section, and if I remove it the simulation seems to run fine. But I can't tell if the periodicity works as it should. I also read a lot of topics in the mailing list, but I couldn't figure things completely out. My input files are attached, and about any comments I would be grateful. 

Also, there is one other problem. After a while of simulation time, the QM water molecules are slowly falling apart. I.e. a hydrogen would detach from one water molecule and goes to the oxygen of another one. This happens almost always when MM water molecules interact with QM water molecules, but very few times I saw also a QM water molecule falling apart which was not directly near an MM water. I suspect that the electron spillout problem (as it is called in your 2005 article) is the cause. In the article is written that a Gaussian charge distribution is used to solve the electron spillout problem. Is this mechanism always active in QM/MM force evaluations in CP2K, or do I need to active it? (I could not find much about that topic.) Or do you have other recommendations on how to solve the water instability problem?

If you see anything other unusual in my input file I would be glad to as well. I have checked things as best as possible, but don't have much experience with CP2K. 

Best wishes,

Jadzia

[Matt W]

Hi,

you want to keep the periodic section, if you want a periodic system. This _should_* make the QM region feel a potential that comes from the periodically replicated MM region (so essentially Ewald sum). Without that section you just get the potential from the point charges surrounding the region, that needs lots of care about boundary effects otherwise you can get stray fields and other artifacts.

Multipole correction is (normally a fairly small correction) that allows you to use a smaller QM box but still be periodic. You can ignore it if your QM and QMMM cells are the same size.

* I don't know if the periodic functionality works with DFTB ...

Matt

[Jadzia]

Hi Matt, 

Thank you very much for your answer. 

Since previously I did not have coupled the (DFTB3) QM region (which is the same size of the MM region) to the periodic MM images, if I understand you correctly this could be the cause that some of my water molecules fall apart since I did not take many measures against boundary effects as you have written. The only thing I have done is to add standard LJ pontentials for the MM/QM atoms close to each other. What other measures are there which one could take?

I have now tried to set up fully periodic QM/MM:

  &QMMM

    &CELL

      ABC 41.576 41.369 42.526

      PERIODIC XYZ

    &END CELL

    ECOUPL Coulomb

    &PERIODIC

      GMAX 5.0E-01

      &MULTIPOLE off

      &END MULTIPOLE

      &POISSON

        &EWALD

          EWALD_TYPE spme

          GMAX 1

        &END EWALD

      &END POISSON

    &END PERIODIC

    CENTER setup_only

    CENTER_GRID FALSE

    USE_GEEP_LIB 7

    NOCOMPATIBILITY

    &FORCEFIELD

      &NONBONDED14

        @include ../../../cp2k.in.qmmm.lj.system1

      &END NONBONDED14

    &END FORCEFIELD

  &END QMMM

However, I got:

 *******************************************************************************

 *   ___                                                                       *

 *  /   \                                                                      *

 * [ABORT]                                                                     *

 *  \___/     QM-QM long range correction not possible with COULOMB coupling   *

 *    |                                                                        *

 *  O/|                                                                        *

 * /| |                                                                        *

 * / \                                                    qmmm_per_elpot.F:379 *

 *******************************************************************************

Thus it seems I can either

1) switch to mechanical embedding (i.e. ECOUPL NONE)

2) switch off the QM/MM electrostatic periodicity as I had before, in which case I would need to add more of the special measures you mentioned it seems.

And recommendations?

I used DFTB3 without hydrogen corrections. Could that be the problem, or maybe the dispersion correction enhances the electron spillout problem?

Best wishes,

Jadzia

[Jadzia]

Another thing I wonder about: Matt you said that I should keep the section QMMM periodic section so that the QM region feels the potential of the periodically replicated MM region. Which I would call a "MM-QM long range correction". But the error message of my last post states that a QM-QM long range correction is not possible with Coulomb coupling when I add this section. Therefore I wonder if the periodic section which I set up for the MM-QM long range correction was indeed what you had in mind...

[Matt W]

As I said, I wasn't sure if there was a periodic correction implemented using Coulomb coupling for DFTB/semi-empirical.

Looks like you just have to run without periodic coupling, or switch to DFT to use periodic coupling.

[Jadzia]

Thank you for your answer Matt. 

Since error messages can be very helpful if they are precise, I would like summarize my understanding so far related to the error message I have gotten:

QM-QM long range correction not possible with COULOMB coupling.

Firstly, from what you wrote I implicitly conclude that the periodic section which I have added to my input file went into the direction of what you had in mind when you wrote about including the coupling of the QM region with the periodic MM images. In this case the above error message which I got seems to be not fully suitable since I had disabled the QM-QM long range correction (i.e. the MULTIPOLE section). To my understanding (which is admittedly incomplete) the problem might be of a more general nature and summarized as follows: QM/MM long range correction not possible with COULOMB coupling

because it seems to cover both cases QM-MM and optionally also QM-QM long term corrections (if the multipole section is not switched off) for the QM region. 

Secondly, the error message seems to be too general in another aspect, because you wrote that QM/MM periodic couplings are possible with DFT (and Coulomb coupling I guess), in which case the error message does not hold in its current general form. 

[Matt W]

Hi again,

sorry I don't really have time to give long enough replies to be clear, so perhaps a bit confusing.

QM-QM coupling is not needed in your case - your QM box is the same size as the QMMM box. It calculates multipole moments of the QM charge distribution and subtracts artificial interactions if the QM box is smaller than QMMM and adds back in correct ones.

QM-MM long range interactions as you say are the about including the ewald summed potential from the classical part - I think that this is not possible except using E_COUPL = GAUSS, which in turn requires DFT. See the paper from Laio et al for details.

There are several different QMMM methods available, and they have different ranges of applicability.

Matt

[Jadzia]

Hi everyone, 

@Matt: Thank you very much for the clarification, it was helpful.

Since full QM/MM periodicity is not possible, I now would like to get the QM/MM simulation of water working without QM/MM periodicity, but with electrostatic embedding (E_COUPL = COULOMB). The QM method is still DFTB3. 

And I have tried almost everything to get it working, but my water molecules still fall apart. 

I have tried...

• different GMAX values in the EWALD section
• different ALPHA values in the EWALD section
• different O_SPLINE values in the EWALD section
• the QM part completely without electrostatics (no POISSON section)
• DFTB3 with and without dispersion correction
• different time step sizes
• Standard BOMD and second generation CPMD
• ...

Nothing seems to really solve the problem. I have also tried more things such as mechanical embedding with and without QM/MM periodicity. 

Here is a link to a sample input package which can be started directly, all files included:

https://drive.google.com/open?id=0B-Ij17NhUNZoeUpETFJhS1hpbFE

Any help I would greatly appreciate. 

Many thanks,

Jadzia

[Jadzia]

A short update: If I run the water simulation with the same QM settings but QM only (instead of QM/MM), then it seems to run stable and the water is not falling apart. Therefore DFTB3 and water seems to work generally with my setup, and the problem seems to arise only when QM/MM is switched on. Any suggestions or comments are greatly appreciated on how to stabilize the water molecules when using QM/MM. 

Many thanks,

Jadzia

[Yingchun Zhang]

Hi, Jadzia

  just some suggestion:

  

  1, If you want to avoid the spill out problem as you indicated in you original post, you might want to use E_COUPL = GAUSS;

  

  2, check if your pure MM run of the whole systems works.

Yingchun

[Jadzia]

Hi Yingchun, 

thank you so much for your answer, I'm really glad about any comments.

Before I switched to QM/MM I have run a lot MM simulations of the same system, and they run all fine.

The first E_COUPL which I wanted to try was GAUSS as you had suggest, however CP2K rejected my request telling me E_COUPL=GAUSS is not implemented for DFTB.

As you seemed to indicate, I also guess that the electron spillout problem might be the cause. Closely related I could also imagine that the problem might be caused by some other electrostatic related issue. 

What would you try next to find the cause or a solution?

Many thanks,

Jadzia

[Jadzia]

I think I found the CP2K setting which is causing the problem: When I remove the DFTB3_PARAM settings, then my water molecules are not falling apart anymore. 

I will open a new topic because this issue not anymore about periodic QM/MM.