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