1      real*8 function dia_torsion(x,nx,w,nw,na,i,i2,j,j2,k,k2,l,l2)
2c
3c $Id$
4c
5      implicit none
6#include "global.fh"
7c
8      integer nx,nw,na,i,i2,j,j2,k,k2,l,l2
9      real*8 x(nx,3),w(nw,na,3)
10c
11      real*8 xijx,xkjx,xklx,xikx,xjlx
12      real*8 xijy,xkjy,xkly,xiky,xjly
13      real*8 xijz,xkjz,xklz,xikz,xjlz
14      real*8 xmx,xmy,xmz,xnx,xny,xnz,rm2i,rn2i,rmni,cphi,s,phi
15c
16      real*8 xa(3),xb(3),xc(3),xd(3)
17c
18      if(i2.le.0) then
19      xa(1)=x(i,1)
20      xa(2)=x(i,2)
21      xa(3)=x(i,3)
22      else
23      xa(1)=w(i,i2,1)
24      xa(2)=w(i,i2,2)
25      xa(3)=w(i,i2,3)
26      endif
27      if(j2.le.0) then
28      xb(1)=x(j,1)
29      xb(2)=x(j,2)
30      xb(3)=x(j,3)
31      else
32      xb(1)=w(j,j2,1)
33      xb(2)=w(j,j2,2)
34      xb(3)=w(j,j2,3)
35      endif
36      if(k2.le.0) then
37      xc(1)=x(k,1)
38      xc(2)=x(k,2)
39      xc(3)=x(k,3)
40      else
41      xc(1)=w(k,k2,1)
42      xc(2)=w(k,k2,2)
43      xc(3)=w(k,k2,3)
44      endif
45      if(l2.le.0) then
46      xd(1)=x(l,1)
47      xd(2)=x(l,2)
48      xd(3)=x(l,3)
49      else
50      xd(1)=w(l,l2,1)
51      xd(2)=w(l,l2,2)
52      xd(3)=w(l,l2,3)
53      endif
54c
55c     determine the dihedral angle
56c     ----------------------------
57c
58      xijx=xa(1)-xb(1)
59      xkjx=xc(1)-xb(1)
60      xklx=xc(1)-xd(1)
61      xikx=xijx-xkjx
62      xjlx=xklx-xkjx
63      xijy=xa(2)-xb(2)
64      xkjy=xc(2)-xb(2)
65      xkly=xc(2)-xd(2)
66      xiky=xijy-xkjy
67      xjly=xkly-xkjy
68      xijz=xa(3)-xb(3)
69      xkjz=xc(3)-xb(3)
70      xklz=xc(3)-xd(3)
71      xikz=xijz-xkjz
72      xjlz=xklz-xkjz
73      xmx=xijy*xkjz-xkjy*xijz
74      xmy=xijz*xkjx-xkjz*xijx
75      xmz=xijx*xkjy-xkjx*xijy
76      xnx=xkjy*xklz-xkly*xkjz
77      xny=xkjz*xklx-xklz*xkjx
78      xnz=xkjx*xkly-xklx*xkjy
79      rm2i=1.0d0/(xmx*xmx+xmy*xmy+xmz*xmz)
80      rn2i=1.0d0/(xnx*xnx+xny*xny+xnz*xnz)
81      rmni=sqrt(rm2i*rn2i)
82      cphi=(xmx*xnx+xmy*xny+xmz*xnz)*rmni
83      if(cphi.lt.-1.0d0) cphi=-1.0d0
84      if(cphi.gt. 1.0d0) cphi= 1.0d0
85      phi=acos(cphi)
86      s=xkjx*(xmy*xnz-xmz*xny) +xkjy*(xmz*xnx-xmx*xnz)
87     + +xkjz*(xmx*xny-xmy*xnx)
88      if(s.lt.0.0d0) phi=-phi
89c
90      dia_torsion=phi
91c
92      return
93      end
94