problems of improper dihedral calculation with Fist

[姚懿]

The way I was doing is not appropriate.

The problem I faced is like this. For example

an improper dihedral angle of 178 degree should be very close to -178 
degree. But in the implementation of cp2k, the difference between them is 
356 degree.

I post the correct way I thought to deal with this problem here in case 
someone faced the same problem.
=====
  SUBROUTINE force_imp_torsions (id_type,s32, is32, ism, isn,  dist1, 
dist2, tm,&
       tn, t12, t32, k, phi0, gt1, gt2, gt3, gt4, energy, fscalar )
    INTEGER, INTENT(IN)                      :: id_type
    REAL(KIND=dp), INTENT(IN)                :: s32, is32, ism, isn, dist1, 
&
                                                dist2
    REAL(KIND=dp), DIMENSION(:), INTENT(IN)  :: tm, tn, t12, t32
    REAL(KIND=dp), INTENT(IN)                :: k, phi0
    REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gt1, gt2, gt3, gt4
    REAL(KIND=dp), INTENT(OUT)               :: energy, fscalar

    CHARACTER(len=*), PARAMETER :: routineN = 'force_imp_torsions', &
      routineP = moduleN//':'//routineN
    REAL(KIND=dp), PARAMETER                 :: f12 = 1.0_dp/2.0_dp

    LOGICAL                                  :: failure
    REAL(KIND=dp)                            :: cosphi, phi ,pi, phi1
    pi  = 3.14159265358979323846264338_dp

    failure = .FALSE.
    ! define cosphi
    cosphi = DOT_PRODUCT ( tm, tn ) * ism * isn
    IF(cosphi >  1.0_dp) cosphi = 1.0_dp
    IF(cosphi < -1.0_dp) cosphi = -1.0_dp
    phi = SIGN( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn ))
    !phi = ACOS ( cosphi )
    phi1 = phi0
    IF ( phi - phi0 < -pi ) phi1 = phi0 - 2.0_dp*pi
    IF ( phi - phi0 > pi ) phi1 = phi0 + 2.0_dp*pi
    !write(*,*) phi, phi0, phi1

    SELECT CASE (id_type)
    CASE (do_ff_charmm)
       ! compute energy
       energy = k * ( phi - phi1 )**2

       ! compute fscalar
       fscalar = - 2.0_dp * k * ( phi - phi1 )

    CASE (do_ff_harmonic,do_ff_g87,do_ff_g96)
       ! compute energy
       energy = f12 * k * ( phi - phi1 )**2

       ! compute fscalar
       fscalar = - k * ( phi - phi1 )

    CASE DEFAULT
       CALL stop_program(routineN,moduleN,__LINE__,"Unmatched improper 
kind")
    END SELECT

    ! compute the gradients
    gt1 =  ( s32 * ism * ism )   * tm
    gt4 = -( s32 * isn * isn )   * tn
    gt2 =  ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 - dist2 * is32 ** 2 * gt4
    gt3 =  ( dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1
  END SUBROUTINE force_imp_torsions
======

On Friday, August 14, 2015 at 11:28:30 AM UTC-4, 姚懿 wrote:
>
> In the mol_force.F file. there is one line in force_imp_torsions subroutine
>
> phi = SIGN( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn ))
>
> I suggest to change it to 
>
> phi = ACOS( cosphi )
>
> since improper potential is not periodic. For a dihedral close to 180 
> degree it can varies between -180.0 and 180.0 that will influence the 
> result a lot.
>
> YY
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20150814/3a0e4ec7/attachment.htm>