metadynamics energy convergence

[SAURABH SINGH]

Hii 

I am facing a problem in metadynamics study

I am using normal metadynamics.
In my case I am getting the multiple recrossing, but my energy is not converging 
I am using wall height 1 kcal/mol and sigma as 0.3 

Is it indicating issues with CV?

I will be very grateful if you can give some valuable insights

Thank you

[Quentin Pessemesse]

Hello,

I think you need to provide a bit more information such as the total time of your simulation, the type of system and maybe an input file. Metadynamics can be very long to converge, and it might be impossible to achieve convergence with ab initio depending on your CV for instance. If you observe crossing, it means that your CV is able to "push" the system over the barrier which is already a good sign.

The plumed tutorials and google groups have information on how to understand how to see if a metadynamics simulation is converged, for instance:

https://groups.google.com/g/plumed-users/c/3R9CghiH8EE

https://www.plumed.org/doc-v2.9/user-doc/html/master-_i_s_d_d-2.html

There are other sampling methods that can help you achieve convergence of a free energy barrier if metadynamics does not work for your system.

Take care,

Quentin

[SAURABH SINGH]

Thank you Quentin

My system consists of simple molecules Cl CO2 and H2O

&GLOBAL
  PROJECT 1a
  RUN_TYPE MD
  PRINT_LEVEL MEDIUM
&END GLOBAL

#&EXT_RESTART
 # RESTART_FILE_NAME 1a-1.restart
#&END EXT_RESTART

&MOTION
  &MD
    ENSEMBLE NVT # ENSEMBLE NVE
    STEPS 50000
    TIMESTEP 1
    TEMPERATURE 300.0
    TEMP_TOL   100
    &THERMOSTAT
      &NOSE
        LENGTH 3 # 3 is the default
        MTS 2 # 2 is the default
        TIMECON 300.0
        YOSHIDA 3 # 3 is the default
      &END NOSE
    &END THERMOSTAT
  &END MD

  &FREE_ENERGY
    &METADYN
      DO_HILLS
      NT_HILLS 100
      WW [kcalmol] 0.5
      &METAVAR
        SCALE 0.2
        COLVAR 1
        &WALL
          TYPE QUADRATIC
          POSITION [angstrom] 3.0
          &QUADRATIC
             DIRECTION WALL_PLUS
             K [kcalmol] 100.0
          &END
        &END
      &END METAVAR
      &METAVAR
        SCALE 0.2
        COLVAR 2
        &WALL
          TYPE QUADRATIC
          POSITION [angstrom] 3
          &QUADRATIC
             DIRECTION WALL_PLUS
             K [kcalmol] 100.0
          &END
        &END
      &END METAVAR

      &PRINT
        &COLVAR
           COMMON_ITERATION_LEVELS 3
           &EACH
             MD 1
           &END
        &END
        &HILLS
           COMMON_ITERATION_LEVELS 3
           &EACH
             MD 1
           &END
        &END
      &END

    &END METADYN
  &END

  &PRINT
    &RESTART_HISTORY
      &EACH
        MD 100
      &END
    &END RESTART_HISTORY
    &VELOCITIES
      &EACH
        MD 1
      &END
    &END
  # Normal restart file
  &RESTART
    &EACH
      MD 1
    &END
  &END RESTART
  &END PRINT
&END MOTION

&FORCE_EVAL
  METHOD Quickstep              ! Electronic structure method (DFT,...)
  &DFT
    CHARGE = 0
    MULTIPLICITY = 2
    LSD
    BASIS_SET_FILE_NAME  BASIS_MOLOPT.dat
    POTENTIAL_FILE_NAME  GTH_POTENTIALS.dat
    # if WFN restart is required
    WFN_RESTART_FILE_NAME 1a.wfn
    &MGRID
     #CUTOFF 1000
      CUTOFF 360
      NGRIDS 4
      REL_CUTOFF 60
    &END MGRID
    &QS
      METHOD GPW
      EPS_DEFAULT 1.0E-10 # 1.0E-10 is the default
    &END QS
   &SCF
     SCF_GUESS RESTART
     EPS_SCF 1.0E-5 # 1.0E-5 is the default
     MAX_SCF 600 # 50 is the default
     &OT
     MINIMIZER CG
     PRECONDITIONER FULL_ALL
     LINESEARCH 3PNT
     &END
     &OUTER_SCF
       EPS_SCF 1.0E-5 # 1.0E-5 is the default
       MAX_SCF 20 # Keep > 100 when OT is on
     &END
   &END SCF  
   &XC
     &XC_FUNCTIONAL PBE
     &END XC_FUNCTIONAL
     &vdW_POTENTIAL
       DISPERSION_FUNCTIONAL PAIR_POTENTIAL # POTENTIAL_TYPE is alias to DISPERSION_FUNCTIONAL
       &PAIR_POTENTIAL
         TYPE DFTD3
         PARAMETER_FILE_NAME dftd3.dat
         REFERENCE_FUNCTIONAL PBE
       &END PAIR_POTENTIAL
    &END vdW_POTENTIAL
   &END XC
   &POISSON                    ! Solver requested for non periodic calculations
      PERIODIC NONE        ! Type of solver
      POISSON_SOLVER MT
   &END POISSON
  &END DFT  
  &SUBSYS
    &CELL
  A     20.0000000     0.00000000     0.00000000
  B     0.00000000     20.0000000     0.00000000
  C     0.00000000      0.00000000    20.00000000

      PERIODIC NONE              ! Non periodic calculations. That's why the POISSON section is needed
    &END CELL
    &TOPOLOGY                    ! Section used to center the atomic
      COORD_FILE_FORMAT xyz
      COORD_FILE_NAME  1a.xyz
      &CENTER_COORDINATES
      &END
    &END
    &KIND C
      ELEMENT C
      BASIS_SET DZVP-MOLOPT-SR-GTH-q4
      POTENTIAL GTH-PBE-q4
    &END KIND    
    &KIND H
      ELEMENT H
      BASIS_SET DZVP-MOLOPT-SR-GTH-q1
      POTENTIAL GTH-PBE-q1
    &END KIND    
    &KIND Cl
      ELEMENT Cl
      BASIS_SET DZVP-MOLOPT-SR-GTH-q7
      POTENTIAL GTH-PBE-q7
    &END KIND
    &KIND O
      ELEMENT O
      BASIS_SET DZVP-MOLOPT-SR-GTH-q6
      POTENTIAL GTH-PBE-q6
    &END KIND

    &COLVAR
       &DISTANCE
          ATOMS 5 3
       &END DISTANCE
    &END COLVAR
    &COLVAR
       &DISTANCE
          ATOMS 7 3
       &END DISTANCE
    &END COLVAR
  &END SUBSYS
&END FORCE_EVAL


I have attached my input file

[Quentin Pessemesse]

In my opinion with small hills deposited every 100 fs, you’re unlikely to converge the Metadynamics bias. 

You should have a look at J. Phys. Chem. B 2005, 109, 14, 6714–6721 https://doi.org/10.1021/jp045424k, on choosing the optimal parameters for ab initio metadynamics. However, even with increased bias deposition rate, it is possible that meta dynamics bias convergence could take more than hundreds or thousands of ps depending on the system, although it’s hard to say. Maybe you can have a qualitative intepretation if you are close to convergence and do not require accurate barriers, and rather need orders of magnitude.

Take care :) 

Q.