[CP2K-user] [CP2K:16611] TIP4P/2005f implementation

[mayank...@gmail.com]

Hi,

I was trying to implement TIP4P/2005f force field in CP2K. I have setup the 
relevant forcefield and geometry section in the following manner:

&FORCE_EVAL
  METHOD FIST
  &MM
    &FORCEFIELD
      &SPLINE
        EMAX_SPLINE 1.0E+128
        RCUT_NB 12.0
      &END SPLINE
      &BEND
        ATOMS      HW     OW     HW
        THETA0     [deg] 107.4
        K          [rad^2kjmol] 367.810
        KIND   G87
      &END BEND
      &BEND
        ATOMS      HW     OW     MW
        THETA0     [deg]  53.7
        !K          [rad^2kjmol]  418.400
        K          [rad^2kjmol]  0.000
        KIND   G87
      &END BEND
      &BEND
        ATOMS      MW     HW     OW
        THETA0     [deg] 107.4
        K          [rad^2kjmol] 0.000
        KIND   G87
      &END BEND
      &BEND
        ATOMS      MW     HW     HW
        THETA0     [deg] 107.4
        K          [rad^2kjmol] 0.000
        KIND   G87
      &END BEND
      &BEND
        ATOMS      OW     HW     HW
        THETA0     [deg] 107.4
        K          [rad^2kjmol] 0.000
        KIND   G87
      &END BEND
      &BEND
        ATOMS      OW     MW     HW
        THETA0     [deg] 107.4
        K          [rad^2kjmol] 0.000
        KIND   G87
      &END BEND
      &BEND
        ATOMS      HW     MW     HW
        THETA0     [deg] 107.4
        K          [rad^2kjmol] 0.000
        KIND   G87
      &END BEND
      &BOND
        ATOMS OW HW
        K   [kjmol]     432.581  [nm^-1]   22.87
        R0     [nm] 0.09419
        KIND   MORSE
      &END BOND
      &BOND
        ATOMS OW MW
        K      [nm^-2kjmol]  753120.0
        !K      [nm^-2kjmol]  00.0
        R0     [nm]  0.01750
        KIND   G87
      &END BOND
      &BOND
        ATOMS HW HW
        K      [nm^-2kjmol]  000000.0
        R0     [nm]  0.01750
        KIND   G87
      &END BOND
      &BOND
        ATOMS HW MW
        K      [nm^-2kjmol]  000000.0
        R0     [nm]  0.01750
        KIND   G87
      &END BOND
      &NONBONDED
        &LENNARD-JONES
          ATOMS      OW OW
          EPSILON    [kjmol] 0.77490
          SIGMA      [nm] 3.1644E-1
        &END LENNARD-JONES
        &LENNARD-JONES
          ATOMS      OW MW
          EPSILON    0.0
          SIGMA      0.1
        &END LENNARD-JONES
        &LENNARD-JONES
          ATOMS      OW HW
          EPSILON    0.0
          SIGMA      0.1
        &END LENNARD-JONES
        &LENNARD-JONES
          ATOMS      HW HW
          EPSILON    0.0
          SIGMA      0.1
        &END LENNARD-JONES
        &LENNARD-JONES
          ATOMS      HW MW
          EPSILON    0.0
          SIGMA      0.1
        &END LENNARD-JONES
        &LENNARD-JONES
          ATOMS      MW MW
          EPSILON    0.0
          SIGMA      0.1
        &END LENNARD-JONES
      &END NONBONDED
      &CHARGE
        ATOM OW
        CHARGE 0.000000
      &END CHARGE
      &CHARGE
        ATOM HW
        CHARGE 0.52
      &END CHARGE
      &CHARGE
        ATOM MW
        CHARGE -1.04
      &END CHARGE
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE SPME
        ALPHA .44
        GMAX 11
        O_SPLINE 4
      &END EWALD
    &END POISSON
  &END MM
  &SUBSYS
    &CELL
       ABC 8 8 8
    &END CELL
    &COORD
    OW    2.4532128864011757E+01    2.9058620209822454E+01   
 3.8939754250319247E+01 H2O 1
    HW    2.5226244482849356E+01    2.9202964725171807E+01   
 3.8296636945801339E+01 H2O 1
    HW    2.4993856984368403E+01    2.8952666650621637E+01   
 3.9771507029174259E+01 H2O 1
    MW    2.4684627509816810E+01    2.9063685400405998E+01   
 3.8964642263648159E+01 H2O 1
    &END COORD
    &COLVAR
      &ANGLE
        ATOMS 1 2 3
      &END ANGLE
    &END COLVAR

Additionally I am using VIRTUAL_ATOMS section for the TIP4P virtual site, 
and I had fixed the H-O-H angle to cross check the implementation. When I 
ran the simulation I observed that the O-H distance seemed to vary between 
0.6 to 1.7 angstroms which seems a bit too high for water NVE run at 300K. 
If the remove the H-O-H angle constraint, the water molecule completely 
dissociates. Does anyone have any suggestions or ideas on how to 
effectively setup this forcefield model in CP2K, or what should I correct 
in my input file?

Best Regards

-- 
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+unsubscribe at googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/22d4e3fc-4ae6-4dd4-894d-feb82860057an%40googlegroups.com.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20220217/fc07620d/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test2.inp
Type: chemical/x-gamess-input
Size: 5563 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20220217/fc07620d/attachment.inp>