1 | ! **************************************************************
|
---|
2 | !
|
---|
3 | ! This file contains the subroutines: dihedr,valang
|
---|
4 | !
|
---|
5 | ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
6 | ! Shura Hayryan, Chin-Ku
|
---|
7 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
8 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
9 | !
|
---|
10 | ! **************************************************************
|
---|
11 |
|
---|
12 | real*8 function dihedr(i1,i2,i3,i4)
|
---|
13 |
|
---|
14 | ! .............................................
|
---|
15 | ! PURPOSE: return dihedral angle (i1,i2,i3,i4)
|
---|
16 | ! [in rad.]
|
---|
17 | !
|
---|
18 | ! INPUT: i1,i2,i3,i4 - indices of four atoms
|
---|
19 | !
|
---|
20 | ! CALLS: none
|
---|
21 | ! .............................................
|
---|
22 |
|
---|
23 | include 'INCL.H'
|
---|
24 |
|
---|
25 | integer i1, i2, i3, i4
|
---|
26 |
|
---|
27 | double precision x1, y1, z1, x2, y2, z2, ux1,uy1, uz1, a, u, u1,
|
---|
28 | & u2, ux2, uy2, uz2
|
---|
29 |
|
---|
30 | x1=xat(i2)-xat(i1)
|
---|
31 | y1=yat(i2)-yat(i1)
|
---|
32 | z1=zat(i2)-zat(i1)
|
---|
33 | x2=xat(i3)-xat(i2)
|
---|
34 | y2=yat(i3)-yat(i2)
|
---|
35 | z2=zat(i3)-zat(i2)
|
---|
36 | ux1=y1*z2-z1*y2
|
---|
37 | uy1=z1*x2-x1*z2
|
---|
38 | uz1=x1*y2-y1*x2
|
---|
39 | x1=xat(i4)-xat(i3)
|
---|
40 | y1=yat(i4)-yat(i3)
|
---|
41 | z1=zat(i4)-zat(i3)
|
---|
42 | ux2=z1*y2-y1*z2
|
---|
43 | uy2=x1*z2-z1*x2
|
---|
44 | uz2=y1*x2-x1*y2
|
---|
45 |
|
---|
46 | u1=ux1*ux1+uy1*uy1+uz1*uz1
|
---|
47 | u2=ux2*ux2+uy2*uy2+uz2*uz2
|
---|
48 | u=u1*u2
|
---|
49 |
|
---|
50 | if (u.ne.zero) then
|
---|
51 | a=(ux1*ux2+uy1*uy2+uz1*uz2)/sqrt(u)
|
---|
52 | a=max(a,-one)
|
---|
53 | a=min(a,one)
|
---|
54 | dihedr=acos(a)
|
---|
55 | if (ux1*(uy2*z2-uz2*y2)+uy1*(uz2*x2-ux2*z2)+
|
---|
56 | & uz1*(ux2*y2-uy2*x2).lt.zero) dihedr =-dihedr
|
---|
57 | return
|
---|
58 | else
|
---|
59 | write (*,'(a,4i5)')' dihedr> Error in coordinates of atoms #: '
|
---|
60 | & ,i1,i2,i3,i4
|
---|
61 |
|
---|
62 | write (*,*) 'stored coordinates are xvals :',
|
---|
63 | & xat(i1),xat(i2),xat(i3),xat(i4)
|
---|
64 | write (*,*) 'yvals:', yat(i1),yat(i2),yat(i3),yat(i4)
|
---|
65 | write (*,*) 'zvals:', zat(i1),zat(i2),zat(i3),zat(i4)
|
---|
66 | call outvar(0,'crash.var')
|
---|
67 | stop
|
---|
68 | endif
|
---|
69 |
|
---|
70 | end
|
---|
71 | ! ************************************
|
---|
72 | real*8 function valang(i1,i2,i3)
|
---|
73 |
|
---|
74 | ! .........................................
|
---|
75 | ! PURPOSE: return valence angle (i1,i2,i3)
|
---|
76 | ! [in rad.] with 'i2' as vertex
|
---|
77 | !
|
---|
78 | ! INPUT: i1,i2,i3 - indices of 3 atoms
|
---|
79 | !
|
---|
80 | ! CALLS: none
|
---|
81 | ! .............................................
|
---|
82 |
|
---|
83 | include 'INCL.H'
|
---|
84 |
|
---|
85 | integer i1, i2, i3
|
---|
86 |
|
---|
87 | double precision h1, h2, h3, x1, x2, x3, y1, y2, y3, x, y, u, a
|
---|
88 |
|
---|
89 | h1=xat(i2)
|
---|
90 | h2=yat(i2)
|
---|
91 | h3=zat(i2)
|
---|
92 | x1=xat(i1)-h1
|
---|
93 | x2=yat(i1)-h2
|
---|
94 | x3=zat(i1)-h3
|
---|
95 | y1=xat(i3)-h1
|
---|
96 | y2=yat(i3)-h2
|
---|
97 | y3=zat(i3)-h3
|
---|
98 |
|
---|
99 | x=x1*x1+x2*x2+x3*x3
|
---|
100 | y=y1*y1+y2*y2+y3*y3
|
---|
101 | u=x*y
|
---|
102 |
|
---|
103 | if (u.ne.zero) then
|
---|
104 |
|
---|
105 | a=(x1*y1+x2*y2+x3*y3)/sqrt(u)
|
---|
106 | a=max(a,-one)
|
---|
107 | a=min(a,one)
|
---|
108 | valang=acos(a)
|
---|
109 | return
|
---|
110 |
|
---|
111 | else
|
---|
112 | write (*,'(a,3i5)')' valang> Error in coordinates of atoms #: '
|
---|
113 | & ,i1,i2,i3
|
---|
114 | write (*,*) 'stored coordinates are xvals :',
|
---|
115 | & xat(i1),xat(i2),xat(i3)
|
---|
116 | write (*,*) 'yvals:', yat(i1),yat(i2),yat(i3)
|
---|
117 | write (*,*) 'zvals:', zat(i1),zat(i2),zat(i3)
|
---|
118 | call outvar(0,'crash.var')
|
---|
119 | stop
|
---|
120 | endif
|
---|
121 |
|
---|
122 | end
|
---|
123 |
|
---|