1*----------------------------------------------------------------------*
2      subroutine ucc_grad(ccvec1,ccvec2,civec1,civec2,
3     &                    n_cc_amp,
4     &                    lu_amp,lu_exrf,luhexrf)
5*----------------------------------------------------------------------*
6*     calculate the gradient of the Unitary CC Energy functional using
7*     the Wilson formula and Gauss-Legendre quadrature
8*
9*     lu_amp   contains the amplitudes
10*     lu_exrf  contains exp(sig)|R>
11*     luhexrf  contains H exp(sig)|R>
12*
13*     result is returned on ccvec1
14*
15*     andreas, adapted from gtbce_gradE, August 2004
16*
17*----------------------------------------------------------------------*
18
19      implicit none
20
21      integer, parameter ::
22     &     ntest = 00, npnts = 5, mx_term = -150, iopsym = -1
23      real(8), parameter ::
24     &     thresh = 1d-20
25
26      integer, intent(in) ::
27     &     n_cc_amp,
28     &     lu_amp, lu_exrf, luhexrf
29      real(8), intent(inout) ::
30     &     ccvec1(n_cc_amp), ccvec2(n_cc_amp),
31     &     civec1(*), civec2(*)
32
33      character ::
34     &     cctype*8
35      integer ::
36     &     lblk, mxpnts, ipnt,
37     &     lusc1,lusc2,lusc3,lusc4,lusc5,lusc6,lusc7
38      real(8) ::
39     &     cpu0, cpu, wall0, wall,
40     &     xnrm, etest, etest2, xnrm2, dltalp
41
42      real(8) ::
43     &     alp(npnts), wght(npnts)
44
45      integer, external ::
46     &     iopen_nus
47
48      real(8), external ::
49     &     inprod, inprdd
50
51      call atim(cpu0,wall0)
52
53      lblk = -1
54      if (ntest.ge.5) then
55        write (6,*) '====================='
56        write (6,*) ' This is ucc_grad'
57        write (6,*) '====================='
58        write (6,*)
59        write (6,*) 'on entry: '
60        write (6,*) ' lu_amp:  ', lu_amp
61        write (6,*) ' lu_exrf: ', lu_exrf
62        write (6,*) ' luhexrf: ', luhexrf
63      end if
64
65      if (ntest.ge.1000) then
66        write(6,*) 'on entry:'
67        write(6,*) 'e^G|0> on LU_EXRF'
68        call wrtvcd(civec1,lu_exrf,1,lblk)
69        write(6,*) 'H e^G|0> on LUHEXRF'
70        call wrtvcd(civec1,luhexrf,1,lblk)
71      end if
72
73      ! for I/O
74      lblk = -1
75      ! for expt_ref
76      cctype='GEN_CC'
77
78      lusc1 = iopen_nus('UCCGRDSCR1')
79      lusc2 = iopen_nus('UCCGRDSCR2')
80      lusc3 = iopen_nus('UCCGRDSCR3')
81      lusc4 = iopen_nus('UCCGRDSCR4')
82      lusc5 = iopen_nus('UCCGRDSCR5')
83      lusc6 = iopen_nus('UCCGRDSCR6')
84      lusc7 = iopen_nus('UCCGRDSCR7')
85
86*--------------------------------------------------------------------*
87* set up points and weights
88*--------------------------------------------------------------------*
89      call gl_weights(0d0,1d0,npnts,alp,wght)
90
91      mxpnts=npnts
92      ! if G == 0 ...
93      xnrm2 = inprod(ccvec1,ccvec1,n_cc_amp)
94      ! ... things are trivial and we evaluate the formula only once
95      if (xnrm2.lt.thresh) then
96        mxpnts=1
97        wght(1)=1d0
98        alp(1)=0d0
99        if (ntest.ge.5) then
100          write(6,*) 'Detected zero amplitudes: ',
101     &               'only case alpha = 0 will be processed'
102        end if
103      end if
104
105*--------------------------------------------------------------------*
106* H|0tilde>
107*  init lusc1
108*--------------------------------------------------------------------*
109      call copvcd(luhexrf,lusc1,civec1,1,lblk)
110
111**-------------------------------------------------------------------*
112* loop over quadrature points
113**-------------------------------------------------------------------*
114      do ipnt = 1, mxpnts
115        if (ntest.ge.5) then
116          write(6,*) 'info for quadrature point: ', ipnt,'/',npnts
117          write(6,*) 'point, weight: ', alp(ipnt), wght(ipnt)
118        end if
119
120        if (ipnt.gt.1.and.(alp(ipnt).le.alp(ipnt-1))) then
121          write(6,*) 'quadrature point should be in ascending order!'
122          stop 'ucc_grad > quadrature '
123        end if
124
125        ! read amplitudes
126        call vec_from_disc(ccvec1,n_cc_amp,1,lblk,lu_amp)
127
128        if (ipnt.eq.1) then
129          dltalp = alp(1)
130        else
131          dltalp = alp(ipnt)-alp(ipnt-1)
132          call copvcd(lusc2,lusc1,civec1,1,lblk)
133        end if
134*--------------------------------------------------------------------*
135* |a_i> = exp(a_i G^+) [H exp(G)|0>]
136*       = exp((a_i-a_{i-1})G^+)[exp(a_{i-1}G^+) H exp(G)|0>]
137*       = exp(-(a_i-a_{i-1})G) [exp(-a_{i-1}G) H exp(G)|0>]
138*  result on lusc2
139*--------------------------------------------------------------------*
140        if (ntest.ge.5) then
141          write(6,*)
142     &         'constructing |a_i> = exp(a_i G^+) [H exp(G)|0>]'
143        end if
144
145        if (abs(dltalp).lt.1d-20) then
146          call copvcd(lusc1,lusc2,civec1,1,lblk)
147        else
148          ! scale it (minus sign: G^+ = -G)
149          call scalve(ccvec1,-dltalp,n_cc_amp)
150          call expt_ref2(lusc1,lusc2,lusc4,lusc5,lusc6,
151     &         thresh,mx_term, ccvec1, ccvec2, civec1, civec2,
152     &         n_cc_amp,cctype, iopsym)
153          if (ntest.ge.5) then
154            xnrm = sqrt(inprod(ccvec1,ccvec1,n_cc_amp))
155            write(6,*) '|dlta G^+|, dlta = ',xnrm, dltalp,
156     &                 'for alp(i) = ', alp(ipnt)
157          end if
158
159        end if
160
161*--------------------------------------------------------------------*
162* |b_i> = exp(-a_i G)exp(G)|0> =
163*       = exp(-(a_i-a_{i-1})G) [exp(-a_{i-1}G)exp(G)|0>]
164*  result on lusc3
165*--------------------------------------------------------------------*
166        if (ipnt.eq.1) then
167          call copvcd(lu_exrf,lusc1,civec1,1,lblk)
168        else
169          call copvcd(lusc3,lusc1,civec1,1,lblk)
170        end if
171
172        if (ntest.ge.5) then
173          write(6,*) 'constructing |b_i> = exp(-a_i G) exp(G)|0>]'
174        end if
175
176        if (abs(dltalp).lt.1d-20) then
177          call copvcd(lusc1,lusc3,civec1,1,lblk)
178        else
179          ! scaled vector is already on ccvec1
180          call expt_ref2(lusc1,lusc3,lusc4,lusc5,lusc6,
181     &         thresh,mx_term, ccvec1, ccvec2, civec1, civec2,
182     &         n_cc_amp,cctype, iopsym)
183
184          if (ntest.ge.5) then
185            xnrm = sqrt(inprod(ccvec2,ccvec2,n_cc_amp))
186            etest = inprdd(civec1,civec2,lusc3,lusc3,1,lblk)
187            etest2= inprdd(civec1,civec2,lusc2,lusc3,1,lblk)
188            write(6,*) '|dltaG|, dlta = ',xnrm, dltalp
189            write(6,*) '<b_i|b_i> = ', etest,
190     &           'for alp(i) = ', alp(ipnt)
191            write(6,*) '<a_i|b_i>  = ', etest2,
192     &           'for alp(i) = ', alp(ipnt)
193          end if
194        end if
195
196*--------------------------------------------------------------------*
197* dE_u +=   w_i (<a_i|gamma_u-gamma^+_u|b_i>)
198*         = w_i (<a_i|gamma_u|b_i> - <b_i|gamma_u|a_i>)
199*--------------------------------------------------------------------*
200
201        if (ntest.ge.1000) then
202          write(6,*) 'Before calling sigden_cc:'
203          write(6,*) '|a_i> on lusc2:'
204          call wrtvcd(civec1,lusc2,1,lblk)
205          write(6,*) '|b_i> on lusc3:'
206          call wrtvcd(civec1,lusc3,1,lblk)
207        end if
208
209        call den_gcc_s(civec1,civec2,lusc2,lusc3,ccvec1,ccvec2,-1)
210
211        if (ipnt.gt.1) then
212          call vec_from_disc(ccvec2,n_cc_amp,1,lblk,lusc7)
213          call vecsum(ccvec1,ccvec1,ccvec2,2d0*wght(ipnt),1d0,n_cc_amp)
214        else
215          call scalve(ccvec1,2d0*wght(ipnt),n_cc_amp)
216        end if
217        if (ipnt.lt.mxpnts)
218     &       call vec_to_disc(ccvec1,n_cc_amp,1,lblk,lusc7)
219
220      end do ! loop over quadrature points
221
222
223      if (ntest.ge.5) then
224        xnrm = sqrt(inprod(ccvec1,ccvec1,n_cc_amp))
225        write(6,*) ' ucc_grad > '
226        write(6,*) '     n_cc_amp,norm of grad: ',n_cc_amp,xnrm
227      end if
228      if (ntest.ge.100) then
229        call wrt_cc_vec2(ccvec1,6,'GEN_CC')
230      end if
231
232c prelim: divide by 2d0 again
233      call scalve(ccvec1,0.5d0,n_cc_amp)
234
235      call relunit(lusc1,'delete')
236      call relunit(lusc2,'delete')
237      call relunit(lusc3,'delete')
238      call relunit(lusc4,'delete')
239      call relunit(lusc5,'delete')
240      call relunit(lusc6,'delete')
241      call relunit(lusc7,'delete')
242
243      call atim(cpu,wall)
244      call prtim(6,'time in ucc_grad',cpu-cpu0,wall-wall0)
245
246      return
247      end
248c $Id$
249