<div dir="ltr">The way I was doing is not appropriate.<div><br></div><div>The problem I faced is like this. For example</div><div><br></div><div>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.</div><div><br></div><div>I post the correct way I thought to deal with this problem here in case someone faced the same problem.</div><div>=====</div><div><div>  SUBROUTINE force_imp_torsions (id_type,s32, is32, ism, isn,  dist1, dist2, tm,&</div><div>       tn, t12, t32, k, phi0, gt1, gt2, gt3, gt4, energy, fscalar )</div><div>    INTEGER, INTENT(IN)                      :: id_type</div><div>    REAL(KIND=dp), INTENT(IN)                :: s32, is32, ism, isn, dist1, &</div><div>                                                dist2</div><div>    REAL(KIND=dp), DIMENSION(:), INTENT(IN)  :: tm, tn, t12, t32</div><div>    REAL(KIND=dp), INTENT(IN)                :: k, phi0</div><div>    REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gt1, gt2, gt3, gt4</div><div>    REAL(KIND=dp), INTENT(OUT)               :: energy, fscalar</div><div><br></div><div>    CHARACTER(len=*), PARAMETER :: routineN = 'force_imp_torsions', &</div><div>      routineP = moduleN//':'//routineN</div><div>    REAL(KIND=dp), PARAMETER                 :: f12 = 1.0_dp/2.0_dp</div><div><br></div><div>    LOGICAL                                  :: failure</div><div>    REAL(KIND=dp)                            :: cosphi, phi ,pi, phi1</div><div>    pi  = 3.14159265358979323846264338_dp</div><div><br></div><div>    failure = .FALSE.</div><div>    ! define cosphi</div><div>    cosphi = DOT_PRODUCT ( tm, tn ) * ism * isn</div><div>    IF(cosphi >  1.0_dp) cosphi = 1.0_dp</div><div>    IF(cosphi < -1.0_dp) cosphi = -1.0_dp</div><div>    phi = SIGN( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn ))</div><div>    !phi = ACOS ( cosphi )</div><div>    phi1 = phi0</div><div>    IF ( phi - phi0 < -pi ) phi1 = phi0 - 2.0_dp*pi</div><div>    IF ( phi - phi0 > pi ) phi1 = phi0 + 2.0_dp*pi</div><div>    !write(*,*) phi, phi0, phi1</div></div><div><br></div><div><div>    SELECT CASE (id_type)</div><div>    CASE (do_ff_charmm)</div><div>       ! compute energy</div><div>       energy = k * ( phi - phi1 )**2</div><div><br></div><div>       ! compute fscalar</div><div>       fscalar = - 2.0_dp * k * ( phi - phi1 )</div><div><br></div><div>    CASE (do_ff_harmonic,do_ff_g87,do_ff_g96)</div><div>       ! compute energy</div><div>       energy = f12 * k * ( phi - phi1 )**2</div><div><br></div><div>       ! compute fscalar</div><div>       fscalar = - k * ( phi - phi1 )</div><div><br></div><div>    CASE DEFAULT</div><div>       CALL stop_program(routineN,moduleN,__LINE__,"Unmatched improper kind")</div><div>    END SELECT</div><div><br></div><div>    ! compute the gradients</div><div>    gt1 =  ( s32 * ism * ism )   * tm</div><div>    gt4 = -( s32 * isn * isn )   * tn</div><div>    gt2 =  ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 - dist2 * is32 ** 2 * gt4</div><div>    gt3 =  ( dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1</div><div>  END SUBROUTINE force_imp_torsions</div></div><div>======</div><div><br>On Friday, August 14, 2015 at 11:28:30 AM UTC-4, 姚懿 wrote:<blockquote class="gmail_quote" style="margin: 0;margin-left: 0.8ex;border-left: 1px #ccc solid;padding-left: 1ex;"><div dir="ltr">In the mol_force.F file. there is one line in force_imp_torsions subroutine<div><br><div>phi = SIGN( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn ))<br></div></div><div><br></div><div>I suggest to change it to </div><div><br></div><div>phi = ACOS( cosphi )</div><div><br></div><div>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.</div><div><br></div><div>YY</div></div></blockquote></div></div>