1C $Id$
2      Subroutine drdy_react(jop)
3C
4C Calculate potential and hessian for reactant and product geometries
5C  and write out to file 30
6C
7      Implicit none
8      Integer jop
9      Integer i,iloop,iop,j,nbar
10      Double Precision swf
11C
12#include "drdyP.fh"
13C
14C  Geometries for reactants or products read in subroutine input
15C
16      swf = sreact
17      Do iloop = 1,2
18         iop = 2*iloop - 1
19         if (iop.gt.jop) then
20            do i = 1,n3
21               x(i) = xr(i,iop)
22            end do
23            if (iop.eq.1) then
24               write (fu6,600)
25               write (fu6,610)
26            else
27               write (fu6,620)
28               write (fu6,610)
29            endif
30            write (fu6,640) (j+1,(x(3*j+i),i=1,3),j=0,natom-1)
31C  Get potential and first and second derivatives
32            nbar = (n3*(n3+1))/2
33            Call drdy_pot2 (x,v,dx,f,hess,scr1,amass,natom,n3,n3tm,nbar)
34            if (iloop.eq.1) then
35               vr = v
36            else
37               vp = v
38            endif
39            if(lgs(38).gt.0) call drdy_wfiles(swf,0)
40            if (lgs(39).eq.0) then
41               write (fu6,650) v,v*ckcal
42            else
43               call drdy_potsp(x,vspc,scr1,amass,natom,n3)
44               if (iloop.eq.1) then
45                  vrsp = vspc
46               else
47                  vpsp = vspc
48               endif
49               write (fu6,651) v,v*ckcal,vspc,vspc*ckcal
50               if(lgs(38).eq.4) call drdy_wfiles(swf,1)
51            endif
52            write (fu6,640) (j+1,(dx(3*j+i),i=1,3),j=0,natom-1)
53C  put correct parts of hessian matrix into fsv and call rphwrt
54            Call drdy_react2(iop)
55         endif
56         if ((iop.eq.1 .and. lgs(6).le.2) .or.
57     *    (iop.eq.3 .and. (lgs(6).eq.1.or.lgs(6).eq.3))) then
58            iop = iop+1
59            if (iop.gt.jop) then
60               v = 0.0
61C  put correct parts of hessian matrix into fsv and call rphwrt
62               call drdy_react2(iop)
63            endif
64         end if
65         swf = sprod
66      enddo
67      RETURN
68C
69600   Format(//1x,10(1h*),' Reactants',/)
70610   Format(1x,'Geometry in mass-scaled cartesians (bohrs):',
71     *   //,17x,1hx,15x,1hy,15x,1hz)
72620   Format(//1x,10(1h*),' Products',/)
73640   Format(1x,i3,4x,1p,3e16.6)
74650   Format(/,' Potential energy=',T20,1pe18.10,' Hartree=',e13.5,
75     *   ' kcal/mol',//,' Derivatives in mass-scaled cartesians',
76     *   ' (au):')
77651   Format(/,' Potential energy=',T20,1pe18.10,' Hartree=',e13.5,
78     *   ' kcal/mol',/,5x,'Single point',T20,1pe18.10,' Hartree=',
79     *   e13.5,
80     *   ' kcal/mol',//,' Derivatives in mass-scaled cartesians (au):')
81      END
82