problems of improper dihedral calculation with Fist
姚懿
yao... at gmail.com
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
! 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>
More information about the CP2K-user
mailing list