# problems of improper dihedral calculation with Fist

Fri Aug 14 23:08:02 UTC 2015

```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

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>
```