[gmx-users] tabulated dihedral potential

[Parvez Mh]

(Repost)
Dear all:

I am trying to use tabulated dihedral potential in my simulation. As a test
case, i tried with something known case like octane. I know the
Ryckaert-Bellemans function coefficient of dihedral. So i tried to create a
table from this known case using following fortran code

program table
    implicit none
    double precision::a0,a1,a2,a3,a4,a5
    double precision::x
    double precision::fx,dfx
    character(len=15)::filename
    integer::angle_deg
    a0=8.3945
    a1=16.7787
    a2=1.1332
    a3=-26.3070
    a4=0.000676
    a5=-0.00001816
    filename='table_d0.xvg'
    open(unit=100,file=filename,action='write',status='replace')

    do angle_deg=-180,180,1

        x=angle_deg*0.0174533    !converting degree to radian

        fx=a0-a1*cos(x)+a2*(cos(x))**2-a3*(cos(x))**3+a4*(cos(x))**4 &
        -a5*(cos(x))**5

        dfx=sin(x)*(a1-2*a2*cos(x)+3*a3*(cos(x))**2-4*a4*(cos(x))**3+ &
        5*a5*(cos(x))**4)     !analytical differentiation of
Ryckaert-Bellemans function

        write(100,*)angle_deg,fx,-dfx
    end do
end program table

I used force constant equal to 1 in topology file.
In npt equlibration, my simulation crashed with lot of LINCS error.
Would somebody help me figure out where i did the mistake.

Regards
Masrul