methane in a water box QM/MM

[Marius Retegan]

Hi,

In my continuing saga of using cp2k for running qm/mm calculations I have encounter yet another bump.
While using the input below I get this error "FORCEFIELD| Missing Bond              (  CT,  HC)" and angle ...

&GLOBAL
  PROJECT mebox_qmmm
  RUN_TYPE MD
  PRINT_LEVEL LOW
&END GLOBAL

&MOTION
  &MD
   ENSEMBLE NVT
   STEPS 1000
   TIMESTEP 1.0
   TEMPERATURE 298
   &THERMOSTAT
     REGION GLOBAL
     TYPE CSVR
     &CSVR
       TIMECON 100
     &END CSVR
   &END THERMOSTAT
  &END MD
  &PRINT
    &RESTART OFF
    &END RESTART
    &RESTART_HISTORY OFF
    &END RESTART_HISTORY
  &END PRINT
&END MOTION

&FORCE_EVAL
  METHOD QMMM
 
  &DFT
    CHARGE 0
    BASIS_SET_FILE_NAME         GTH_BASIS_SETS
    POTENTIAL_FILE_NAME         GTH_POTENTIALS
    &MGRID
      CUTOFF 280
      COMMENSURATE
    &END MGRID
    &SCF
      MAX_SCF     20
      EPS_SCF     1.0E-6
      &OT T
        PRECONDITIONER FULL_ALL
        MINIMIZER DIIS
      &END OT
      &PRINT
        &RESTART OFF
        &END RESTART
        &RESTART_HISTORY OFF
        &END RESTART_HISTORY
      &END PRINT
    &END SCF
    &QS
      EPS_DEFAULT 1.0E-12
      MAP_CONSISTENT
      EXTRAPOLATION ASPC
      EXTRAPOLATION_ORDER 3
    &END QS
    &XC
      &XC_GRID
        XC_SMOOTH_RHO NN10
        XC_DERIV NN10_SMOOTH
      &END XC_GRID
      &XC_FUNCTIONAL BLYP
      &END XC_FUNCTIONAL
    &END XC
    &POISSON
      POISSON_SOLVER MULTIPOLE
      PERIODIC NONE
    &END POISSON
  &END DFT
 
  &MM
    &FORCEFIELD
      VDW_SCALE14 0.5
      EI_SCALE14  0.8333333
      SCALE_CUTOFF .FALSE.
      PARM_FILE_NAME mebox.prmtop
      PARMTYPE AMBER   
      &SPLINE
        EMAX_SPLINE 1.0
        RCUT_NB 10
      &END
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE SPME
        ALPHA .36
        GMAX 23 23 23
      &END EWALD
    &END POISSON
  &END MM
 
  &QMMM
    E_COUPL GAUSS
    USE_GEEP_LIB  7
    &CELL
      ABC 8.0 8.0 8.0
      PERIODIC NONE
    &END CELL
    &QM_KIND C
      MM_INDEX 1
    &END QM_KIND
    &QM_KIND H
      MM_INDEX 2 3 4 5
    &END QM_KIND
  &END QMMM

  &SUBSYS
    &CELL
      ABC [angstrom] 23.2782522  22.7141099  23.2326917 
      PERIODIC XYZ
    &END CELL
    &TOPOLOGY
      COORDINATE PDB
      COORD_FILE_NAME mebox_npt.pdb
      CONNECTIVITY AMBER
      CONN_FILE_NAME mebox.prmtop
    &END
    &KIND C
      BASIS_SET DZVP-GTH
      POTENTIAL GTH-BLYP-q4
    &END KIND
    &KIND H
      BASIS_SET DZVP-GTH
      POTENTIAL GTH-BLYP-q1
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

If I change QMMM with FIST the input runs.

I also have some additional questions, but on somewhat different matters.

In the latest papers that I could find that use cp2k to do qm/mm calculations (http://dx.doi.org/10.1021/ja076081h and http://dx.doi.org/10.1021/jp074858n) the authors used a integration step of 0.5 and 0.48. Unfortunately, they did not discussed their choice. I'm used to running qm/mm md with other packages, but I've always used an integration step of 1 fs. Is cp2k different :-)?

Is there any place where I can find information regarding the parameters xc_deriv and xc_smooth_rho and are "good" choices?

Thanks,
Marius

[Teodoro Laino]

Hi Marius,

you catched another bug connected with the new AMBER parsing feature
when using QMMM.
I'm working right now on a patch.
I will be back to you soon!
Thanks
Teo

[Juerg Hutter]

Hi

0.5 fs is an often used timestep for all atom simulations.
It is adequate for all vibrations including C-H and O-H
stretching. A 1 fs timestep will lead for those high frequency
modes to problems and energy drifts.

regards

Juerg Hutter

----------------------------------------------------------
Juerg Hutter Phone : ++41 44 635 4491
Physical Chemistry Institute FAX : ++41 44 635 6838
University of Zurich E-mail: hut...@pci.uzh.ch
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
----------------------------------------------------------

[Axel]

On Mar 18, 11:12 am, Marius Retegan <marius.s.rete...@gmail.com>
wrote:
> Hi,

>
> If I change QMMM with FIST the input runs.
>
> I also have some additional questions, but on somewhat different matters.
>
> In the latest papers that I could find that use cp2k to do qm/mm

> calculations (http://dx.doi.org/10.1021/ja076081handhttp://dx.doi.org/10.1021/jp074858n) the authors used a integration step of

> 0.5 and 0.48. Unfortunately, they did not discussed their choice. I'm used
> to running qm/mm md with other packages, but I've always used an integration
> step of 1 fs. Is cp2k different :-)?

unless you make _all_ bonds (including the ones in the QM region) that
contain
hydrogens rigid, an MD with a time step larger than 0.5fs will not run
stable
in any code without doing unphysical tricks, or switching to langevin
dynamics
with a significant friction coefficient.
for non-deuterated, unconstrained hydrogens already 0.5fs is pretty
much at
the limit in my experience and if, for example, you have a proton
transferred
across a large barrier, you may get instabilities because of the
proton accelerating
too much. to also give a lower barrier, i found 0.25fs to be a very
safe and
conservative choice (in more than one way ;)).

cheers,
axel.

[Teodoro Laino]

The bug (which was depending on the way of reading natively AMBER prmtop
files) is now fixed in the CVS (thanks for reporting)!
Just few comments about your QMMM setup.
1) The way it is setup you will hardly have a conserved energy (you have
a PERIODIC MM) but you don't use the periodic QMMM version (see PERIODIC
in QMMM). Alternatively you may think about using a spherical cutoff for
the electrostatic (that is tapering the interaction after a certain cutoff).
2) you may want to specify proper covalent radius for the QMMM
electrostatic interactions (at the moment you rely on the default which
is 0.78 (if I'm not wrong) for all atoms).

Have a nice QMMM.
Teo

p.s.: about the xc_deriv and xc_smooth_rho I would not suggest you to
use them unless you really know what you are doing. Most of the time you
don't need to specify those keywords (only for few pathological cases.
CH4 or H2O are definitely not one of them).

[Marius Retegan]

Thanks a lot for all your comments.

Marius

[Marius Retegan]

Hello again,

I was a little hasty with my previous email, since some things are a little unclear.
I looked at the PERIODIC section in QMMM but sincerely that was of very little help.
So if I want to do periodic QMMM, I should remove the decoupling from the QM part and add CELL information in the QMMM part. Additionally in the QMMM part I should also add PERIODIC. If this is correct, I'm stuck since I don't know what parameters to chose next.

Usually, I'm not used to asking so many questions regarding the softwares that I used, since I try as much as possible to find the answer alone. With cp2k (especially the QM/MM part) this is very difficult, mostly due to the its particular nature (plainwaves, original techniques for computing QM/MM electrostatic potential, etc...).
I saw many times being said on the list, that the test are hardly something to take as viable production runs, since they are made for speed. The manual is OK for constructing your input, but it provides very little help in choosing the "right" parameters for your system. I understand that parameters depend on the system studied, but some heads-up on what is safe to use would be nice. What happend with the wiki based documentation proposal? I think that having a great software is important, but giving new users a place to start is equally important.

Please don't take this email as hostility towards cp2k (very far from my intentions), but more as lines coming from a "frustrated" user who tries to see the light at the end of the tunnel :).

Thanks
Marius

[Axel]

On Mar 18, 2:29 pm, Marius Retegan <marius.s.rete...@gmail.com> wrote:

> Usually, I'm not used to asking so many questions regarding the softwares
> that I used, since I try as much as possible to find the answer alone. With

marius,

if you can demonstrate that you made a reasonable effort to
find out things on your own, nobody is going to going to complain,
and most are happy to help. the irritating questions are questions
that say something along the lines of "can somebody give me a step-by-
step
documentation for how to do <insert-a-horribly-complicated-method-
here>
and by the way i have never ever used cp2k before, so please explain
everything to me.

> cp2k (especially the QM/MM part) this is very difficult, mostly due to the
> its particular nature (plainwaves, original techniques for computing QM/MM
> electrostatic potential, etc...).

it _is_ difficult. the only way to make it easy is to hide all the
details from you and have the software use hardcoded defaults. even
though it is tough at the beginning, there is some advantage to having
to fight your way through it, you get a much deeper understanding of
what is going on and can use cp2k much better afterwards. if you look
at very old messages in this forum, you can find some similarly
frustrated
messages written by me as well.

> I saw many times being said on the list, that the test are hardly something
> to take as viable production runs, since they are made for speed. The manual
> is OK for constructing your input, but it provides very little help in
> choosing the "right" parameters for your system. I understand that
> parameters depend on the system studied, but some heads-up on what is safe
> to use would be nice. What happend with the wiki based documentation
> proposal? I think that having a great software is important, but giving new
> users a place to start is equally important.

yes, but who is in the best position to do this?
the person who is currently struggling with it, i.e. you!

there are always a lot of good ideas and sometimes even people
get started, but it is difficult to keep the momentum, particularly
if there is nobody else stepping up to get involved, too.
i can give you an example from my personal experience.
a lot of what was done in the last few years in updating the
manual of the CPMD package and the related tutorials were written
while i was learning how to use and program CPMD myself. now
that i've also contributed a few features, i learned that if you
program something for your own needs, or to help somebody you
are working with, you rarely have the energy left to document
everything so that a beginner can deal with it, because for
yourself it is all very obvious. but if you are trying to learn
it yourself, it can be even helpful to write down what you understood
and have more experienced people point out your mistakes to you.

every problem is also an opportunity.

cheers,
axel.

> Please don't take this email as hostility towards cp2k (very far from my
> intentions), but more as lines coming from a "frustrated" user who tries to
> see the light at the end of the tunnel :).

>
> Thanks
> Marius
>
> On Wed, Mar 18, 2009 at 5:18 PM, Marius Retegan

> <marius.s.rete...@gmail.com>wrote:

>
> > Thanks a lot for all your comments.
>
> > Marius
>

[Teodoro Laino]

Hi Marius,

we discussed so many time about the documentation that is almost an
exhausted thread :-).
Anyway.. Axel started with the pages on the google forum and they are
exactly in the state he left them..
Why don't you take over and contribute to them?

Let's just go back now to your question..
When you want to do a PERIODIC QMMM calculation, conceptually speaking
you have to decouple and recouple the QM images with the proper periodicity.
have a look at:
QMMM/QS/regtest-4/crys_per_qmmm.inp, which is true that it is a regtest
but is still ok as a starting point..

First don't forget in all QMMM runs that use quickstep to add
COMMENSURATE to the MGRID section. It is mandatory!

Second if you do fully periodic QMMM you don't have to specify any
poisson section in the DFT section (since the decoupling is handled
directly with the re-coupling).
In the CELL section of the QMMM section you have to enable fully
PERIODIC (i.e. delete possible PERIODIC NONE in both cell sections) and
add to the QMMM section the PERIODIC section:


&PERIODIC
GMAX 0.5 (is a reasonable value.. but depends on the GEEP expansion ..
USE_GEEP_LIB 10 is good enough)
&MULTIPOLE
EWALD_PRECISION 1.0E-6
RCUT CELL/2 (is good enough)
&END
&END

at the moment the only way to decouple and recouple is with the
MULTIPOLE scheme.
You can start from these hints and check the set of parameters better
for you.
And in case you have still time and willingness put down a couple of
pages on how to use the QMMM in cp2k.
Future users will be extremely glad. ;-)
Cheers
Teo

[Marius Retegan]

Hello,

Thank you for the help. At the moment I don't feel sufficiently comfortable with writing a qm/mm tutorial, since I just started using cp2k. As soon as things start to crystallize, I'm going to write a draft, that other, more experience people, can polish. Until then it's possible that I'm going to spam the list with other questions :-).

Marius