1      subroutine rotate(v,w,g,x,y)
2c
3c $Id$
4c
5      implicit none
6c
7      real*8 v(3),w(3),x(3),y(3),xx(3),t(3,3)
8      real*8 small,pi,r,a,b,ca,cb,cg,sa,sb,sg,g
9      integer i
10      parameter (small=1.0d-24)
11c
12c     rotation with angle g around vector from v to w
13c     of point x giving result in y
14c
15      if(abs(w(2)-v(2)).lt.small) then
16      if(abs(w(1)-v(1)).lt.small) then
17      a=0.0d0
18      else
19      if(w(1)-v(1).gt.0.0d0) then
20      a=2.0d0*datan(1.0d0)
21      else
22      a=(-2.0d0)*datan(1.0d0)
23      endif
24      endif
25      else
26      a=atan(abs(w(1)-v(1))/abs(w(2)-v(2)))
27      pi=4.0d0*atan(1.0d0)
28      if(w(1)-v(1).gt.0.0d0.and.w(2)-v(2).lt.0.0d0) a=pi-a
29      if(w(1)-v(1).lt.0.0d0.and.w(2)-v(2).gt.0.0d0) a=-a
30      if(w(1)-v(1).lt.0.0d0.and.w(2)-v(2).lt.0.0d0) a=pi+a
31      endif
32      r=0.0d0
33      do 1 i=1,3
34      r=r+(w(i)-v(i))**2
35      xx(i)=x(i)-v(i)
36    1 continue
37      if(r.lt.small) then
38      y(1)=x(1)
39      y(2)=x(2)
40      y(3)=x(3)
41      return
42      endif
43      b=acos((w(3)-v(3))/sqrt(r))
44      sa=sin(a)
45      ca=cos(a)
46      sb=sin(b)
47      cb=cos(b)
48      sg=sin(g)
49      cg=cos(g)
50      t(1,1)=ca*ca*cg-sa*ca*cb*sg+sa*ca*cb*sg
51     +       +sa*sa*cb*cb*cg+sa*sa*sb*sb
52      t(1,2)=(-sa)*ca*cg-ca*ca*cb*sg-sa*sa*cb*sg
53     +       +sa*ca*cb*cb*cg+sa*ca*sb*sb
54      t(1,3)=ca*sb*sg-sa*sb*cb*cg+sa*sb*cb
55      t(2,1)=(-sa)*ca*cg+sa*sa*cb*sg+ca*ca*cb*sg
56     +       +sa*ca*cb*cb*cg+sa*ca*sb*sb
57      t(2,2)=sa*sa*cg+sa*ca*cb*sg-sa*ca*cb*sg
58     +       +ca*ca*cb*cb*cg+ca*ca*sb*sb
59      t(2,3)=(-sa)*sb*sg-ca*sb*cb*cg+ca*sb*cb
60      t(3,1)=(-ca)*sb*sg-sa*sb*cb*cg+sa*sb*cb
61      t(3,2)=sa*sb*sg-ca*sb*cb*cg+ca*sb*cb
62      t(3,3)=sb*sb*cg+cb*cb
63      do 2 i=1,3
64      y(i)=xx(1)*t(i,1)+xx(2)*t(i,2)+xx(3)*t(i,3)+v(i)
65    2 continue
66      return
67      end
68      subroutine super(x,nx,mx,y,w,ws,sdev,ny,my,mod,rms0,rms1,
69     + xw,mw,nw,ma,na,lpara)
70c
71c     superimpose x(1:n,1:3) onto y(1:n,1:3)
72c
73      implicit none
74c
75#include "msgids.fh"
76#include "global.fh"
77c
78      real*8 zero
79      parameter(zero=0.0d0)
80c
81      integer nx,mx,ny,my,mw,nw,ma,na
82      real*8 x(mx,3),y(my,3),w(my),ws(my),sdev(my),xw(mw,ma,3),rms0,rms1
83      logical mod,lpara
84c
85      integer i,j,k,l,nr
86      real*8 u(3,3),c(4,4),q(4),b(4),v(4,4),wnorm,wnorms
87      real*8 xt(3),yt(3),sd
88c
89c      write(*,'(a,6f12.6)') 'Super in ',(x(1,i),i=1,3),(x(nx,i),i=1,3)
90      wnorm=zero
91      wnorms=zero
92      do 1 j=1,3
93      xt(j)=zero
94      yt(j)=zero
95    1 continue
96      do 2 k=1,nx
97      do 3 j=1,3
98      xt(j)=xt(j)+w(k)*x(k,j)
99    3 continue
100    2 continue
101      if(lpara) call ga_dgop(mag_d08,xt,3,'+')
102      do 4 k=1,ny
103      wnorm=wnorm+w(k)
104      wnorms=wnorms+ws(k)
105      do 5 j=1,3
106      yt(j)=yt(j)+w(k)*y(k,j)
107    5 continue
108    4 continue
109c
110      rms0=zero
111      do 6 j=1,3
112      xt(j)=xt(j)/wnorm
113      yt(j)=yt(j)/wnorm
114      do 7 k=1,nx
115      rms0=rms0+ws(k)*(x(k,j)-xt(j)-y(k,j)+yt(j))**2
116    7 continue
117    6 continue
118      if(lpara) call ga_dgop(mag_d09,rms0,1,'+')
119      rms0=sqrt(rms0/wnorms)
120c
121      do 8 i=1,3
122      do 9 j=1,3
123      u(i,j)=zero
124      do 10 k=1,nx
125      u(i,j)=u(i,j)+w(k)*(x(k,i)-xt(i))*(y(k,j)-yt(j))
126   10 continue
127    9 continue
128    8 continue
129      if(lpara) call ga_dgop(mag_d10,u,9,'+')
130c
131      c(1,1)=u(1,1)+u(2,2)+u(3,3)
132      c(1,2)=u(3,2)-u(2,3)
133      c(1,3)=u(1,3)-u(3,1)
134      c(1,4)=u(2,1)-u(1,2)
135      c(2,2)=u(1,1)-u(2,2)-u(3,3)
136      c(2,3)=u(1,2)+u(2,1)
137      c(2,4)=u(3,1)+u(1,3)
138      c(3,3)=u(2,2)-u(3,3)-u(1,1)
139      c(3,4)=u(2,3)+u(3,2)
140      c(4,4)=u(3,3)-u(1,1)-u(2,2)
141c
142      do 11 j=1,3
143      do 12 i=j+1,4
144      c(i,j)=c(j,i)
145   12 continue
146   11 continue
147c
148      call md_jacobi(c,4,4,b,v,nr)
149c
150      do 13 i=1,4
151      q(i)=v(i,4)
152   13 continue
153c
154      u(1,1)=q(1)*q(1)+q(2)*q(2)-q(3)*q(3)-q(4)*q(4)
155      u(1,2)=2.0d0*(q(3)*q(2)+q(1)*q(4))
156      u(1,3)=2.0d0*(q(4)*q(2)-q(1)*q(3))
157      u(2,1)=2.0d0*(q(2)*q(3)-q(1)*q(4))
158      u(2,2)=q(1)*q(1)-q(2)*q(2)+q(3)*q(3)-q(4)*q(4)
159      u(2,3)=2.0d0*(q(4)*q(3)+q(1)*q(2))
160      u(3,1)=2.0d0*(q(2)*q(4)+q(1)*q(3))
161      u(3,2)=2.0d0*(q(3)*q(4)-q(1)*q(2))
162      u(3,3)=q(1)*q(1)-q(2)*q(2)-q(3)*q(3)+q(4)*q(4)
163c
164      rms1=zero
165      do 14 k=1,nx
166      do 15 i=1,3
167      q(i)=0.0d0
168      do 16 j=1,3
169      q(i)=q(i)+u(i,j)*(x(k,j)-xt(j))
170   16 continue
171   15 continue
172      if(mod) then
173      do 17 i=1,3
174      x(k,i)=q(i)+yt(i)
175   17 continue
176      endif
177      sd=0.0d0
178      do 18 j=1,3
179      sd=sd+(q(j)-y(k,j)+yt(j))**2
180   18 continue
181      sdev(k)=sdev(k)+sd
182      rms1=rms1+ws(k)*sd
183   14 continue
184      if(lpara) call ga_dgop(mag_d11,rms1,1,'+')
185      rms1=sqrt(rms1/wnorms)
186c
187      if(mod.and.nw.gt.0) then
188      do 19 l=1,nw
189      do 20 k=1,na
190      do 21 i=1,3
191      q(i)=0.0d0
192      do 22 j=1,3
193      q(i)=q(i)+u(i,j)*(xw(l,k,j)-xt(j))
194   22 continue
195   21 continue
196      do 23 i=1,3
197      xw(l,k,i)=q(i)+yt(i)
198   23 continue
199   20 continue
200   19 continue
201      endif
202c
203c      write(*,'(a,6f12.6)') 'Super out',(x(1,i),i=1,3),(x(nx,i),i=1,3)
204      return
205      end
206      subroutine super2(x,ix,nx,mx,y,w,ws,sdev,ny,my,mod,rms0,rms1,
207     + xw,mw,nw,ma,na,lpara)
208c
209c     superimpose x(1:n,1:3) onto y(1:n,1:3)
210c
211      implicit none
212c
213#include "msgids.fh"
214#include "global.fh"
215c
216      real*8 zero
217      parameter(zero=0.0d0)
218c
219      integer nx,mx,ny,my,mw,nw,ma,na
220      integer ix(mx)
221      real*8 x(mx,3),y(my,3),w(my),ws(my),sdev(my),xw(mw,ma,3),rms0,rms1
222      logical mod,lpara
223c
224      integer i,j,k,l,nr
225      real*8 u(3,3),c(4,4),q(4),b(4),v(4,4),wnorm,wnorms
226      real*8 xt(3),yt(3),sd
227c
228c      write(*,'(a,6f12.6)') 'Super in ',(x(1,i),i=1,3),(x(nx,i),i=1,3)
229      wnorm=zero
230      wnorms=zero
231      do 1 j=1,3
232      xt(j)=zero
233      yt(j)=zero
234    1 continue
235      do 2 k=1,nx
236      do 3 j=1,3
237      xt(j)=xt(j)+w(ix(k))*x(k,j)
238    3 continue
239    2 continue
240      if(lpara) call ga_dgop(mag_d08,xt,3,'+')
241      do 4 k=1,ny
242      wnorm=wnorm+w(k)
243      wnorms=wnorms+ws(k)
244      do 5 j=1,3
245      yt(j)=yt(j)+w(k)*y(k,j)
246    5 continue
247    4 continue
248c
249      rms0=zero
250      do 6 j=1,3
251      xt(j)=xt(j)/wnorm
252      yt(j)=yt(j)/wnorm
253      do 7 k=1,nx
254      rms0=rms0+ws(ix(k))*(x(k,j)-xt(j)-y(ix(k),j)+yt(j))**2
255    7 continue
256    6 continue
257      if(lpara) call ga_dgop(mag_d09,rms0,1,'+')
258      rms0=sqrt(rms0/wnorms)
259c
260      do 8 i=1,3
261      do 9 j=1,3
262      u(i,j)=zero
263      do 10 k=1,nx
264      u(i,j)=u(i,j)+w(ix(k))*(x(k,i)-xt(i))*(y(ix(k),j)-yt(j))
265   10 continue
266    9 continue
267    8 continue
268      if(lpara) call ga_dgop(mag_d10,u,9,'+')
269c
270      c(1,1)=u(1,1)+u(2,2)+u(3,3)
271      c(1,2)=u(3,2)-u(2,3)
272      c(1,3)=u(1,3)-u(3,1)
273      c(1,4)=u(2,1)-u(1,2)
274      c(2,2)=u(1,1)-u(2,2)-u(3,3)
275      c(2,3)=u(1,2)+u(2,1)
276      c(2,4)=u(3,1)+u(1,3)
277      c(3,3)=u(2,2)-u(3,3)-u(1,1)
278      c(3,4)=u(2,3)+u(3,2)
279      c(4,4)=u(3,3)-u(1,1)-u(2,2)
280c
281      do 11 j=1,3
282      do 12 i=j+1,4
283      c(i,j)=c(j,i)
284   12 continue
285   11 continue
286c
287      call md_jacobi(c,4,4,b,v,nr)
288c
289      do 13 i=1,4
290      q(i)=v(i,4)
291   13 continue
292c
293      u(1,1)=q(1)*q(1)+q(2)*q(2)-q(3)*q(3)-q(4)*q(4)
294      u(1,2)=2.0d0*(q(3)*q(2)+q(1)*q(4))
295      u(1,3)=2.0d0*(q(4)*q(2)-q(1)*q(3))
296      u(2,1)=2.0d0*(q(2)*q(3)-q(1)*q(4))
297      u(2,2)=q(1)*q(1)-q(2)*q(2)+q(3)*q(3)-q(4)*q(4)
298      u(2,3)=2.0d0*(q(4)*q(3)+q(1)*q(2))
299      u(3,1)=2.0d0*(q(2)*q(4)+q(1)*q(3))
300      u(3,2)=2.0d0*(q(3)*q(4)-q(1)*q(2))
301      u(3,3)=q(1)*q(1)-q(2)*q(2)-q(3)*q(3)+q(4)*q(4)
302c
303      rms1=zero
304      do 14 k=1,nx
305      do 15 i=1,3
306      q(i)=0.0d0
307      do 16 j=1,3
308      q(i)=q(i)+u(i,j)*(x(k,j)-xt(j))
309   16 continue
310   15 continue
311      if(mod) then
312      do 17 i=1,3
313      x(k,i)=q(i)+yt(i)
314   17 continue
315      endif
316      sd=0.0d0
317      do 18 j=1,3
318      sd=sd+(q(j)-y(ix(k),j)+yt(j))**2
319   18 continue
320      sdev(ix(k))=sdev(ix(k))+sd
321      rms1=rms1+ws(ix(k))*sd
322   14 continue
323      if(lpara) call ga_dgop(mag_d11,rms1,1,'+')
324      rms1=sqrt(rms1/wnorms)
325c
326      if(mod.and.nw.gt.0) then
327      do 19 l=1,nw
328      do 20 k=1,na
329      do 21 i=1,3
330      q(i)=0.0d0
331      do 22 j=1,3
332      q(i)=q(i)+u(i,j)*(xw(l,k,j)-xt(j))
333   22 continue
334   21 continue
335      do 23 i=1,3
336      xw(l,k,i)=q(i)+yt(i)
337   23 continue
338   20 continue
339   19 continue
340      endif
341c
342c      write(*,'(a,6f12.6)') 'Super out',(x(1,i),i=1,3),(x(nx,i),i=1,3)
343      return
344      end
345      real*8 function angl(x,y,z)
346c
347      implicit none
348c
349      real*8 x(3),y(3),z(3)
350c
351      real*8 xy(3),zy(3),rxy,rzy,phi
352      integer i
353c
354      rxy=0.0d0
355      rzy=0.0d0
356      do 1 i=1,3
357      xy(i)=x(i)-y(i)
358      zy(i)=z(i)-y(i)
359    1 continue
360      rxy=xy(1)*xy(1)+xy(2)*xy(2)+xy(3)*xy(3)
361      rzy=zy(1)*zy(1)+zy(2)*zy(2)+zy(3)*zy(3)
362      phi=(xy(1)*zy(1)+xy(2)*zy(2)+xy(3)*zy(3))/sqrt(rxy*rzy)
363      if(phi.lt.-1.0d0) phi=-1.0d0
364      if(phi.gt.1.0d0) phi=1.0d0
365      angl=acos(phi)
366c
367      return
368      end
369      real*8 function atom_radius(number)
370c
371      implicit none
372      integer number
373c
374      real*8 radius(0:105)
375c
376      data radius / 99999.99,
377     + 0.35, 1.22, 1.23, 0.89, 0.88, 0.77, 0.70, 0.66, 0.58, 1.60,
378     + 1.40, 1.36, 1.25, 1.17, 1.10, 1.04, 0.99, 1.91, 2.03, 1.74,
379     + 1.44, 1.32, 1.22, 1.19, 1.17,1.165, 1.16, 1.15, 1.17, 1.25,
380     + 1.25, 1.22, 1.21, 1.17, 1.14, 1.98, 2.22, 1.92, 1.62, 1.45,
381     + 1.34, 1.29, 1.27, 1.24, 1.25, 1.28, 1.34, 1.41, 1.50, 1.40,
382     + 1.41, 1.37, 1.33, 2.09, 2.35, 1.98, 1.69, 1.65, 1.65, 1.64,
383     + 1.65, 1.66, 1.65, 1.61, 1.59, 1.59, 1.58, 1.57, 1.56, 1.56,
384     + 1.56, 1.44, 1.34, 1.30, 1.28, 1.26, 1.26, 1.29, 1.34, 1.44,
385     + 1.55, 1.54, 1.52, 1.53, 1.50, 2.20, 3.24, 2.68, 2.25, 2.16,
386     + 1.93, 1.66, 1.57, 1.81, 2.21, 1.43, 1.42, 1.40, 1.39, 1.38,
387     + 1.37, 1.36, 1.34, 1.30, 1.30 /
388c
389      atom_radius=1.1d-01*radius(number)
390c
391      return
392      end
393      subroutine povinc(lfn,xmin,xmax,ymin,ymax,zmin,zmax)
394c
395      implicit none
396c
397      integer lfn
398      real*8 xmin,xmax,ymin,ymax,zmin,zmax
399c
400      open(unit=lfn,file='camera.inc',form='formatted',status='new',
401     + err=9)
402      write(lfn,1000)
403     + 0.0d0,0.0d0,-1.0d1*max(abs(xmax),abs(ymax),abs(zmax),
404     + abs(xmin),abs(ymin),abs(zmin)),
405     + 5.0d-1*(xmax+xmin),5.0d-1*(ymax+ymin),5.0d-1*(zmax+zmin)
406c     + 0.0d0, 0.0d0, 0.0d0,
407c     + 1.3d0, 0.0d0, 0.0d0,
408c     + 0.0d0, 0.0d0, 1.0d0,
409c     + 0.0d0, 1.0d0, 0.0d0
410 1000 format('camera {',/,
411     + ' location <',f12.6,',',f12.6,',',f12.6,'>',/,
412     + ' look_at <',f12.6,',',f12.6,',',f12.6,'>',/,
413     + ' angle 20.0',/,'}')
414      write(lfn,1001)  0.0d1, 0.0d1,-5.0d1, 1.0d0,1.0d0,1.0d0
415      write(lfn,1001) -1.0d1, 2.0d1,-1.0d1, 1.0d0,1.0d0,1.0d0
416      write(lfn,1001)  1.0d1, 2.0d1,-1.0d1, 1.0d0,1.0d0,1.0d0
417      write(lfn,1001)  0.0d1, 1.0d1,-2.0d1, 1.0d0,1.0d0,1.0d0
418 1001 format('light_source { <',f12.6,',',f12.6,',',f12.6,
419     + '> color rgb <',f4.2,',',f4.2,',',f4.2,'> }')
420      write(lfn,1002) 0.0d0,0.0d0,0.0d0
421 1002 format('background { color rgb <',f4.2,',',f4.2,',',f4.2,'> }')
422      close(unit=lfn)
423c
424    9 continue
425      open(unit=lfn,file='colors.inc',form='formatted',status='new',
426     + err=99)
427      write(lfn,2000)
428 2000 format('#declare Colors_Inc_Temp = version ;',/,'#version 2.0 ;')
429      write(lfn,2001) 'Gray05',0.05,0.05,0.05
430      write(lfn,2001) 'Gray05',0.05,0.05,0.05
431      write(lfn,2001) 'Gray10',0.10,0.10,0.10
432      write(lfn,2001) 'Gray15',0.15,0.15,0.15
433      write(lfn,2001) 'Gray20',0.20,0.20,0.20
434      write(lfn,2001) 'Gray25',0.25,0.25,0.25
435      write(lfn,2001) 'Gray30',0.30,0.30,0.30
436      write(lfn,2001) 'Gray35',0.35,0.35,0.35
437      write(lfn,2001) 'Gray40',0.40,0.40,0.40
438      write(lfn,2001) 'Gray45',0.45,0.45,0.45
439      write(lfn,2001) 'Gray50',0.50,0.50,0.50
440      write(lfn,2001) 'Gray55',0.55,0.55,0.55
441      write(lfn,2001) 'Gray60',0.60,0.60,0.60
442      write(lfn,2001) 'Gray65',0.65,0.65,0.65
443      write(lfn,2001) 'Gray70',0.70,0.70,0.70
444      write(lfn,2001) 'Gray75',0.75,0.75,0.75
445      write(lfn,2001) 'Gray80',0.80,0.80,0.80
446      write(lfn,2001) 'Gray85',0.85,0.85,0.85
447      write(lfn,2001) 'Gray90',0.90,0.90,0.90
448      write(lfn,2001) 'Gray95',0.95,0.95,0.95
449      write(lfn,2001) 'DimGray',0.329412,0.329412,0.329412
450      write(lfn,2001) 'DimGrey',0.329412,0.329412,0.329412
451      write(lfn,2001) 'Gray',0.752941,0.752941,0.752941
452      write(lfn,2001) 'Grey',0.752941,0.752941,0.752941
453      write(lfn,2001) 'LightGray',0.658824,0.658824,0.658824
454      write(lfn,2001) 'LightGrey',0.658824,0.658824,0.658824
455      write(lfn,2001) 'VLightGrey',0.80,0.80,0.80
456      write(lfn,2001) 'White',1.0,1.0,1.0
457      write(lfn,2001) 'Red',1.0,0.0,0.0
458      write(lfn,2001) 'Green',0.0,1.0,0.0
459      write(lfn,2001) 'Blue',0.0,0.0,1.0
460      write(lfn,2001) 'Yellow',1.0,1.0,0.0
461      write(lfn,2001) 'Cyan',0.0,1.0,1.0
462      write(lfn,2001) 'Magenta',1.0,0.0,1.0
463      write(lfn,2001) 'Black',0.0,0.0,0.0
464      write(lfn,2001) 'Aquamarine',0.439216,0.858824,0.576471
465      write(lfn,2001) 'BlueViolet',0.62352,0.372549,0.623529
466      write(lfn,2001) 'Brown',0.647059,0.164706,0.164706
467      write(lfn,2001) 'CadetBlue',0.372549,0.623529,0.623529
468      write(lfn,2001) 'Coral',1.0,0.498039,0.0
469      write(lfn,2001) 'CornflowerBlue',0.258824,0.258824,0.435294
470      write(lfn,2001) 'DarkGreen',0.184314,0.309804,0.184314
471      write(lfn,2001) 'DarkOliveGreen',0.309804,0.309804,0.184314
472      write(lfn,2001) 'DarkOrchid',0.6,0.196078,0.8
473      write(lfn,2001) 'DarkSlateBlue',0.419608,0.137255,0.556863
474      write(lfn,2001) 'DarkSlateGray',0.184314,0.309804,0.309804
475      write(lfn,2001) 'DarkSlateGrey',0.184314,0.309804,0.309804
476      write(lfn,2001) 'DarkTurquoise',0.439216,0.576471,0.858824
477      write(lfn,2001) 'Firebrick',0.556863,0.137255,0.137255
478      write(lfn,2001) 'ForestGreen',0.137255,0.556863,0.137255
479      write(lfn,2001) 'Gold',0.8,0.498039,0.196078
480      write(lfn,2001) 'Goldenrod',0.858824,0.858824,0.439216
481      write(lfn,2001) 'GreenYellow',0.576471,0.858824,0.439216
482      write(lfn,2001) 'IndianRed',0.309804,0.184314,0.184314
483      write(lfn,2001) 'Khaki',0.623529,0.623529,0.372549
484      write(lfn,2001) 'LightBlue',0.74902,0.847059,0.847059
485      write(lfn,2001) 'LightSteelBlue',0.560784,0.560784,0.737255
486      write(lfn,2001) 'LimeGreen',0.196078,0.8,0.196078
487      write(lfn,2001) 'Maroon',0.556863,0.137255,0.419608
488      write(lfn,2001) 'MediumAquamarine',0.196078,0.8,0.6
489      write(lfn,2001) 'MediumBlue',0.196078,0.196078,0.8
490      write(lfn,2001) 'MediumForestGreen',0.419608,0.556863,0.137255
491      write(lfn,2001) 'MediumGoldenrod',0.917647,0.917647,0.678431
492      write(lfn,2001) 'MediumOrchid',0.576471,0.439216,0.858824
493      write(lfn,2001) 'MediumSeaGreen',0.258824,0.435294,0.258824
494      write(lfn,2001) 'MediumSlateBlue',0.498039,1.0,0.0
495      write(lfn,2001) 'MediumSpringGreen',0.498039,1.0,0.0
496      write(lfn,2001) 'MediumTurquoise',0.439216,0.858824,0.858824
497      write(lfn,2001) 'MediumVioletRed',0.858824,0.439216,0.576471
498      write(lfn,2001) 'MidnightBlue',0.184314,0.184314,0.309804
499      write(lfn,2001) 'Navy',0.137255,0.137255,0.556863
500      write(lfn,2001) 'NavyBlue',0.137255,0.137255,0.556863
501      write(lfn,2001) 'Orange',1,0.5,0.0
502      write(lfn,2001) 'OrangeRed',1.0,0.498039,0.0
503      write(lfn,2001) 'Orchid',0.858824,0.439216,0.858824
504      write(lfn,2001) 'PaleGreen',0.560784,0.737255,0.560784
505      write(lfn,2001) 'Pink',0.737255,0.560784,0.560784
506      write(lfn,2001) 'Plum',0.917647,0.678431,0.917647
507      write(lfn,2001) 'Salmon',0.435294,0.258824,0.258824
508      write(lfn,2001) 'SeaGreen',0.137255,0.556863,0.419608
509      write(lfn,2001) 'Sienna',0.556863,0.419608,0.137255
510      write(lfn,2001) 'SkyBlue',0.196078,0.6,0.8
511      write(lfn,2001) 'SlateBlue',0.0,0.498039,1.0
512      write(lfn,2001) 'SpringGreen',0.0,1.0,0.498039
513      write(lfn,2001) 'SteelBlue',0.137255,0.419608,0.556863
514      write(lfn,2001) 'Tan',0.858824,0.576471,0.439216
515      write(lfn,2001) 'Thistle',0.847059,0.74902,0.847059
516      write(lfn,2001) 'Turquoise',0.678431,0.917647,0.917647
517      write(lfn,2001) 'Violet',0.309804,0.184314,0.309804
518      write(lfn,2001) 'VioletRed',0.8,0.196078,0.6
519      write(lfn,2001) 'Wheat',0.847059,0.847059,0.74902
520      write(lfn,2001) 'YellowGreen',0.6,0.8,0.196078
521      write(lfn,2001) 'SummerSky',0.22,0.69,0.87
522      write(lfn,2001) 'RichBlue',0.35,0.35,0.67
523      write(lfn,2001) 'Brass',0.71,0.65,0.26
524      write(lfn,2001) 'Copper',0.72,0.45,0.20
525      write(lfn,2001) 'Bronze',0.55,0.47,0.14
526      write(lfn,2001) 'Bronze2',0.65,0.49,0.24
527      write(lfn,2001) 'Silver',0.90,0.91,0.98
528      write(lfn,2001) 'BrightGold',0.85,0.85,0.10
529      write(lfn,2001) 'OldGold',0.81,0.71,0.23
530      write(lfn,2001) 'Feldspar',0.82,0.57,0.46
531      write(lfn,2001) 'Quartz',0.85,0.85,0.95
532      write(lfn,2001) 'Mica',0.0,0.0,0.0
533      write(lfn,2001) 'NeonPink',1.00,0.43,0.78
534      write(lfn,2001) 'DarkPurple',0.53,0.12,0.47
535      write(lfn,2001) 'NeonBlue',0.30,0.30,1.00
536      write(lfn,2001) 'CoolCopper',0.85,0.53,0.10
537      write(lfn,2001) 'MandarinOrange',0.89,0.47,0.20
538      write(lfn,2001) 'LightWood',0.91,0.76,0.65
539      write(lfn,2001) 'MediumWood',0.65,0.50,0.39
540      write(lfn,2001) 'DarkWood',0.52,0.37,0.26
541      write(lfn,2001) 'SpicyPink',1.00,0.11,0.68
542      write(lfn,2001) 'SemiSweetChoc',0.42,0.26,0.15
543      write(lfn,2001) 'BakersChoc',0.36,0.20,0.09
544      write(lfn,2001) 'Flesh',0.96,0.80,0.69
545      write(lfn,2001) 'NewTan',0.92,0.78,0.62
546      write(lfn,2001) 'NewMidnightBlue',0.00,0.00,0.61
547      write(lfn,2001) 'VeryDarkBrown',0.35,0.16,0.14
548      write(lfn,2001) 'DarkBrown',0.36,0.25,0.20
549      write(lfn,2001) 'DarkTan',0.59,0.41,0.31
550      write(lfn,2001) 'GreenCopper',0.32,0.49,0.46
551      write(lfn,2001) 'DkGreenCopper',0.29,0.46,0.43
552      write(lfn,2001) 'DustyRose',0.52,0.39,0.39
553      write(lfn,2001) 'HuntersGreen',0.13,0.37,0.31
554      write(lfn,2001) 'Scarlet',0.55,0.09,0.09
555      write(lfn,2002) 'Clear',1.0,1.0,1.0,1.0
556 2001 format('#declare ',a,' = color red ',f8.6,' green ',f8.6,
557     + ' blue ',f8.6,' ;')
558 2002 format('#declare ',a,' = color red ',f8.6,' green ',f8.6,
559     + ' blue ',f8.6,' filter ',f8.6,' ;')
560      write(lfn,2003)
561 2003 format('#declare Plane_Map = 0 ;',/,'#declare Sphere_Map = 1 ;',/,
562     + '#declare Cylinder_Map = 2 ;',/,'#declare Torus_Map = 5 ;',/,
563     + '#declare Bi   = 2 ;',/,'#declare Norm = 4 ;',/,
564     + '#version Colors_Inc_Temp ;')
565      close(unit=lfn)
566c
567   99 continue
568      open(unit=lfn,file='plane.inc',form='formatted',status='new',
569     + err=999)
570      write(lfn,3001)
571 3001 format('plane{ <0,1,0>,-1 pigment { White } }')
572      close(unit=lfn)
573c
574  999 continue
575c
576      return
577      end
578      real*8 function torsion(a,b,c,d)
579c
580      implicit none
581c
582      real*8 a(3),b(3),c(3),d(3)
583c
584      torsion=0.0d0
585c
586      return
587      end
588      character*255 function atom_color(number)
589      implicit none
590      integer number
591c
592      integer indexc(0:105)
593c
594      data indexc / 12,
595     +  7,  5, 14, 12, 13,  0,  1,  2,  6, 12,
596     + 13, 15,  9,  6,  8,  3, 13, 12, 16, 16,
597     + 12,  9, 12,  9,  9,  8, 12, 10, 10, 10,
598     + 12, 12, 12, 12, 10, 12, 12, 12, 12, 12,
599     + 12, 12, 12, 12, 12, 12, 16, 12, 12, 12,
600     + 12, 12, 11, 12, 12,  8, 12, 12, 12, 12,
601     + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
602     + 12, 12, 12, 12, 12, 12, 12, 16, 17, 16,
603     + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
604     + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
605     + 12, 12, 12, 12, 12 /
606c
607      atom_color='white '
608      if(indexc(number).eq.0) atom_color='LightGrey '
609      if(indexc(number).eq.1) atom_color='SkyBlue '
610      if(indexc(number).eq.2) atom_color='Red '
611      if(indexc(number).eq.3) atom_color='Yellow '
612      if(indexc(number).eq.4) atom_color='White '
613      if(indexc(number).eq.5) atom_color='Pink '
614      if(indexc(number).eq.6) atom_color='Goldenrod '
615      if(indexc(number).eq.7) atom_color='Blue '
616      if(indexc(number).eq.8) atom_color='Orange '
617      if(indexc(number).eq.9) atom_color='Gray85 '
618      if(indexc(number).eq.10) atom_color='Brown '
619      if(indexc(number).eq.11) atom_color='Purple '
620      if(indexc(number).eq.12) atom_color='SpicyPink '
621      if(indexc(number).eq.13) atom_color='Green '
622      if(indexc(number).eq.14) atom_color='Firebrick '
623      if(indexc(number).eq.15) atom_color='DarkGreen '
624      if(indexc(number).eq.16) atom_color='Silver '
625      if(indexc(number).eq.17) atom_color='Gold '
626c
627      return
628      end
629      subroutine md_jacobi(a,n,na,d,v,nrot)
630c
631c     compute eigenvectors and eigenvalues for real symmetric
632c     matrix using the Jacobi diagonalization
633c
634      implicit none
635c
636      integer nrmax
637      parameter(nrmax=100)
638c
639      real*8 zero,half,one,two
640      parameter(zero=0.0d0)
641      parameter(half=0.5d0)
642      parameter(one=1.0d0)
643      parameter(two=2.0d0)
644c
645      integer n,na,nrot
646      real*8 a(na,na),d(na),v(na,na)
647      real*8 at,b,dma,q
648c
649      integer i,j,k,l
650      real*8 c,s,t,sum,temp
651c
652      do 1 i=1,n
653      do 2 j=1,n
654      v(i,j)=zero
655    2 continue
656      v(i,i)=one
657      d(i)=a(i,i)
658    1 continue
659c
660      nrot=0
661      do 3 l=1,nrmax
662      nrot=nrot+1
663      sum=zero
664      do 4 i=1,n-1
665      do 5 j=i+1,n
666      sum=sum+abs(a(i,j))
667    5 continue
668    4 continue
669      if(sum.eq.zero) then
670      do 6 i=1,n-1
671      do 7 j=i+1,n
672      if(d(i).gt.d(j)) then
673      temp=d(i)
674      d(i)=d(j)
675      d(j)=temp
676      do 8 k=1,n
677      temp=v(k,i)
678      v(k,i)=v(k,j)
679      v(k,j)=temp
680    8 continue
681      endif
682    7 continue
683    6 continue
684      return
685      endif
686      do 9 j=2,n
687      do 10 i=1,j-1
688      b=a(i,j)
689      if(abs(b).gt.zero) then
690      dma=d(j)-d(i)
691      if(abs(dma)+abs(b).le.abs(dma)) then
692      t=b/dma
693      else
694      q=half*dma/b
695      t=sign(one/(abs(q)+sqrt(one+q*q)),q)
696      endif
697      c=one/sqrt(t*t+one)
698      s=t*c
699      a(i,j)=zero
700      do 11 k=1,i-1
701      at=c*a(k,i)-s*a(k,j)
702      a(k,j)=s*a(k,i)+c*a(k,j)
703      a(k,i)=at
704   11 continue
705      do 12 k=i+1,j-1
706      at=c*a(i,k)-s*a(k,j)
707      a(k,j)=s*a(i,k)+c*a(k,j)
708      a(i,k)=at
709   12 continue
710      do 13 k=j+1,n
711      at=c*a(i,k)-s*a(j,k)
712      a(j,k)=s*a(i,k)+c*a(j,k)
713      a(i,k)=at
714   13 continue
715      do 14 k=1,n
716      at=c*v(k,i)-s*v(k,j)
717      v(k,j)=s*v(k,i)+c*v(k,j)
718      v(k,i)=at
719   14 continue
720      at=c*c*d(i)+s*s*d(j)-two*c*s*b
721      d(j)=s*s*d(i)+c*c*d(j)+two*c*s*b
722      d(i)=at
723      endif
724   10 continue
725    9 continue
726    3 continue
727c
728      call md_abort('md_jacobi: maximum iterations reached',0)
729c
730      return
731      end
732      subroutine swatch(today,now)
733c
734      implicit none
735c
736      character*10 today,now
737c
738#if defined(LINUX)
739      character*26 string
740#endif
741#if defined(IBM)
742      character*26 string
743#endif
744#if defined(SP1) || defined(SOLARIS)
745      character*26 string
746#endif
747#if defined(SGI)
748      character*9 string
749#endif
750c
751      today='00/00/00  '
752      now='00:00:00'
753c
754#if defined(IBM)
755      call fdate(string)
756      if(string(4:6).eq.'Jan') today(1:2)='01'
757      if(string(4:6).eq.'Feb') today(1:2)='02'
758      if(string(4:6).eq.'Mar') today(1:2)='03'
759      if(string(4:6).eq.'Apr') today(1:2)='04'
760      if(string(4:6).eq.'May') today(1:2)='05'
761      if(string(4:6).eq.'Jun') today(1:2)='06'
762      if(string(4:6).eq.'Jul') today(1:2)='07'
763      if(string(4:6).eq.'Aug') today(1:2)='08'
764      if(string(4:6).eq.'Sep') today(1:2)='09'
765      if(string(4:6).eq.'Oct') today(1:2)='10'
766      if(string(4:6).eq.'Nov') today(1:2)='11'
767      if(string(4:6).eq.'Dec') today(1:2)='12'
768      today(7:8)=string(8:9)
769      today(4:5)=string(1:2)
770      now=string(11:20)
771#endif
772#if defined(SP1) || defined(SOLARIS)
773      call util_date(string)
774      if(string(5:7).eq.'Jan') today(1:2)='01'
775      if(string(5:7).eq.'Feb') today(1:2)='02'
776      if(string(5:7).eq.'Mar') today(1:2)='03'
777      if(string(5:7).eq.'Apr') today(1:2)='04'
778      if(string(5:7).eq.'May') today(1:2)='05'
779      if(string(5:7).eq.'Jun') today(1:2)='06'
780      if(string(5:7).eq.'Jul') today(1:2)='07'
781      if(string(5:7).eq.'Aug') today(1:2)='08'
782      if(string(5:7).eq.'Sep') today(1:2)='09'
783      if(string(5:7).eq.'Oct') today(1:2)='10'
784      if(string(5:7).eq.'Nov') today(1:2)='11'
785      if(string(5:7).eq.'Dec') today(1:2)='12'
786      today(7:8)=string(23:24)
787      today(4:5)=string(9:10)
788      now=string(11:20)
789#endif
790#if defined(LINUX)
791      call util_date(string)
792      if(string(5:7).eq.'Jan') today(1:2)='01'
793      if(string(5:7).eq.'Feb') today(1:2)='02'
794      if(string(5:7).eq.'Mar') today(1:2)='03'
795      if(string(5:7).eq.'Apr') today(1:2)='04'
796      if(string(5:7).eq.'May') today(1:2)='05'
797      if(string(5:7).eq.'Jun') today(1:2)='06'
798      if(string(5:7).eq.'Jul') today(1:2)='07'
799      if(string(5:7).eq.'Aug') today(1:2)='08'
800      if(string(5:7).eq.'Sep') today(1:2)='09'
801      if(string(5:7).eq.'Oct') today(1:2)='10'
802      if(string(5:7).eq.'Nov') today(1:2)='11'
803      if(string(5:7).eq.'Dec') today(1:2)='12'
804      today(7:8)=string(23:24)
805      today(4:5)=string(9:10)
806      now=string(11:20)
807#endif
808#if defined(SGI)
809      call date(string)
810      if(string(4:6).eq.'Jan') today(1:2)='01'
811      if(string(4:6).eq.'Feb') today(1:2)='02'
812      if(string(4:6).eq.'Mar') today(1:2)='03'
813      if(string(4:6).eq.'Apr') today(1:2)='04'
814      if(string(4:6).eq.'May') today(1:2)='05'
815      if(string(4:6).eq.'Jun') today(1:2)='06'
816      if(string(4:6).eq.'Jul') today(1:2)='07'
817      if(string(4:6).eq.'Aug') today(1:2)='08'
818      if(string(4:6).eq.'Sep') today(1:2)='09'
819      if(string(4:6).eq.'Oct') today(1:2)='10'
820      if(string(4:6).eq.'Nov') today(1:2)='11'
821      if(string(4:6).eq.'Dec') today(1:2)='12'
822      today(7:8)=string(8:9)
823      today(4:5)=string(1:2)
824      call time(now(1:8))
825      now(9:10)='  '
826#endif
827      if(today(4:4).eq.' ') today(4:4)='0'
828      return
829      end
830      subroutine matinv(a,n,ndim)
831c
832      implicit none
833
834      integer maxdim
835      real*8 zero,small,one
836      parameter(maxdim=3)
837      parameter(zero=0.0d0)
838      parameter(small=1.0d-6)
839      parameter(one=1.0d0)
840c
841      integer n,ndim
842      real*8 a(ndim,ndim)
843      integer ia(2,maxdim),ib(maxdim),ic(maxdim)
844      real*8 d(maxdim)
845c
846      integer idim,i,j,k,l,m
847      real*8 b,e
848c
849      if(ndim.gt.maxdim) call md_abort('matinv dimension error',0)
850c
851      do 1 idim=1,n
852      ia(1,idim)=0
853      ia(2,idim)=0
854    1 continue
855c
856      do 9 idim=1,n
857      b=zero
858      do 3 l=1,n
859      do 4 m=1,n
860      if(ia(1,l).ne.1.and.ia(2,m).ne.1) then
861      e=dabs(a(l,m))
862      if(e.ge.b) then
863      i=l
864      k=m
865      endif
866    8 b=dmax1(b,e)
867      endif
868    4 continue
869    3 continue
870      ia(1,i)=1
871      ia(2,k)=1
872      ib(k)=i
873      ic(i)=k
874      b=a(i,k)
875c
876      if(dabs(b).lt.small) call md_abort('arg_matinv singular matrix',0)
877      a(i,k)=one/b
878      do 6 l=1,n
879      if(l.ne.k) a(i,l)=-a(i,l)/b
880  6   continue
881      do 5 l=1,n
882      do 7 m=1,n
883      if(l.ne.i.and.m.ne.k) a(l,m)=a(l,m)+a(l,k)*a(i,m)
884    7 continue
885  5   continue
886      do 11 l=1,n
887      if(l.ne.i) a(l,k)=a(l,k)/b
888  11  continue
889  9   continue
890c
891      do 15 l=1,n
892      do 13 j=1,n
893      k=ib(j)
894      d(j)=a(k,l)
895   13 continue
896      do 14 j=1,n
897      a(j,l)=d(j)
898   14 continue
899  15  continue
900c
901      do 16 l=1,n
902      do 17 j=1,n
903      k=ic(j)
904      d(j)=a(l,k)
905   17 continue
906      do 18 j=1,n
907      a(l,j)=d(j)
908   18 continue
909  16  continue
910c
911      return
912      end
913      logical function frequency(istep,nstep)
914c
915      implicit none
916c
917      integer istep,nstep
918c
919      if(nstep.le.0) then
920      frequency=.false.
921      else
922      frequency=mod(istep,nstep).eq.0
923      endif
924c
925      return
926      end
927      subroutine rolex(elaps,cputim)
928c
929      implicit none
930c
931      real*8 elaps,cputim
932c
933      real*8 util_wallsec,util_cpusec
934      external util_wallsec,util_cpusec
935c
936      elaps=util_wallsec()
937      cputim=util_cpusec()
938c
939      return
940      end
941      subroutine timer_init()
942c
943      implicit none
944c
945      call timer(0,0)
946c
947      return
948      end
949      subroutine timer_reset(itime)
950c
951      implicit none
952c
953      integer itime
954c
955      call timer(itime,-2)
956c
957      return
958      end
959      subroutine timer_start(itime)
960c
961      implicit none
962c
963      integer itime
964c
965      call timer(itime,0)
966c
967      return
968      end
969      subroutine timer_stop(itime)
970c
971      implicit none
972c
973      integer itime
974c
975      call timer(itime,1)
976c
977      return
978      end
979      subroutine timer(itime,iopt)
980c
981      implicit none
982c
983      integer itime,iopt
984c
985      real*8 elaps,cputim
986c
987      integer mtime
988      parameter(mtime=250)
989      integer ncall(250)
990      real*8 ttime(250,3),ctime(250,3)
991      common/tim/ncall,ttime,ctime
992c
993      integer i
994c
995      call rolex(elaps,cputim)
996c
997      if(itime.eq.0) then
998      do 1 i=1,mtime
999      ncall(i)=0
1000      ctime(i,1)=0.0d0
1001      ttime(i,1)=0.0d0
1002      ctime(i,2)=0.0d0
1003      ttime(i,2)=0.0d0
1004      ctime(i,3)=1.0d9
1005      ttime(i,3)=1.0d9
1006    1 continue
1007      elseif(itime.le.0.or.itime.gt.mtime) then
1008      call md_abort('Timer index out of range',0)
1009      elseif(iopt.eq.-2) then
1010      ncall(itime)=0
1011      ctime(itime,1)=0.0d0
1012      ttime(itime,1)=0.0d0
1013      ctime(itime,2)=0.0d0
1014      ttime(itime,2)=0.0d0
1015      ctime(itime,3)=1.0d9
1016      ttime(itime,3)=1.0d9
1017      elseif(iopt.eq.-1) then
1018      ncall(itime)=0
1019      ctime(itime,1)=-cputim
1020      ttime(itime,1)=-elaps
1021      ctime(itime,2)=-cputim
1022      ttime(itime,2)=-elaps
1023      ctime(itime,3)=1.0d9
1024      ttime(itime,3)=1.0d9
1025      elseif(iopt.eq.0) then
1026      ctime(itime,1)=ctime(itime,1)-cputim
1027      ttime(itime,1)=ttime(itime,1)-elaps
1028      ctime(itime,2)=-cputim
1029      ttime(itime,2)=-elaps
1030      elseif(iopt.eq.1) then
1031      ncall(itime)=ncall(itime)+1
1032      ctime(itime,1)=ctime(itime,1)+cputim
1033      ttime(itime,1)=ttime(itime,1)+elaps
1034      ctime(itime,2)=ctime(itime,2)+cputim
1035      ttime(itime,2)=ttime(itime,2)+elaps
1036      ctime(itime,3)=min(ctime(itime,2),ctime(itime,3))
1037      ttime(itime,3)=min(ttime(itime,2),ttime(itime,3))
1038      else
1039      call md_abort('Unimplemented timer option',0)
1040      endif
1041c
1042      return
1043      end
1044      real*8 function timer_cpu(itime)
1045c
1046      implicit none
1047c
1048      integer itime
1049c
1050      integer mtime
1051      parameter(mtime=250)
1052      integer ncall(250)
1053      real*8 ttime(250,3),ctime(250,3)
1054      common/tim/ncall,ttime,ctime
1055c
1056      if(itime.le.0.or.itime.gt.mtime)
1057     + call md_abort('Illegal timer index',0)
1058c
1059      if(ncall(itime).le.0) then
1060      timer_cpu=0.0d0
1061      else
1062      timer_cpu=ctime(itime,2)
1063      endif
1064c
1065      return
1066      end
1067      real*8 function timer_wall(itime)
1068c
1069      implicit none
1070c
1071      integer itime
1072c
1073      integer mtime
1074      parameter(mtime=250)
1075      integer ncall(250)
1076      real*8 ttime(250,3),ctime(250,3)
1077      common/tim/ncall,ttime,ctime
1078c
1079      if(itime.le.0.or.itime.gt.mtime)
1080     + call md_abort('Illegal timer index',0)
1081c
1082      if(ncall(itime).le.0) then
1083      timer_wall=0.0d0
1084      else
1085      timer_wall=ttime(itime,2)
1086      endif
1087c
1088      return
1089      end
1090      integer function timer_calls(itime)
1091c
1092      integer itime
1093c
1094      integer mtime
1095      parameter(mtime=250)
1096      integer ncall(250)
1097      real*8 ttime(250,3),ctime(250,3)
1098      common/tim/ncall,ttime,ctime
1099c
1100      if(itime.le.0.or.itime.gt.mtime)
1101     + call md_abort('Illegal timer index',0)
1102c
1103      timer_calls=ncall(itime)
1104c
1105      return
1106      end
1107      real*8 function timer_cpu_minimum(itime)
1108c
1109      implicit none
1110c
1111      integer itime
1112c
1113      integer mtime
1114      parameter(mtime=250)
1115      integer ncall(250)
1116      real*8 ttime(250,3),ctime(250,3)
1117      common/tim/ncall,ttime,ctime
1118c
1119      if(itime.le.0.or.itime.gt.mtime)
1120     + call md_abort('Illegal timer index',0)
1121c
1122      if(ncall(itime).le.0) then
1123      timer_cpu_minimum=0.0d0
1124      else
1125      timer_cpu_minimum=ctime(itime,3)
1126      endif
1127c
1128      return
1129      end
1130      real*8 function timer_wall_minimum(itime)
1131c
1132      implicit none
1133c
1134      integer itime
1135c
1136      integer mtime
1137      parameter(mtime=250)
1138      integer ncall(250)
1139      real*8 ttime(250,3),ctime(250,3)
1140      common/tim/ncall,ttime,ctime
1141c
1142      if(itime.le.0.or.itime.gt.mtime)
1143     + call md_abort('Illegal timer index',0)
1144c
1145      if(ncall(itime).le.0) then
1146      timer_wall_minimum=0.0d0
1147      else
1148      timer_wall_minimum=ttime(itime,3)
1149      endif
1150c
1151      return
1152      end
1153      real*8 function timer_cpu_average(itime)
1154c
1155      implicit none
1156c
1157      integer itime
1158c
1159      integer mtime
1160      parameter(mtime=250)
1161      integer ncall(250)
1162      real*8 ttime(250,3),ctime(250,3)
1163      common/tim/ncall,ttime,ctime
1164c
1165      if(itime.le.0.or.itime.gt.mtime)
1166     + call md_abort('Illegal timer index',0)
1167c
1168      if(ncall(itime).le.0) then
1169      timer_cpu_average=0.0d0
1170      else
1171      timer_cpu_average=ctime(itime,1)/dble(ncall(itime))
1172      endif
1173c
1174      return
1175      end
1176      real*8 function timer_wall_average(itime)
1177c
1178      implicit none
1179c
1180      integer itime
1181c
1182      integer mtime
1183      parameter(mtime=250)
1184      integer ncall(250)
1185      real*8 ttime(250,3),ctime(250,3)
1186      common/tim/ncall,ttime,ctime
1187c
1188      if(itime.le.0.or.itime.gt.mtime)
1189     + call md_abort('Illegal timer index',0)
1190c
1191      if(ncall(itime).le.0) then
1192      timer_wall_average=0.0d0
1193      else
1194      timer_wall_average=ttime(itime,1)/dble(ncall(itime))
1195      endif
1196c
1197      return
1198      end
1199      real*8 function timer_cpu_total(itime)
1200c
1201      implicit none
1202c
1203      integer itime
1204c
1205      integer mtime
1206      parameter(mtime=250)
1207      integer ncall(250)
1208      real*8 ttime(250,3),ctime(250,3)
1209      common/tim/ncall,ttime,ctime
1210c
1211      if(itime.le.0.or.itime.gt.mtime)
1212     + call md_abort('Illegal timer index',0)
1213c
1214      if(ncall(itime).le.0) then
1215      timer_cpu_total=0.0d0
1216      else
1217      timer_cpu_total=ctime(itime,1)
1218      endif
1219c
1220      return
1221      end
1222      real*8 function timer_wall_total(itime)
1223c
1224      implicit none
1225c
1226      integer itime
1227c
1228      integer mtime
1229      parameter(mtime=250)
1230      integer ncall(250)
1231      real*8 ttime(250,3),ctime(250,3)
1232      common/tim/ncall,ttime,ctime
1233c
1234      if(itime.le.0.or.itime.gt.mtime)
1235     + call md_abort('Illegal timer index',0)
1236c
1237      if(ncall(itime).le.0) then
1238      timer_wall_total=0.0d0
1239      else
1240      timer_wall_total=ttime(itime,1)
1241      endif
1242c
1243      return
1244      end
1245      subroutine error(lauto,lapprox,maxacf,data,ndata,
1246     + aver,drift,stderr,corerr,ratio)
1247c
1248      implicit none
1249c
1250      real*8 zero,one
1251      parameter(zero=0.0d0)
1252      parameter(one=1.0d0)
1253c
1254#include "mafdecls.fh"
1255c
1256      logical lauto,lapprox
1257      integer ndata,lenacf,maxacf
1258      real*8 data(ndata),aver,drift,stderr,corerr,ratio
1259c
1260      integer i,i_acf,l_acf
1261      real*8 dsum,ddsum,dtsum,tsum,ttsum,dstep
1262c
1263c      integer nacf,kapprx(15),iapprx,klarge
1264c      real*8 warg
1265c      real*8 data(ndata),acf(ndata),approx(15),cdac(15),weight
1266c
1267c      integer i,j,k,l,m,nacfa,nfunc
1268c      real*8 dsum,ddsum,dtsum,tsum,ttsum,dstep,dfsum,dvar
1269c      real*8 cdap(15,15),cdaq(15,15),cdad(15),rfact
1270c      real*8 xappm,xappi,xappj,wterm,wsum1,wsum2,cdawgt
1271c
1272c     arg_error iopt   = 0 : average
1273c                            drift
1274c                            standard error
1275c                      = 1 : average
1276c                            drift
1277c                            standard error
1278c                            autocorrelation function
1279c                            correlation error from actual acf
1280c                            sampling ration from actual acf
1281c                      = 2 : average
1282c                            drift
1283c                            standard error
1284c                            autocorrelation function
1285c                            correlation error from approximated acf
1286c                            sampling ration from approximated acf
1287c               data   : 1-dimensional array with data
1288c               ndata  : number of data
1289c               nacf   : length of autocorrelation function
1290c
1291c               aver   : average
1292c               stderr : standard error
1293c               drift  : drift
1294c               acf    : autocorrelation function
1295c               corerr : correlated error
1296c               ratio  : sampling ratio
1297c
1298      lenacf=maxacf
1299c
1300      dsum=zero
1301      ddsum=zero
1302      dtsum=zero
1303      tsum=zero
1304      ttsum=zero
1305      do 1 i=1,ndata
1306      dstep=dble(i)
1307      dsum=dsum+data(i)
1308      ddsum=ddsum+data(i)*data(i)
1309      dtsum=dtsum+dstep*data(i)
1310      tsum=tsum+dstep
1311      ttsum=ttsum+dstep*dstep
1312    1 continue
1313c
1314c     average, drift and standard error
1315c
1316      aver=dsum/dble(ndata)
1317      drift=(dble(ndata)*dtsum-dsum*tsum)/(dble(ndata)*ttsum-tsum*tsum)
1318      stderr=sqrt(abs((ddsum/dble(ndata)-dsum*dsum/dble(ndata*ndata)))/
1319     + (dble(ndata-1)))
1320c
1321      corerr=stderr
1322      ratio=one
1323c
1324      if(.not.lauto) return
1325c
1326      if(.not.ma_push_get(mt_dbl,lenacf,'acf',l_acf,i_acf))
1327     + call md_abort('Failed to allocate memory for acf',0)
1328c
1329      call auto_corr(data,ndata,aver,dbl_mb(i_acf),lenacf,ratio)
1330c
1331      if(.not.ma_pop_stack(l_acf))
1332     + call md_abort('Failed to deallocate memory for acf',0)
1333c
1334      corerr=ratio*stderr
1335c
1336      return
1337      end
1338      subroutine auto_corr(data,ndata,aver,acf,lacf,ratio)
1339c
1340      implicit none
1341c
1342      real*8 zero,half,one,two
1343      parameter(zero=0.0d0)
1344      parameter(half=5.0d-1)
1345      parameter(one=1.0d0)
1346      parameter(two=2.0d0)
1347c
1348      integer ndata,lacf,nacf
1349      real*8 data(ndata),acf(lacf),aver,ratio
1350c
1351      integer i,j
1352      real*8 dsum,dvar
1353c
1354      nacf=min(ndata,lacf)
1355c
1356      dsum=zero
1357      do 1 i=1,ndata
1358      dsum=dsum+(data(i)-aver)**2
1359    1 continue
1360      dvar=dble(ndata)/dsum
1361      do 2 i=1,nacf-1
1362      dsum=zero
1363      do 3 j=1,ndata-i
1364      dsum=dsum+(data(j)-aver)*(data(i+j)-aver)
1365    3 continue
1366      acf(i)=dvar*(dsum/dble(ndata-i))*(one-dble(i-1)/dble(ndata-2))
1367    2 continue
1368c
1369      ratio=half
1370      do 4 i=1,nacf-1
1371      ratio=ratio+(one-(dble(i)/dble(ndata))**2)*abs(acf(i))
1372    4 continue
1373      ratio=sqrt(two*abs(ratio))
1374      if(ratio.lt.one) ratio=one
1375c
1376      lacf=nacf
1377c
1378      return
1379      end
1380      subroutine acf_approx(acf,acfapp,lacf)
1381c
1382      implicit none
1383c
1384      integer lacf
1385      real*8 acf(lacf),acfapp(lacf)
1386c
1387c      integer kapprox(15)
1388c
1389c      rfact=two*sqrt(dble(klarge))/dble(nacfa-1)
1390c      wsum1=zero
1391c      wsum2=zero
1392c      warg=zero
1393c      do 5 i=1,nacfa-1
1394c      acf(i)=abs(acf(i))
1395c      if(acf(i).gt.zero) then
1396c      wterm=exp(weight*dble(nacfa-i)/dble(nacfa-1))
1397c      wsum1=wsum1+wterm
1398c      wsum2=wsum2+wterm*log(acf(i))/(rfact*dble(i))
1399c      endif
1400c    5 continue
1401c      if(abs(wsum1).gt.small) warg=wsum2/wsum1
1402c      do 55 i=1,nacfa-1
1403c      acf(i)=acf(i)-exp(warg*dble(i)*rfact)
1404c   55 continue
1405c      nfunc=iapprx
1406c      if(nfunc.le.0) nfunc=15
1407c      do 6 k=1,nfunc
1408c      do 7 l=1,nfunc
1409c      cdap(k,l)=zero
1410c      do 8 m=1,nacfa-1
1411c      xappm=dble(m)*rfact
1412c      cdawgt=exp(weight*dble(nacfa-m)/dble(nacfa-1))
1413c      cdap(k,l)=cdap(k,l)+cdawgt*exp((-xappm)*xappm)*
1414c     + approx(k)*approx(l)*(xappm**(kapprx(k)+kapprx(l)))
1415c    8 continue
1416c      cdaq(k,l)=cdap(k,l)
1417c    7 continue
1418c    6 continue
1419c      do 9 i=1,nfunc
1420c      cdad(i)=zero
1421c      do 10 j=1,nacfa-1
1422c      xappj=dble(j)*rfact
1423c      cdawgt=exp(weight*dble(nacfa-j)/dble(nacfa-1))
1424c      cdad(i)=cdad(i)+cdawgt*exp((-half)*xappj*xappj)*
1425c     + acf(j)*approx(i)*(xappj**kapprx(i))
1426c   10 continue
1427c    9 continue
1428c      do 11 k=1,nfunc
1429c      do 12 i=k,nfunc
1430c      cdaq(i,k)=cdap(i,k)
1431c      do 13 l=1,k-1
1432c      cdaq(i,k)=cdaq(i,k)-cdaq(i,l)*cdaq(l,k)
1433c   13 continue
1434c   12 continue
1435c      do 14 i=k+1,nfunc
1436c      cdaq(k,i)=cdap(k,i)
1437c      do 15 l=1,k-1
1438c      cdaq(k,i)=cdaq(k,i)-cdaq(k,l)*cdaq(l,i)
1439c   15 continue
1440c      cdaq(k,i)=cdaq(k,i)/cdaq(k,k)
1441c   14 continue
1442c   11 continue
1443c      do 16 j=1,nfunc
1444c      cdac(j)=cdad(j)
1445c      do 17 i=1,j-1
1446c      cdac(j)=cdac(j)-cdaq(j,i)*cdac(i)
1447c   17 continue
1448c      cdac(j)=cdac(j)/cdaq(j,j)
1449c   16 continue
1450c      do 18 k=1,nfunc
1451c      j=nfunc+1-k
1452c      do 19 i=j+1,nfunc
1453c      cdac(j)=cdac(j)-cdaq(j,i)*cdac(i)
1454c   19 continue
1455c   18 continue
1456c      do 20 i=1,nacfa-1
1457c      xappi=dble(i)*rfact
1458c      acf(i)=exp(warg*xappi)
1459c      do 21 j=1,nfunc
1460c      acf(i)=acf(i)+exp((-half)*xappi*xappi)*
1461c     + approx(j)*cdac(j)*(xappi**kapprx(j))
1462c   21 continue
1463c   20 continue
1464c
1465c     autocorrelation function upto lag min(nacf,ndata)
1466c
1467c      if(iopt.gt.0) then
1468c      warg=zero
1469c      nacfa=nacf
1470c      if(nacfa.gt.ndata) nacfa=ndata
1471c      dfsum=zero
1472c      do 2 i=1,ndata
1473c      dfsum=dfsum+(data(i)-aver)**2
1474c    2 continue
1475c      dvar=dble(ndata)/dfsum
1476c      do 3 i=1,nacfa-1
1477c      dfsum=zero
1478c      do 4 j=1,ndata-i
1479c      dfsum=dfsum+(data(j)-aver)*(data(i+j)-aver)
1480c    4 continue
1481c      acf(i)=dvar*(dfsum/dble(ndata-i))*(one-dble(i-1)/dble(ndata-2))
1482c    3 continue
1483c      endif
1484cc
1485cc     approximate acf here
1486cc
1487c      if(iopt.gt.1) then
1488c      rfact=two*sqrt(dble(klarge))/dble(nacfa-1)
1489c      wsum1=zero
1490c      wsum2=zero
1491c      warg=zero
1492c      do 5 i=1,nacfa-1
1493c      acf(i)=abs(acf(i))
1494c      if(acf(i).gt.zero) then
1495c      wterm=exp(weight*dble(nacfa-i)/dble(nacfa-1))
1496c      wsum1=wsum1+wterm
1497c      wsum2=wsum2+wterm*log(acf(i))/(rfact*dble(i))
1498c      endif
1499c    5 continue
1500c      if(abs(wsum1).gt.small) warg=wsum2/wsum1
1501c      do 55 i=1,nacfa-1
1502c      acf(i)=acf(i)-exp(warg*dble(i)*rfact)
1503c   55 continue
1504c      nfunc=iapprx
1505c      if(nfunc.le.0) nfunc=15
1506c      do 6 k=1,nfunc
1507c      do 7 l=1,nfunc
1508c      cdap(k,l)=zero
1509c      do 8 m=1,nacfa-1
1510c      xappm=dble(m)*rfact
1511c      cdawgt=exp(weight*dble(nacfa-m)/dble(nacfa-1))
1512c      cdap(k,l)=cdap(k,l)+cdawgt*exp((-xappm)*xappm)*
1513c     + approx(k)*approx(l)*(xappm**(kapprx(k)+kapprx(l)))
1514c    8 continue
1515c      cdaq(k,l)=cdap(k,l)
1516c    7 continue
1517c    6 continue
1518c      do 9 i=1,nfunc
1519c      cdad(i)=zero
1520c      do 10 j=1,nacfa-1
1521c      xappj=dble(j)*rfact
1522c      cdawgt=exp(weight*dble(nacfa-j)/dble(nacfa-1))
1523c      cdad(i)=cdad(i)+cdawgt*exp((-half)*xappj*xappj)*
1524c     + acf(j)*approx(i)*(xappj**kapprx(i))
1525c   10 continue
1526c    9 continue
1527c      do 11 k=1,nfunc
1528c      do 12 i=k,nfunc
1529c      cdaq(i,k)=cdap(i,k)
1530c      do 13 l=1,k-1
1531c      cdaq(i,k)=cdaq(i,k)-cdaq(i,l)*cdaq(l,k)
1532c   13 continue
1533c   12 continue
1534c      do 14 i=k+1,nfunc
1535c      cdaq(k,i)=cdap(k,i)
1536c      do 15 l=1,k-1
1537c      cdaq(k,i)=cdaq(k,i)-cdaq(k,l)*cdaq(l,i)
1538c   15 continue
1539c      cdaq(k,i)=cdaq(k,i)/cdaq(k,k)
1540c   14 continue
1541c   11 continue
1542c      do 16 j=1,nfunc
1543c      cdac(j)=cdad(j)
1544c      do 17 i=1,j-1
1545c      cdac(j)=cdac(j)-cdaq(j,i)*cdac(i)
1546c   17 continue
1547c      cdac(j)=cdac(j)/cdaq(j,j)
1548c   16 continue
1549c      do 18 k=1,nfunc
1550c      j=nfunc+1-k
1551c      do 19 i=j+1,nfunc
1552c      cdac(j)=cdac(j)-cdaq(j,i)*cdac(i)
1553c   19 continue
1554c   18 continue
1555c      do 20 i=1,nacfa-1
1556c      xappi=dble(i)*rfact
1557c      acf(i)=exp(warg*xappi)
1558c      do 21 j=1,nfunc
1559c      acf(i)=acf(i)+exp((-half)*xappi*xappi)*
1560c     + approx(j)*cdac(j)*(xappi**kapprx(j))
1561c   21 continue
1562c   20 continue
1563c      endif
1564cc
1565cc     sampling ratio
1566cc
1567c      ratio=one
1568c      if(iopt.gt.0) then
1569c      ratio=half
1570c      do 22 i=1,nacfa-1
1571c      ratio=ratio+(one-(dble(i)/dble(ndata))**2)*abs(acf(i))
1572c   22 continue
1573c      ratio=sqrt(two*abs(ratio))
1574c      if(ratio.lt.one) ratio=one
1575c      endif
1576cc
1577cc
1578c      corerr=ratio*stderr
1579cc
1580      return
1581      end
1582      logical function md_zmat(x1,x2,x3,x4,d,a,t)
1583c
1584      implicit none
1585c
1586      real*8 x1(3),x2(3),x3(3),x4(3),d,a,t
1587c
1588      real*8 p(3),q(3),s(3),v(3),w(3)
1589      real*8 r
1590      integer i
1591c
1592      do 1 i=1,3
1593      p(i)=x3(i)-x2(i)
1594      q(i)=x4(i)-x2(i)
1595    1 continue
1596      v(1)=p(2)*q(3)-p(3)*q(2)+x2(1)
1597      v(2)=p(3)*q(1)-p(1)*q(3)+x2(2)
1598      v(3)=p(1)*q(2)-p(2)*q(1)+x2(3)
1599c
1600      r=sqrt(p(1)*p(1)+p(2)*p(2)+p(3)*p(3))
1601      do 2 i=1,3
1602      s(i)=d*p(i)/r+x2(i)
1603    2 continue
1604c
1605      call rotate(x2,v,a,s,w)
1606      call rotate(x3,x2,t,w,s)
1607c
1608      do 3 i=1,3
1609      x1(i)=s(i)
1610    3 continue
1611c
1612      md_zmat=.true.
1613      return
1614      end
1615      subroutine md_abort(string, icode)
1616      implicit none
1617#include "global.fh"
1618#include "stdio.fh"
1619      character*(*) string
1620      character*255 card
1621      integer icode
1622      if(ga_nodeid().eq.0) then
1623      write(luout,1000) 0,string,icode
1624 1000 format(/,1x,10('*'),/,' * ',i3,': ',a,i5,/,1x,10('*'))
1625      card=' '
1626      else
1627      write(card,1001) ga_nodeid(),string,icode
1628 1001 format(' * ',i3,': ',a,i5,' * ')
1629      endif
1630      call ga_error(card,icode)
1631      return
1632      end
1633