How to do metadynamics with CV1 and CV2 and set up restraint using CV3 in cp2k?

[hongxi...@gmail.com]

Dear cp2k users and developers,

I am trying to do a metadynamics in cp2k. I set 3 collective variables. I want to use CV1 and CV2 as collective variables for metadynamics, while at the same time using CV3 for a restraint to constraint the motion of two molecules not too far away from each other. Could anyone help me set this up appropriately? 

My input is as folllows.

I set up the three COLVAR in SUBSYS section first. 

&COLVAR 1

      &COORDINATION

        ATOMS_FROM 2 3

        ATOMS_TO 1 4

        R_0 [angstrom] 1.2

        NN 8

        ND 14

      &END COORDINATION

    &END COLVAR

    &COLVAR 2

      &COORDINATION

        ATOMS_FROM 5 6

        ATOMS_TO 1 4

        R_0 [angstrom] 1.2

        NN 8

        ND 14

      &END COORDINATION

    &END COLVAR

    &COLVAR 3

      &DISTANCE

        POINTS 1 2

        &POINT 1

          ATOMS 2 3

          TYPE GEO_CENTER

        &END POINT

        &POINT 2

          ATOMS 5 6

          TYPE GEO_CENTER

        &END POINT

      &END DISTANCE

    &END COLVAR

And then I add the constraint section in the MOTION section.

&MOTION

 &CONSTRAINT

  &COLLECTIVE

    COLVAR 3

    INTERMOLECULAR .TRUE.

    TARGET [angstrom] 5.0

    &RESTRAINT

      K 100000

    &END RESTRAINT

  &END COLLECTIVE

 &END CONSTRAINT

 &MD

    ENSEMBLE NVT

    STEPS  400000

    TIMESTEP 0.5

    TEMPERATURE 298.0

    &THERMOSTAT

      TYPE CSVR

      REGION MASSIVE

      &CSVR

        TIMECON [fs] 200.0

      &END

    &END

  &END MD

  &FREE_ENERGY

    &METADYN

      DO_HILLS

      WELL_TEMPERED

      DELTA_T 50000

      NT_HILLS 30

      WW 0.005

      &METAVAR

        SCALE 0.2

        COLVAR 1

      &END METAVAR

      &METAVAR

        SCALE 0.2

        COLVAR 2

      &END METAVAR

      &PRINT

        &COLVAR

           COMMON_ITERATION_LEVELS 2

           &EACH

             MD 1

           &END

        &END

        &HILLS

           COMMON_ITERATION_LEVELS 2

           &EACH

             MD 1

           &END

        &END

      &END

    &END METADYN

  &END FREE_ENERGY

Does anyone know what I did wrong and how to set it up appropriately instead? With this input, I got super huge potential energy and then the geometry fly away.

Thanks in advance!

Sincerely

Hongxia

[Matt D]

Hi Hongxia,

I guess because you have set DeltaT to 50000K... I'm not surprised the atoms fly apart in this case. 

If this is your first metadynamics run, I'd recommend running an unbiased simulation first (i.e. not WELL_TEMPERED and no DELTA_T set) and see how the CV evolves. You can then set the width of the Gaussians (SCALE) according to variations in the value of the CV with time, for example by taking the standard deviation of the CV and setting this to be SCALE. Once you have set an appropriately sized gaussian, you can then introduce the bias. The bias factor is given by (T+deltaT)/T. Given 298K is your starting temperature and you were using 50000 as deltaT, your bias factor is rather large. A good initial range might be somewhere in the region of 5 or 10 for the bias factor, so deltaT around 1000 - 2000 K. Of course this is system dependent so you need to try several values, the higher the value the more of free-energy space you will explore, though you want to make sure it's comparable to the free energy barriers you want to traverse to avoid getting situations where you explode the simulation as you have done by using 50000K for deltaT.

Best wishes,

Matt