Bug in &CONSTRAINT&COLLECTIVE&RESTRAINT / COLVAR?

[caozhe... at gmail.com]

P.S.
The force should be equally distributed onto each atom according to the 
colvar_types.f

  SUBROUTINE eval_point_der(points, i, dsdr, f)
    TYPE(point_type), DIMENSION(:), POINTER  :: points
    INTEGER, INTENT(IN)                      :: i
    REAL(KIND=dp), DIMENSION(:, :), POINTER  :: dsdr
    REAL(KIND=dp), DIMENSION(3), INTENT(IN)  :: f

    INTEGER                                  :: ind, j
    REAL(KIND=dp)                            :: fac

    SELECT CASE(points(i)%type_id)
    CASE(do_clv_geo_center)
       ind = 0
       DO j = 1, i-1
          IF (ASSOCIATED(points(j)%atoms)) THEN
             ind = ind + SIZE(points(j)%atoms)
          END IF
       END DO
       DO j = 1,SIZE(points(i)%atoms)
          fac = points(i)%weights(j)
          dsdr(:,ind+j) = dsdr(:,ind+j) + f * fac
       END DO
    CASE(do_clv_fix_point )
       ! Do nothing if it's a fixed point in space
    END SELECT

  END SUBROUTINE eval_point_der

the force is weighted by the weighting factors :
         fac = points(i)%weights(j)
          dsdr(:,ind+j) = dsdr(:,ind+j) + f * fac









On Sunday, November 1, 2015 at 2:05:43 PM UTC+3, caoz... at gmail.com wrote:
>
> 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/49fa2d5a/attachment.htm>