1* 2* $Id$ 3* 4 subroutine mcscf_hessv_2e_ao( geom, basis, nbf, nclosed, nact, 5 $ tol2e, oskel, dm1, g_movecs, 6 $ g_x, g_tmp, g_ax ) 7#include "global.fh" 8#include "errquit.fh" 9#include "mafdecls.fh" 10#include "mcscfprof.fh" 11c 12c 13 integer geom, basis ! [input] Handles 14 integer nbf ! [input] Basis functions 15 integer nclosed ! [input] Closed shells 16 integer nact ! [input] Active shells 17 double precision tol2e ! [input] Integral tolerance 18 logical oskel ! [input] Symmetry selection 19 double precision dm1(nact,nact) ! [input] 1PDM 20 integer g_movecs ! [input] MO coefficients 21 integer g_tmp ! [input] Temporary (nbf * nbf) 22 integer g_x ! [input] Argument parameter matrix 23 integer g_ax ! [output] Hessian product (in matrix rep) 24c 25c 26c 27 integer g_f1, g_f2, g_f3, g_f4, g_f5, g_f6 ! Atom-blocked Fock matrices 28 integer g_d1, g_d2, g_d3, g_d4, g_d5, g_d6 ! Atom-blocked densities 29 integer nvir, vlen, voff, aoff, aend, nn 30 integer g_tmp1, g_tmp2 31c 32c 33 integer nsets 34 parameter(nsets=6) 35 integer iv_dens(nsets), iv_fock(nsets) 36 double precision jfac(nsets),kfac(nsets) 37c 38c 39 integer ga_create_atom_blocked 40 external ga_create_atom_blocked 41c 42 data jfac/6*1.d0/ 43 data kfac/6*-0.5d0/ 44c 45c 46 if (omcscfprof) call pstat_on(ps_hv2eao) 47 nvir = nbf - nclosed - nact 48 vlen = (nclosed+nact)*nvir + nclosed*nact 49 voff = nclosed + nact + 1 50 aoff = nclosed + 1 51 aend = nclosed + nact 52c 53c 54 if (.not.ga_duplicate(g_tmp,g_tmp1,'temp1')) 55 $ call errquit('mcscf_hessv_2e_ao: cannot duplicate',0, 56 & GA_ERR) 57 if (.not.ga_duplicate(g_tmp,g_tmp2,'temp2')) 58 $ call errquit('mcscf_hessv_2e_ao: cannot duplicate',0, 59 & GA_ERR) 60 call ga_zero(g_tmp1) 61 call ga_zero(g_tmp2) 62 call ga_zero(g_tmp) 63 if (ga_nodeid().eq.0) 64 $ call ga_put(g_tmp,aoff,aend,aoff,aend,dm1,nact) 65 call ga_sync() 66 g_d1 = ga_create_atom_blocked( geom, basis, 'Density 1') 67 g_d2 = ga_create_atom_blocked( geom, basis, 'Density 2') 68 g_d3 = ga_create_atom_blocked( geom, basis, 'Density 3') 69 g_d4 = ga_create_atom_blocked( geom, basis, 'Density 4') 70 g_d5 = ga_create_atom_blocked( geom, basis, 'Density 5') 71 g_d6 = ga_create_atom_blocked( geom, basis, 'Density 6') 72 g_f1 = ga_create_atom_blocked( geom, basis, 'Fock 1') 73 g_f2 = ga_create_atom_blocked( geom, basis, 'Fock 2') 74 g_f3 = ga_create_atom_blocked( geom, basis, 'Fock 3') 75 g_f4 = ga_create_atom_blocked( geom, basis, 'Fock 4') 76 g_f5 = ga_create_atom_blocked( geom, basis, 'Fock 5') 77 g_f6 = ga_create_atom_blocked( geom, basis, 'Fock 6') 78c 79c Inactive-Virtual "density" 80c iv iv t 81c D = C.X .C 82c 83 call ga_matmul_patch( 'n', 't', 1.d0, 0.d0, 84 $ g_x, voff, nbf, 1, nclosed, 85 $ g_movecs, 1, nclosed, 1, nbf, 86 $ g_tmp1, voff, nbf, 1, nbf ) 87 call ga_matmul_patch( 'n', 'n', 1.d0, 0.d0, 88 $ g_movecs, 1, nbf, voff, nbf, 89 $ g_tmp1, voff, nbf, 1, nbf, 90 $ g_tmp2, 1, nbf, 1, nbf ) 91 call ga_symmetrize(g_tmp2) 92 call ga_copy( g_tmp2, g_d1 ) 93 call ga_copy( g_d1, g_d2 ) 94 call ga_copy( g_d1, g_d3 ) 95c 96c Inactive-Virtual, Active-Virtual "density" 97c ,av av t 98c D = C.X .d.C 99c 100 call ga_matmul_patch( 'n', 'n', 1.d0, 0.d0, 101 $ g_x, voff, nbf, aoff, aend, 102 $ g_tmp, aoff, aend, aoff, aend, 103 $ g_tmp1, voff, nbf, aoff, aend ) 104 call ga_matmul_patch( 'n', 't', 1.d0, 0.d0, 105 $ g_tmp1, voff, nbf, aoff, aend, 106 $ g_movecs, aoff, aend, 1, nbf, 107 $ g_tmp2, voff, nbf, 1, nbf ) 108 call ga_matmul_patch( 'n', 'n', 1.d0, 0.d0, 109 $ g_movecs, 1, nbf, voff, nbf, 110 $ g_tmp2, voff, nbf, 1, nbf, 111 $ g_tmp1, 1, nbf, 1, nbf ) 112 call ga_symmetrize(g_tmp1) 113 call ga_copy( g_tmp1, g_d4 ) 114 call ga_dadd( 1.d0, g_d4, 2.d0, g_d1, g_d1 ) 115c 116c Inactive-Virtual, Inactive-Active "density" 117c ,ia t ia 118c D = C.(2-d)X .C 119c 120 call ga_zero(g_tmp2) 121 call diagfill_patch(g_tmp2, 2.d0, aoff, aend ) 122 call ga_dadd(1.d0, g_tmp2, -1.d0, g_tmp, g_tmp2 ) 123 call ga_zero(g_tmp1) 124 call ga_matmul_patch('n', 'n', 1.d0, 0.d0, 125 $ g_tmp2, aoff, aend, aoff, aend, 126 $ g_x, aoff, aend, 1, nclosed, 127 $ g_tmp1, aoff, aend, 1, nclosed ) 128 call ga_matmul_patch('n', 't', 1.d0, 0.d0, 129 $ g_tmp1, aoff, aend, 1, nclosed, 130 $ g_movecs, 1, nclosed, 1, nbf, 131 $ g_tmp2, aoff, aend, 1, nbf ) 132 call ga_matmul_patch('n', 'n', 1.d0, 0.d0, 133 $ g_movecs, 1, nbf, aoff, aend, 134 $ g_tmp2, aoff, aend, 1, nbf, 135 $ g_tmp1, 1, nbf, 1, nbf ) 136 call ga_symmetrize(g_tmp1) 137 call ga_copy( g_tmp1, g_d5 ) 138 call ga_dadd( 1.d0, g_d5, 1.d0, g_d1, g_d1 ) 139c 140c Inactive-Active density 141c ia t ia 142c D = C.X .C 143c 144 call ga_zero(g_tmp1) 145 call ga_zero(g_tmp2) 146 call ga_matmul_patch('n', 't', 1.d0, 0.d0, 147 $ g_x, aoff, aend, 1, nclosed, 148 $ g_movecs, 1, nclosed, 1, nbf, 149 $ g_tmp1, aoff, aend, 1, nbf ) 150 call ga_matmul_patch('n', 'n', 1.d0, 0.d0, 151 $ g_movecs, 1, nbf, aoff, aend, 152 $ g_tmp1, aoff, aend, 1, nbf, 153 $ g_tmp2, 1, nbf, 1, nbf ) 154 call ga_symmetrize(g_tmp2) 155 call ga_copy( g_tmp2, g_d5 ) 156 call ga_dadd( 1.d0, g_d5, 1.d0, g_d2, g_d2 ) 157c 158c 159c CA, CA density 160c ,,ia t ia 161c D = -C.(1-d).X .C 162c 163 call ga_zero(g_tmp1) 164 call diagfill_patch(g_tmp1, 1.d0, aoff, aend ) 165 call ga_dadd( 1.d0, g_tmp1, -1.d0, g_tmp, g_tmp1) 166 call ga_matmul_patch( 'n', 'n', 1.d0, 0.d0, 167 $ g_tmp1, aoff, aend, aoff, aend, 168 $ g_x, aoff, aend, 1, nclosed, 169 $ g_tmp2, aoff, aend, 1, nclosed ) 170 call ga_matmul_patch( 'n', 't', 1.d0, 0.d0, 171 $ g_tmp2, aoff, aend, 1, nclosed, 172 $ g_movecs, 1, nclosed, 1, nbf, 173 $ g_tmp1, aoff, aend, 1, nbf ) 174 call ga_matmul_patch( 'n', 'n', 1.d0, 0.d0, 175 $ g_movecs, 1, nbf, aoff, aend, 176 $ g_tmp1, aoff, aend, 1, nbf, 177 $ g_tmp2, 1, nbf, 1, nbf ) 178 call ga_symmetrize(g_tmp2) 179 call ga_copy(g_tmp2, g_d6) 180c 181c 182c 183c Density summary 184c 185c iv ,av ,ia ,av 186c d1 : 2D + D + D d4 : D 187c 188c iv ia ia 189c d2 : D + D d5 : D 190c 191c iv ,,ia 192c d3 : D d6 : D 193c 194c 195 iv_dens(1) = g_d1 196 iv_dens(2) = g_d2 197 iv_dens(3) = g_d3 198 iv_dens(4) = g_d4 199 iv_dens(5) = g_d5 200 iv_dens(6) = g_d6 201c 202c Fock build 203c 204 nn = 6 205 call ga_zero(g_f1) 206 call ga_zero(g_f2) 207 call ga_zero(g_f3) 208 call ga_zero(g_f4) 209 call ga_zero(g_f5) 210 call ga_zero(g_f6) 211 iv_fock(1) = g_f1 212 iv_fock(2) = g_f2 213 iv_fock(3) = g_f3 214 iv_fock(4) = g_f4 215 iv_fock(5) = g_f5 216 iv_fock(6) = g_f6 217 call fock_2e( geom, basis, nn, jfac, kfac, tol2e, oskel, 218 $ iv_dens, iv_fock, .false. ) 219c 220c Symmetrize AO Fock matrices 221c 222 if (oskel) then 223 do i=1,nn 224 call sym_symmetrize(geom, basis, .false., iv_fock(i)) 225 enddo 226 endif 227c 228c Inactive-Virtual contribution 229c 230c 231 call ga_zero(g_tmp2) 232 call ga_matmul_patch( 't', 'n', 1.d0, 0.d0, 233 $ g_movecs, voff, nbf, 1, nbf, 234 $ g_f1, 1, nbf, 1, nbf, 235 $ g_tmp1, voff, nbf, 1, nbf ) 236 call ga_matmul_patch( 'n', 'n', 8.d0, 0.d0, ! Where does this factor come from? 237 $ g_tmp1, voff, nbf, 1, nbf, 238 $ g_movecs, 1, nbf, 1, nclosed, 239 $ g_tmp2, voff, nbf, 1, nclosed ) 240 call ga_dadd_patch( 1.d0, g_tmp2, voff, nbf, 1, nclosed, 241 $ 1.d0, g_ax, voff, nbf, 1, nclosed, 242 $ g_ax, voff, nbf, 1, nclosed ) 243c 244c Active-Virtual contribution 245c 246c 247 call ga_matmul_patch( 't', 'n', 1.d0, 0.d0, 248 $ g_movecs, voff, nbf, 1, nbf, 249 $ g_f2, 1, nbf, 1, nbf, 250 $ g_tmp1, voff, nbf, 1, nbf ) 251 call ga_zero(g_tmp2) 252 call ga_matmul_patch( 'n', 'n', 1.d0, 0.d0, 253 $ g_tmp1, voff, nbf, 1, nbf, 254 $ g_movecs, 1, nbf, aoff, aend, 255 $ g_tmp2, voff, nbf, aoff, aend ) 256 call ga_matmul_patch( 'n', 'n', 8.d0, 0.d0, 257 $ g_tmp2, voff, nbf, aoff, aend, 258 $ g_tmp, aoff, aend, aoff, aend, 259 $ g_tmp1, voff, nbf, aoff, aend ) 260 call ga_dadd_patch( 1.d0, g_tmp1, voff, nbf, aoff, aend, 261 $ 1.d0, g_ax, voff, nbf, aoff, aend, 262 $ g_ax, voff, nbf, aoff, aend ) 263c 264c Inactive-Active contributions 265c t iv 266c (2 - d) C.F[D ].C 267c 268 call ga_zero(g_tmp1) 269 call ga_zero(g_tmp2) 270 call ga_matmul_patch( 't', 'n', 1.d0, 0.d0, 271 $ g_movecs, aoff, aend, 1, nbf, 272 $ g_f3, 1, nbf, 1, nbf, 273 $ g_tmp1, aoff, aend, 1, nbf ) 274 call ga_matmul_patch( 'n', 'n', 1.d0, 0.d0, 275 $ g_tmp1, aoff, aend, 1, nbf, 276 $ g_movecs, 1, nbf, 1, nclosed, 277 $ g_tmp2, aoff, aend, 1, nclosed ) 278 call ga_zero(g_tmp1) 279 call diagfill_patch( g_tmp1, 2.d0, aoff, aend ) 280 call ga_dadd( 1.d0, g_tmp1, -1.d0, g_tmp, g_tmp1 ) 281 call ga_matmul_patch( 'n', 'n', 8.d0, 0.d0, 282 $ g_tmp1, aoff, aend, aoff, aend, 283 $ g_tmp2, aoff, aend, 1, nclosed, 284 $ g_tmp1, aoff, aend, 1, nclosed ) 285c 286c t ,av 287c -d C.F[D ].C 288c 289 call ga_zero(g_tmp2) 290 call ga_matmul_patch( 't', 'n', 1.d0, 0.d0, 291 $ g_movecs, aoff, aend, 1, nbf, 292 $ g_f4, 1, nbf, 1, nbf, 293 $ g_tmp2, aoff, aend, 1, nbf ) 294 call ga_matmul_patch( 'n', 'n', 8.d0, 1.d0, 295 $ g_tmp2, aoff, aend, 1, nbf, 296 $ g_movecs, 1, nbf, 1, nclosed, 297 $ g_tmp1, aoff, aend, 1, nclosed ) 298 call ga_dadd_patch( 1.d0, g_tmp1, aoff, aend, 1, nclosed, 299 $ 1.d0, g_ax, aoff, aend, 1, nclosed, 300 $ g_ax, aoff, aend, 1, nclosed ) 301c 302c 303c Inactive-Active, Inactive-Active contribution 304c (note this last section has zero contribution 305c in ROHF theory and has not been debugged against the ROHF 306c Hessian product) 307c 308c 309 call ga_matmul_patch( 't', 'n', 1.d0, 0.d0, 310 $ g_movecs, aoff, aend, 1, nbf, 311 $ g_f6, 1, nbf, 1, nbf, 312 $ g_tmp1, aoff, aend, 1, nbf ) 313 call ga_matmul_patch( 'n', 'n', 8.d0, 0.d0, 314 $ g_tmp1, aoff, aend, 1, nbf, 315 $ g_movecs, 1, nbf, 1, nclosed, 316 $ g_tmp2, aoff, aend, 1, nclosed ) 317 call ga_dadd_patch( 1.d0, g_tmp2, aoff, aend, 1, nclosed, 318 $ 1.d0, g_ax, aoff, aend, 1, nclosed, 319 $ g_ax, aoff, aend, 1, nclosed ) 320c 321c 322c 323 call ga_matmul_patch( 't', 'n', 1.d0, 0.d0, 324 $ g_movecs, aoff, aend, 1, nbf, 325 $ g_f5, 1, nbf, 1, nbf, 326 $ g_tmp1, aoff, aend, 1, nbf ) 327 call ga_matmul_patch( 'n', 'n', 8.d0, 0.d0, 328 $ g_tmp1, aoff, aend, 1, nbf, 329 $ g_movecs, 1, nbf, 1, nclosed, 330 $ g_tmp2, aoff, aend, 1, nclosed ) 331 call diagfill_patch( g_tmp1, 1.d0, aoff, aend ) 332 call ga_dadd( 1.d0, g_tmp1, -1.d0, g_tmp, g_tmp1 ) 333 call ga_matmul_patch( 'n', 'n', 1.d0, 0.d0, 334 $ g_tmp1, aoff, aend, aoff, aend, 335 $ g_tmp2, aoff, aend, 1, nclosed, 336 $ g_tmp1, aoff, aend, 1, nclosed ) 337 call ga_dadd_patch( 1.d0, g_tmp1, aoff, aend, 1, nclosed, 338 $ 1.d0, g_ax, aoff, aend, 1, nclosed, 339 $ g_ax, aoff, aend, 1, nclosed ) 340c 341c 342c 343 if (.not.ga_destroy(g_tmp1)) 344 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 345 if (.not.ga_destroy(g_tmp2)) 346 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 347 if (.not.ga_destroy(g_d1)) 348 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 349 if (.not.ga_destroy(g_d2)) 350 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 351 if (.not.ga_destroy(g_d3)) 352 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 353 if (.not.ga_destroy(g_d4)) 354 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 355 if (.not.ga_destroy(g_d5)) 356 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 357 if (.not.ga_destroy(g_d6)) 358 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 359 if (.not.ga_destroy(g_f1)) 360 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 361 if (.not.ga_destroy(g_f2)) 362 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 363 if (.not.ga_destroy(g_f3)) 364 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 365 if (.not.ga_destroy(g_f4)) 366 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 367 if (.not.ga_destroy(g_f5)) 368 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 369 if (.not.ga_destroy(g_f6)) 370 $ call errquit('mcscf_hessv_2e_ao: cannot destroy',0, GA_ERR) 371c 372c 373c 374 if (omcscfprof) call pstat_off(ps_hv2eao) 375 return 376 end 377 378 379