Fixing a bond in DFT

[Pankaj Mishra]

Dear All,

I am wondering whether it is possible to fix  a particular bond distance in molecule while using quickstep.
For example: I have a H2O molecule and i want to fix 1 OH bond length to some particular value. Then treat the molecule with DFT (rather than MM with force field).

Please suggest if you have some idea.

Thank you so much

[Marco]

Hello,

You can use the FIXED_ATOMS functionality. See link below:

http://manual.cp2k.org/cp2k-2_4-branch/CP2K_INPUT/MOTION/CONSTRAINT/FIXED_ATOMS.html

Regards,
Marco

[Ari Paavo Seitsonen]

Dear Pankaj Mishra,

  You can do it during molecular dynamics using section

CP2K_INPUT / MOTION / CONSTRAINT / COLLECTIVE

and specifying it in

CP2K_INPUT / FORCE_EVAL / SUBSYS / COLVAR

  During geometry optimisation ("GEO_OPT") it is not possible to constraint a bond exactly, but one can work with RESTRAIN and make the condition be _almost_ fulfilled.

  The FIXED_ATOMS works only if the bond to be fixed is along one Cartesian axis and the specific coordinate along the direction of the bond can be fixed; this is usually not possible, but it depends on your system.

    Greetings from Zurich,

       apsi

[Pankaj Mishra]

Dear Ari and Marco,

Thank you so much. It was helpful.

with regards--

Pankaj

[Wei]

Dear Pankaj, Ari and Marco,

I'd also like to fix a bond during a MD (DFT) run (not the absolute atomic positions). I tried to specify CP2K_INPUT / MOTION / CONSTRAINT / COLLECTIVE and CP2K_INPUT / FORCE_EVAL / SUBSYS / COLVAR by following an expample from regtest /Fist/regtest-12/water_3_dist_2.inp.

However, the calculation always crashes with an error:

 *
 *** 12:16:30 ERRORL2 in topology_input:read_constraints_section processor  ***
 ***      0  err=-300  condition FAILED at line 561                         ***
 *

 ===== Routine Calling Stack =====

            3 topology_control
            2 quickstep_create_force_env
            1 CP2K
 CP2K| condition FAILED at line 561
 CP2K| Abnormal program termination, stopped by process number 0


Could you please give me some hint on this issue? I put my input file below. Thanks in advance!

Best Regards

Wei

input

&FORCE_EVAL
  METHOD QUICKSTEP
  &DFT
    BASIS_SET_FILE_NAME ./GST.qbs
    POTENTIAL_FILE_NAME ./GTH_POTENTIALS
    &MGRID
      CUTOFF 300
    &END MGRID
    &QS
      EPS_DEFAULT 1.0e-12
      EXTRAPOLATION ASPC
      EXTRAPOLATION_ORDER 0
      MAP_CONSISTENT
    &END QS
    &SCF
      EPS_SCF 1.0e-7
      MAX_SCF 40
      MAX_SCF_HIST 2
      &OT ON
        MINIMIZER DIIS
        PRECONDITIONER FULL_KINETIC
        ENERGY_GAP 0.002
        STEPSIZE 0.0875
      &END OT
      &OUTER_SCF ON
        MAX_SCF 800
        EPS_SCF 1.0e-7
      &END OUTER_SCF
      SCF_GUESS RANDOM
#     SCF_GUESS HISTORY_RESTART
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC   18.62463  18.62463  18.62463
    &END CELL
    &COORD
 Ge 10.0 10.0 10.0  bond
 Te 13.0 10.0 10.0  bond
    &END COORD
    &COLVAR
      &DISTANCE
        ATOMS 1 2
      &END DISTANCE
    &END COLVAR
    &KIND Ge
      BASIS_SET DZVP-MOLOPT-SR-GTH-q4
      POTENTIAL GTH-PBE-q4
    &END KIND
    &KIND Te
      BASIS_SET DZVP-MOLOPT-SR-GTH-q6
      POTENTIAL GTH-PBE-q6
    &END KIND
  &END SUBSYS
&END FORCE_EVAL
&GLOBAL
  PREFERRED_FFT_LIBRARY FFTW
  PRINT_LEVEL LOW
  PROJECT GeTe.md
  RUN_TYPE MD
&END GLOBAL

&MOTION
  &MD
    ENSEMBLE LANGEVIN
    STEPS    20
    TEMPERATURE  1100
    TIMESTEP 2.0
    &LANGEVIN
      GAMMA 4.0E-3
      NOISY_GAMMA 3.5E-5
    &END LANGEVIN
  &END MD
  &CONSTRAINT
    &COLLECTIVE
      MOLNAME  bond
      COLVAR 1
    &END COLLECTIVE
    &COLLECTIVE
      MOLNAME  bond
      COLVAR 2
    &END COLLECTIVE
  &END CONSTRAINT
&END MOTION

[Andreas Funk]

Dear Wei,

The first thing I see in your input is that you specify one colvar in SUBSYS, but in MOTION you call two colvars. This cannot work. As it is only one bond you specify as colvar, you should work with one colvar in the CONSTRAINT section. Then you should keep it "fixed" like this in section MOTION:

&CONSTRAINT
  &COLLECTIVE
    COLVAR 1
    INTERMOLECULAR
    TARGET 2.399706
    &RESTRAINT
      K 1.0
    &END RESTRAINT
  &END COLLECTIVE
&END CONSTRAINT

Please keep in mind that the distance you want it to be fixed (TARGET) is in bohrs, not angstroms. You need to experiment a little with the force constant for the restraint. Check it in your runs if it is sufficient to keep the bond length fixed within your desired accuracy. In principle, it should do the job. I have used this kind of constraint successfully for both, geometry optimizations and MD.

Best regards

Andreas

[Wei]

Dear Andreas,

thanks a lot for your fast reply! Your suggestion helps a lot. Without RESTRAIN section, the bond is strictly fixed at the target value (my purpose).

If I turn on the RESTRAINT section, the bond length follows a harmonic oscillation, but strength is quite large. Moreover, the temperature fluctuates too wildly, from 300K to 30000K.

Best Regards,

Wei

[Pankaj Mishra]

Dear Andreas, Ari and Marco

Your suggestion is helpful. But i still i am missing one point.
Suppose i am fixing  for larger system (suppose 32 water molecules), then how can i fix the  distance between oxygen and hydrogen of water. I mean do i need to define bond for each bonded atom pair separately for distance in  COLVAR ? Is there some way to define molecule? For DFT calculation, i am just using atom ID (O or H) and its positions in random order.
I hope for some suggestion.

Thank you so much