Bug in &CONSTRAINT&COLLECTIVE&RESTRAINT / COLVAR?

[caozhe... at gmail.com]

Dear CP2K Developers/Users,

There seems to be a wired situation when using CP2K  
&CONSTRAINT&COLLECTIVE&RESTRAINT, which I wonder there might be a bug in 
the code...

I have done several tests to identify the problem, and the tests are listed 
as follows:

1) I used following commands to define the COLVAR, while I have turn off 
all the other interactions such as intra/inter molecular interactions
    &COLVAR
      &DISTANCE_POINT_PLANE
        &POINT
          ATOMS 4
          TYPE GEO_CENTER
        &END POINT
        ATOM_POINT 1
        &POINT
          ATOMS 1
        &END
        &POINT
          ATOMS 2
        &END
        &POINT
          ATOMS 3
        &END
        ATOMS_PLANE 2 3 4
        PBC TRUE
      &END DISTANCE_POINT_PLANE
    &END COLVAR

and the following lines to define the constraint:

  &CONSTRAINT
    &FIXED_ATOMS
      LIST 1..3
    &END FIXED_ATOMS
    &COLLECTIVE
      COLVAR               1
      INTERMOLECULAR  T
      TARGET  [angstrom] 0.
      &RESTRAINT
        K   1.
      &END RESTRAINT
    &END COLLECTIVE
  &END CONSTRAINT

The code can give very good results
i.e.
coordinates:
Pt      0.0     0.0     0.0
Pt      1.0     0.0     0.0
Pt      0.0     1.0     0.0
O       0.0     0.0     1.
forces:
       4
 i =        0, time =        0.000, E =         3.5710648573
 PT         0.0000000000        0.0000000000        0.0000000000
 PT         0.0000000000        0.0000000000        0.0000000000
 PT         0.0000000000        0.0000000000        0.0000000000
  O         0.0000000000        0.0000000000       -3.7794522658
It is very clear that the force comes from -2*k*(distance between point and 
plane)
 
However, the results become wired when I used two particles to define the 
point:
2)
    &COLVAR
      &DISTANCE_POINT_PLANE
        &POINT
          ATOMS 4 5
          TYPE GEO_CENTER
        &END POINT
        ATOM_POINT 1
        &POINT
          ATOMS 1
        &END
        &POINT
          ATOMS 2
        &END
        &POINT
          ATOMS 3
        &END
        ATOMS_PLANE 2 3 4
        PBC TRUE
      &END DISTANCE_POINT_PLANE
    &END COLVAR

same constraint
  &CONSTRAINT
    &FIXED_ATOMS
      LIST 1..3
    &END FIXED_ATOMS
    &COLLECTIVE
      COLVAR               1
      INTERMOLECULAR  T
      TARGET  [angstrom] 0.
      &RESTRAINT
        K   1.
      &END RESTRAINT
    &END COLLECTIVE
  &END CONSTRAINT

i.e.
coordinates:
Pt      0.0     0.0     0.0
Pt      1.0     0.0     0.0
Pt      0.0     1.0     0.0
O       0.0     0.0     1.
O       1.0     0.0     1.

forces:
       5
 i =        0, time =        0.000, E =         3.5710648573
 PT         0.0000000000        0.0000000000        0.0000000000
 PT         0.0000000000        0.0000000000        0.0000000000
 PT         0.0000000000        0.0000000000        0.0000000000
  O         0.0000000000        0.0000000000        0.0000000000
  O         0.0000000000        0.0000000000       -3.7794522658

I expected the total force should be distributed equally onto the 4th and 
5th atoms, while the code allocates all the force onto the 5th atom. In 
this case, at least the total force is correct, though the distribution of 
the total force is wired....

Then, I tried the third testing case:
3)
    &COLVAR
      &DISTANCE_POINT_PLANE
        &POINT
          ATOMS 4 5 6
          TYPE GEO_CENTER
        &END POINT
        ATOM_POINT 1
        &POINT
          ATOMS 1
        &END
        &POINT
          ATOMS 2
        &END
        &POINT
          ATOMS 3
        &END
        ATOMS_PLANE 2 3 4
        PBC TRUE
      &END DISTANCE_POINT_PLANE
    &END COLVAR

with the same constraint:
  &CONSTRAINT
    &FIXED_ATOMS
      LIST 1..3
    &END FIXED_ATOMS
    &COLLECTIVE
      COLVAR               1
      INTERMOLECULAR  T
      TARGET  [angstrom] 0.
      &RESTRAINT
        K   1.
      &END RESTRAINT
    &END COLLECTIVE
  &END CONSTRAINT

The output is hard to be understood...
coordinates:
Pt      0.0     0.0     0.0
Pt      1.0     0.0     0.0
Pt      0.0     1.0     0.0
O       0.0     0.0     1.
O       1.0     0.0     1.
O       0.0     1.0     1.
force:
       6
 i =        0, time =        0.000, E =         3.5710648573
 PT         0.0000000000        0.0000000000        0.0000000000
 PT         0.0000000000        0.0000000000        0.0000000000
 PT         0.0000000000        0.0000000000        0.0000000000
  O         0.0000000000        0.0000000000        1.2598174219
  O         0.0000000000        0.0000000000        1.2598174219
  O         0.0000000000        0.0000000000       -3.7794522658

Though the code gives the exact same potential energy (E), the total force 
is wrong...


Can anyone give a hint what on earth happened: is this a bug of the code or 
there is something wrong with my input file?

Thank you in advance!

Best Regards,
Cao, Zhen
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20151101/61c57c98/attachment.htm>