1 subroutine dft_gradients(rtdb) 2c 3c calculate energy gradients with respect to nuclear coordinates 4c modified from scf version for use in DFT gradients 5c 6c------------------------------------------------------------------------------ 7c ___ ___ ___ 8c dE \ dh(i,j) \ d(mn|op) \ dS(i,j) 9c -- = 2 > D(i,j) ------- + > P(m,n,o,p) -------- - 2 > Dw(i,j) ------- 10c dA / dA / dA / dA 11c --- --- --- 12c i,j ijkl i,j 13c 14c 15c dV(nuc-nuc) 16c + ---------- + exchange-correlation terms 17c dA 18c 19c------------------------------------------------------------------------------ 20c 1 21c P(i,j,k,l) = [2 D(i,j)D(k,l) - - (D(i,k)D(j,l) + D(i,l)D(j,k)) 22c 2 23c------------------------------------------------------------------------------ 24c 25c This version computes the pieces specific to DFT (XC on grid 26c and CD-fit) and call the standard grad_force() to do the rest 27c including writing out the results. 28c 29* 30* $Id$ 31* 32 implicit none 33#include "errquit.fh" 34c 35 integer rtdb 36c 37#include "mafdecls.fh" 38#include "global.fh" 39#include "rtdb.fh" 40#include "bas.fh" 41#include "geom.fh" 42#include "stdio.fh" 43#include "msgids.fh" 44#include "sym.fh" 45#include "cdft.fh" 46#include "util.fh" 47#include "dftps.fh" 48c 49c!!! BGJ test !!! 50#include "bgj.fh" 51 integer l_hess, k_hess, g_rhs(3,nw_max_atom), j 52c!!! BGJ test !!! 53 double precision zero, one, two 54 Parameter (zero=0.d0, one=1.d0, two=2.d0) 55c 56 integer ga_create_atom_blocked 57 external ga_create_atom_blocked 58 logical movecs_read_header, movecs_read,xc_rep_close 59 external movecs_read_header, movecs_read,xc_rep_close 60c integer noc(2) 61 integer nmo(2) 62 integer iga_dens(2), g_vecs(2), g_force 63 integer idum(4), ndum 64 double precision edum 65 Integer k_evals(2), l_evals(2) 66 double precision grad_norm, grad_max 67 external grad_norm, grad_max 68 character*255 title_vecs, basis_vecs 69 character*20 scftype_vecs 70 character*80 scftype 71 integer ifocc 72 logical status,frac_occ 73 integer me, nproc, max_sh_bf, max_at_bf, nat, max_sh_bfcd, 74 $ lforce, nactive, i, nbf_vecs, nsets, ispin, 75 $ max1e, max2e, mscratch_1e, mscratch_2e, 76 $ max2e3c, mscratch_2e3c, lbuf, lscratch, lsqa 77 integer l_force, k_force, l_occ, k_occ, l_act, k_act, 78 $ l_buf, k_buf, l_scr, k_scr, l_wdens, k_wdens, 79 $ l_cdcoef, i_cdcoef, ippp, isvec, lsvec, 80 $ ilo, ihi, 81 $ k_frc_2el, k_frc_xc, 82 $ l_frc_2el, l_frc_xc 83 integer lcntoce, icntoce, lcntobfr, icntobfr, 84 $ lcetobfr, icetobfr, lrdens_atom, irdens_atom, 85 $ nscr, lscr, iscr 86 integer ipoint, itmpm,ltmpm,g_tmp(2),lenvec,maavail 87 double precision charge, charge_nuc, rhffact, tol2e, onem, 88 , toll 89c 90c vdw 91 double precision dum 92 logical cgmin 93 logical disp 94 logical xc_chkdispauto 95 external xc_chkdispauto 96c 97c xdm 98 integer xdmdisp 99 integer nxdm 100 integer ixdm_v, ixdm_ml, lxdm_v, lxdm_ml 101 common /xdmd/ nxdm, ixdm_v, lxdm_v, ixdm_ml, lxdm_ml 102c 103 double precision fant_a,fant_d 104 parameter(toll=1.d-9) 105 logical oprint_force_comps 106c 107 logical has_frac_occ 108 external has_frac_occ 109c 110 nproc = ga_nnodes() 111 me=ga_nodeid() 112 oprint_force_comps = util_print('force components', print_debug) 113c 114c Print options 115c 116 if (.not. geom_ncent(geom, nat)) 117 $ call errquit('dft_gradient: could not get natoms',0, 118 & GEOM_ERR) 119c 120 if (.not. bas_nbf_cn_max(ao_bas_han, max_sh_bf)) 121 $ call errquit('dft_gradient: could not get max_sh_bf',0, 122 & BASIS_ERR) 123 max_at_bf = 0 124 do i = 1, nat 125 if (.not. bas_ce2bfr(ao_bas_han, i, ilo, ihi)) 126 $ call errquit('dft_gradient: bas_ce2bfr failed', i, 127 & BASIS_ERR) 128 max_at_bf = max(max_at_bf, ihi-ilo+1) 129 enddo 130c 131c use of scratch array in cdfit ... needs (3,max_at_bf) 132c 133 max_at_bf = max(max_at_bf,3) 134c 135 charge = rcharge 136 status = geom_nuc_charge(geom, charge_nuc) 137 if (.not.status)then 138 call errquit('dft_gradient: no nuclear charge',0, GEOM_ERR) 139 endif 140c 141c check for cgmin since it breaks cdfit 142c 143 if (.not.rtdb_get(rtdb,'dft:cgmin', mt_log, 1, cgmin)) 144 & cgmin=.false. 145 146c if (.not. rtdb_get(rtdb, 'dft:noc', mt_int, 2, noc)) 147c $ call errquit('dft_gradient: rtdb_get of noc failed', 0, 148c & RTDB_ERR) 149c 150c check if fractional occupation is on 151c 152 frac_occ = .false. 153 if (has_frac_occ(rtdb)) frac_occ = .true. 154c 155c allocate and initialize global and local memory 156c 157c mo-vectors 158c 159 if (ipol .eq. 1)then 160 g_vecs(1) = ga_create_atom_blocked (geom, ao_bas_han, 161 $ 'mo vectors') 162 else 163 g_vecs(1) = ga_create_atom_blocked (geom, ao_bas_han, 164 $ 'alpha mo vectors') 165 g_vecs(2) = ga_create_atom_blocked (geom, ao_bas_han, 166 $ 'beta mo vectors') 167 endif 168c 169c global density 170c 171 if (ipol .eq. 1)then 172 rhffact = two 173 else 174 rhffact = one 175 endif 176 if (ipol .eq. 1)then 177 iga_dens(1) = ga_create_atom_blocked (geom, ao_bas_han, 178 $ 'density matrix') 179 else 180 iga_dens(1) = ga_create_atom_blocked (geom, ao_bas_han, 181 $ 'alpha density matrix') 182 iga_dens(2) = ga_create_atom_blocked (geom, ao_bas_han, 183 $ 'beta density matrix') 184 endif 185c 186c forces on atoms (3xnat) 187c 188*ga:1:0 189 status = ga_create(mt_dbl, 3, nat, 'forces', 3, 0, g_force) 190 call ga_zero (g_force) 191c 192c local replication (separate for the different pieces) 193c 194 lforce = nat * 3 195 if (.not.ma_alloc_get(mt_dbl, lforce, 'forces',l_force, k_force)) 196 $ call errquit('could not allocate l_force',1, MA_ERR) 197 call dfill(lforce, 0.0d0, dbl_mb(k_force), 1) 198c 199 if (.not.ma_alloc_get(mt_dbl,lforce,'forces',l_frc_2el,k_frc_2el)) 200 $ call errquit('could not allocate l_frc_2el',1, MA_ERR) 201 call dfill(lforce, 0.0d0, dbl_mb(k_frc_2el), 1) 202c 203 if (.not.ma_alloc_get(mt_dbl,lforce,'forces',l_frc_xc,k_frc_xc)) 204 $ call errquit('could not allocate l_frc_xc',1, MA_ERR) 205 call dfill(lforce, 0.0d0, dbl_mb(k_frc_xc), 1) 206c 207c eigenvalues 208c 209 if (ipol .eq. 1)then 210 if (.not. ma_alloc_get(mt_dbl, nbf_ao, 'MO evals', l_evals(1), 211 $ k_evals(1))) 212 $ call errquit('dft_gradient: could not allocate l_evals',1, 213 & MA_ERR) 214 else 215 status = ma_alloc_get(mt_dbl, nbf_ao, 'alpha MO evals', 216 $ l_evals(1), k_evals(1)) 217 status = status .and. 218 $ ma_alloc_get(mt_dbl, nbf_ao, 'beta MO evals', 219 $ l_evals(2), k_evals(2)) 220 if (.not. status)then 221 call errquit('dft_gradient: could not allocate l_evals',1, 222 & MA_ERR) 223 endif 224 endif 225c 226c occupation numbers (not used, but necessary for movecs_read) 227c 228c should do k_occ for both spins, in case used at some point... 229c 230 if (.not. ma_alloc_get(mt_dbl, nbf_ao*ipol, 'occ. numbers', 231 $ l_occ, k_occ)) 232 $ call errquit('dft_gradient: could not allocate l_occ',1, 233 & MA_ERR) 234c 235c lookup table and list of active atoms 236c 237 if (.not. ma_alloc_get(MT_LOG, nat, 'active atoms', 238 $ l_act, k_act)) 239 $ call errquit('grad: could not allocate l_act',1, MA_ERR) 240 call grad_active_atoms(rtdb, nat, log_mb(k_act), nactive) 241c 242c get MO vectors from file 243c 244 if (.not. rtdb_cget(rtdb, 'dft:input vectors', 1, movecs_in)) 245 $ call errquit('dft_gradient: DFT MO vectors not defined',0, 246 & RTDB_ERR) 247 status = movecs_read_header(movecs_in, title_vecs, basis_vecs, 248 $ scftype_vecs, nbf_vecs, nsets, nmo, 2) 249c 250c ipol - number of spin channels (RKS=1, ROKS=2, UKS=2) 251c nsets - number of sets of vectors (RKS=1, ROKS=1, UKS=2) 252c 253 if (.not. rtdb_cget(rtdb, 'dft:scftype', 1, scftype)) 254 $ call errquit('dft_gradient: DFT scftype not defined',0, 255 & RTDB_ERR) 256 if (scftype.eq.'RHF') then 257 if (ipol .ne. nsets .or. ipol.ne.1)then 258 write (LuOut,*) 'dft_gradient: ERROR ipol, nsets:',ipol,nsets 259 call errquit('dft_gradient: ERROR ipol, nsets disagree',2, 260 & INPUT_ERR) 261 endif 262 elseif (scftype.eq.'ROHF') then 263 if (nsets.ne.1.or.ipol.ne.2) then 264 write (LuOut,*) 'dft_gradient: ERROR ipol, nsets:',ipol,nsets 265 call errquit('dft_gradient: ERROR ipol, nsets disagree',2, 266 & INPUT_ERR) 267 endif 268 elseif (scftype.eq.'UHF') then 269 if (nsets.ne.2.or.ipol.ne.2) then 270 write (LuOut,*) 'dft_gradient: ERROR ipol, nsets:',ipol,nsets 271 call errquit('dft_gradient: ERROR ipol, nsets disagree',2, 272 & INPUT_ERR) 273 endif 274 else 275 call errquit('dft_gradient: illegal scftype',0,UERR) 276 endif 277c 278c Should check much more info than just nbf for consistency 279c 280c 281c get mo eigenvectors 282c 283 if (nbf_ao .ne. nbf_vecs)then 284 write(LuOut,*)'dft_gradient movecs output = ',movecs_in 285 call errquit('dft_gradient: could not read mo vectors',911, 286 & DISK_ERR) 287 else 288 status = .true. 289 do ispin = 1, ipol 290 status = status .and. 291 $ movecs_read(movecs_in, min(ispin,nsets), 292 & dbl_mb(k_occ+(ispin-1)*nbf_ao), 293 $ dbl_mb(k_evals(ispin)), g_vecs(ispin)) 294 enddo 295 endif 296c 297 if (.not.status)then 298 write(LuOut,*)'dft_gradient movecs output = ',movecs_in 299 call errquit('dft_gradient: could not read mo vectors',917, 300 & DISK_ERR) 301 endif 302c 303 if(frac_occ) then 304c 305c fractional occupation, therefore check new nocs 306c 307 if (.not. MA_Push_Get(MT_Dbl, nbf_ao, 'tmpm', ltmpm, itmpm)) 308 & call errquit('dftgforce: failed to alloc tmpm',0, MA_ERR) 309 rhffact = one 310c 311 do ispin=1,ipol 312 g_tmp(ispin) = ga_create_atom_blocked(geom, ao_bas_han, 313 & 'frac vecs') 314 call ga_zero(g_tmp(ispin)) 315 ipoint=k_occ+(ispin-1)*nbf_ao-1 316c 317 do i = ga_nodeid()+1, nbf_ao, ga_nnodes() 318 call get_col(g_vecs(ispin), nbf_ao, i, DBL_MB(itmpm)) 319 call dscal(nbf_ao, dbl_mb(i+ipoint), DBL_MB(itmpm), 1) 320 call put_col(g_tmp(ispin), nbf_ao, i, DBL_MB(itmpm)) 321 enddo 322 do i=1,nbf_ao 323 if(dbl_mb(ipoint+i).ge.toll) noc(ispin)=i 324 enddo 325 enddo 326 if (.not.ma_pop_stack(ltmpm)) 327 & call errquit('dftg_force: cannot pop stack',0, MA_ERR) 328 329 else 330 do ispin=1,ipol 331 g_tmp(ispin)=g_vecs(ispin) 332 enddo 333 endif 334c 335 do ispin = 1, ipol 336c 337c dens = vecs*vecs 338c 339 if (odftps) call pstat_on(ps_dgemm) 340 call ga_dgemm('n', 't', nbf_ao, nbf_ao, noc(ispin), rhffact, 341 $ g_tmp(ispin), g_vecs(ispin), 0.0d0, iga_dens(ispin)) 342 if (odftps) call pstat_off(ps_dgemm) 343 call ga_symmetrize(iga_dens(ispin)) 344c 345c free temporary arrays 346c 347 if(frac_occ) then 348 if(.not.ga_destroy (g_tmp(ispin))) call 349 * errquit('dftg_force: could not gadestr gtmp',ispin, 350 & GA_ERR) 351 endif 352 enddo !ispin 353c 354c Pre-compute mapping vectors 355c 356 if (.not.ma_push_get 357 $ (mt_int,nat*2,'cntoce map',lcetobfr,icetobfr)) 358 $ call errquit('dft_scf:push_get failed', 13, MA_ERR) 359 if (.not.ma_push_get 360 $ (mt_int,nshells_ao,'cntoce map',lcntoce,icntoce)) 361 $ call errquit('dft_scf:push_get failed', 13, MA_ERR) 362 if (.not.ma_push_get 363 $ (mt_int,nshells_ao*2,'cntoce map',lcntobfr,icntobfr)) 364 $ call errquit('dft_scf:push_get failed', 13, MA_ERR) 365c 366 call build_maps(ao_bas_han, int_mb(icntoce), int_mb(icntobfr), 367 $ int_mb(icetobfr), nat, nshells_ao) 368 if (.not.ma_chop_stack(lcntoce)) 369 $ call errquit('dft_gradient: cannot pop stack',0, MA_ERR) 370c 371c Pre-compute reduced density matrices over atoms 372c 373 if (.not.ma_push_get(mt_dbl,ipol*nat*nat,'rdens_atom', 374 $ lrdens_atom,irdens_atom)) 375 $ call errquit('dft_scf: cannot allocate rdens_atom',0, MA_ERR) 376 call dfill(ipol*nat*nat, 0.0d0, dbl_mb(irdens_atom), 1) 377 nscr = nbf_ao_mxnbf_ce*nbf_ao_mxnbf_ce 378 if (.not.ma_push_get(mt_dbl,nscr,'scr',lscr,iscr)) 379 $ call errquit('dft_scf: cannot allocate scr',0, MA_ERR) 380 call util_ga_mat_reduce(nbf_ao, nat, int_mb(icetobfr), 381 $ iga_dens, ipol, dbl_mb(irdens_atom), 382 $ 'absmax', dbl_mb(iscr), nbf_ao_mxnbf_ce,.true.) 383 if (.not.ma_pop_stack(lscr)) 384 $ call errquit('dft_scf: cannot pop stack',0, MA_ERR) 385c 386 if (ipol .eq. 2)status = ga_destroy(g_vecs(2)) 387 status = ga_destroy(g_vecs(1)) 388c 389 if (.not.status)then 390 call errquit('dft_gradient: could not destroy g_eigen_diag',1, 391 & GA_ERR) 392 endif 393c 394 status = ma_free_heap(l_occ) 395 if (ipol .eq. 2)then 396 status = ma_free_heap (l_evals(2)) 397 endif 398 status = ma_free_heap (l_evals(1)) 399c 400 if (CDFIT.and.(.not.cgmin)) then 401c 402c determine memory requirements for integral gradients 403c 404 call int_mem(max1e, max2e, mscratch_1e, mscratch_2e) 405 call int_mem_2e3c(max2e3c, mscratch_2e3c) 406 lbuf = max(max1e, max2e) 407 lbuf = max(lbuf, max2e3c) + 500 408 lscratch = max(mscratch_1e, mscratch_2e) 409 lscratch = max(lscratch, mscratch_2e3c) 410c 411c one-electron contribution 412c buffers for one electron integral derivatives 413c 414 status = ma_push_get(mt_dbl, lbuf, 'deriv buffer', l_buf, k_buf) 415 if(.not.status) then 416 maavail=MA_inquire_avail(mt_dbl) 417 call errquit('dft_gradient: could not allocate buffer', 418 & 8*(lbuf-maavail), MA_ERR) 419 endif 420c 421 status = ma_push_get(mt_dbl, lscratch, 'deriv scratch', 422 $ l_scr, k_scr) 423 if (.not.status) 424 $ call errquit('dft_gradient: could not allocate scratch',1, 425 & MA_ERR) 426c 427c allocate local density matrix block 428c 429 lsqa = max_at_bf * max_at_bf 430c 431 status = ma_push_get(mt_dbl, lsqa, 'local_w_density', 432 $ l_wdens, k_wdens) 433 if (.not.status)call errquit('could not allocate l_dens',1, 434 & MA_ERR) 435c 436c store total DM in ga_dens(1) 437c 438 if (ipol .eq. 2)then 439 call ga_dadd (one,iga_dens(1),one,iga_dens(2),iga_dens(1)) 440 endif 441c 442c define threshold for Schwarz screening (same as in DFT) 443c 444 tol2e=10.d0**(-itol2e) 445c 446c 447c compute 3 center coulomb derivative 448c 449c Determine the characteristics of the AO and CD Gaussian basis sets. 450c 451 if (.not. ma_push_get(mt_dbl, nbf_cd, 'CD coef', 452 $ l_cdcoef, i_cdcoef)) 453 $ call errquit('dft_gradient: could not alloc CD coef',0, 454 & MA_ERR) 455c 456 if(.not.bas_nbf_cn_max(cd_bas_han, max_sh_bfcd)) 457 $ call errquit('dftg_force: basnbfcdmax broken?',0, 458 & BASIS_ERR) 459 lenvec=max(3*max_sh_bfcd,max_sh_bf*max_sh_bf) 460 if (.not. MA_Push_get(MT_DBL, lenvec, 'svec', 461 $ lsvec, isvec)) 462 $ call errquit('dftg_force: could not alloc svec',0, MA_ERR) 463 ippp=k_wdens 464c 465 if (odftps) call pstat_on(ps_vcoul) 466 call xc_rep_init(rtdb, geom, ao_bas_han,iga_dens,iga_dens, 467 & nbf_ao,ipol,.true.,.true.) 468 call dftg_cdfit(geom,ao_bas_han, cd_bas_han, 469 $ nbf_cd, nat, tol2e, dbl_mb(k_scr), 470 $ lscratch, dbl_mb(k_buf), lbuf, 471 $ dbl_mb(isvec), dbl_mb(ippp), max_sh_bf, 472 $ iga_dens, dbl_mb(k_frc_2el), 473 $ DBL_MB(i_cdcoef), oskel) 474 if(.not.xc_rep_close(rtdb, nbf_ao,ipol,ipol, 475 D iga_dens,iga_dens,.true.)) call 476 . errquit(' dftggrad: xcrepclose failed ',0, 0) 477c 478 call ga_dgop(msg_grad_2el, dbl_mb(k_frc_2el), lforce, '+') 479 if (odftps) call pstat_off(ps_vcoul) 480c 481 if (.not.ma_chop_stack(l_buf)) 482 $ call errquit('dft_gradient: cannot chop stack',0, MA_ERR) 483c 484c restore alpha DM in g_dens(1) 485c 486 if (ipol .eq. 2)then 487 onem = -1.d0 488 call ga_dadd(one, iga_dens(1), onem, iga_dens(2), 489 & iga_dens(1)) 490 endif 491 endif ! cdfit 492c 493c get exchange-correlation contribution to the gradient 494c 495c 496c$$$ write(LuOut,*) ' BEFORE CALL TO GETXC' 497c$$$ call ga_print(iga_dens(1)) 498c$$$ call ga_print(iga_dens(2)) 499c$$$ write(LuOut,*) ' nactive ' 500c$$$ call output(dbl_mb(irdens_atom), 1, nat, 1, nat, nat, nat, 1) 501c$$$ write(LuOut,*) (int_mb(icntoce+i),i=0,nshells_ao-1) 502c$$$ write(LuOut,*) (int_mb(icntobfr+i),i=0,2*nshells_ao-1) 503c$$$ write(LuOut,*) (int_mb(icetobfr+i),i=0,2*nat-1) 504 if (odftps) call pstat_on(ps_xc) 505 call dftg_getxc(rtdb, nat, iga_dens, dbl_mb(k_frc_xc), 506 $ log_mb(k_act), nactive, 507 $ dbl_mb(irdens_atom), int_mb(icetobfr)) 508 call ga_dgop(msg_grad_xc, dbl_mb(k_frc_xc), lforce, '+') 509 if(oprint_force_comps.and.me.eq.0)then 510 write(luout,2200) 511 $ 'XC gradient', 512 $ ((dbl_mb(k_frc_xc+i-1+3*(j-1)),i=1,3),j=1,nat) 513 write(luout,2200) 514 $ 'CD gradient', 515 $ ((dbl_mb(k_frc_2el+i-1+3*(j-1)),i=1,3),j=1,nat) 516 endif 517 2200 format(A/,1000(3(1x,F12.6),/)) 518 519 if (odftps) call pstat_off(ps_xc) 520c 521c vdW bit 522c 523 if (.not.rtdb_get(rtdb, 'dft:disp', mt_log, 1, disp)) 524 & disp=.false. 525 if(disp.or.xc_chkdispauto()) 526 & call xc_vdw(rtdb,geom,dum,dbl_mb(k_frc_xc), 'forces') 527c 528c xdm bit 529c 530 if (.not.rtdb_get(rtdb, 'dft:xdm', mt_int, 1, xdmdisp)) 531 & xdmdisp=0 532 if (xdmdisp.ne.0) then 533 call xc_xdm(rtdb,iga_dens,idum,nat,ndum,edum, 534 & dbl_mb(k_frc_xc),dbl_mb(ixdm_v), 535 & dbl_mb(ixdm_ml),'forces') 536 if (.not. rtdb_put(rtdb,'dft:xdmsave', mt_log, 1, .false.)) 537 $ call errquit('dft_gradient: cannot rtdb_put',0, RTDB_ERR) 538 call xc_xdm_cleanup(rtdb) 539 endif 540c 541c add Bonacic-Fantucci repulsive term 542c 543 if (.not.rtdb_get(rtdb, 'dft:fant_d', mt_dbl, 1, 544 & fant_d)) fant_d=-1d0 545 if (.not.rtdb_get(rtdb, 'dft:fant_a', mt_dbl, 1, 546 & fant_a)) fant_a=-1d0 547 if(fant_a.ne.-1d0.and.fant_d.ne.-1d0) 548 A call dftg_fant(geom,nat,fant_a,fant_d, 549 A dbl_mb(k_frc_xc)) 550c 551 if (ga_nodeid() .eq. 0)then 552 status = rtdb_parallel (.false.) 553 do i = 0, lforce-1 554 dbl_mb(k_force+i) = dbl_mb(k_frc_2el+i) + dbl_mb(k_frc_xc+i) 555 enddo 556 if (.not. rtdb_put(rtdb, 'dft:cd+xc gradient', mt_dbl, 557 $ lforce, dbl_mb(k_force))) call errquit 558 $ ('dft_gradient: failed storing cd+xc gradient',0, 559 & UNKNOWN_ERR) 560c 561 status = rtdb_parallel (.true.) 562 endif 563 if(oprint_force_comps.and.me.eq.0)then 564 do 31 i = 1, nat 565 write (luout,2000) i, 566 & (dbl_mb(k_frc_xc+3*(i-1)+j),j=0,2) 567 2000 format (1X,I3,3(1X,F10.6)) 568 31 continue 569 endif 570c 571 call ga_sync() 572c 573 if (.not.ma_pop_stack(lrdens_atom)) 574 $ call errquit('dft_gradient: cannot pop stack',0, MA_ERR) 575 if (.not.ma_pop_stack(lcetobfr)) 576 $ call errquit('dft_gradient: cannot pop stack',0, MA_ERR) 577c 578c!!! BGJ test !!! 579c 580c store total DM in ga_dens(1) 581c 582 if (ipol .eq. 2)then 583 call ga_dadd (1d0,iga_dens(1),1d0,iga_dens(2),iga_dens(1)) 584 endif 585c 586c J hesssian test calculation done by setting bgj:j_hessian 587c to true 588c 589 if (.not. rtdb_get(rtdb, 'bgj:j_hessian', mt_log, 590 & 1, status)) status = .false. 591 if (status) then 592 call schwarz_tidy() 593 call intd_terminate() 594 call int_init(rtdb, 1, ao_bas_han) 595 call schwarz_init (geom, ao_bas_han) 596 call int_terminate() 597 if (CDFIT) then 598 nmo(1) = ao_bas_han 599 nmo(2) = cd_bas_han 600 call intdd_init(rtdb,2,nmo) 601 else 602 call intdd_init(rtdb,1,ao_bas_han) 603 endif 604 605 status = MA_push_get(MT_DBL, 3*nat*3*nat, 606 & 'j hessian', l_hess, k_hess) 607 if (.not.status) 608 & call errquit('dft_gradients: could not alloc j hessian', 609 & 1, MA_ERR) 610 call dfill(9*nat*nat, 0.0d0, dbl_mb(k_hess), 1) 611 if (bgj_print() .gt. 0) 612 & write(LuOut,*)'*** In dft_gradients: calling j_hessian' 613 call j_hessian(iga_dens, log_mb(k_act), nactive, 614 & dbl_mb(k_hess)) 615 status = MA_pop_stack(l_hess) 616 if (.not.status) call 617 & errquit('dft_gradients: could not pop j hessian', 618 & 1, MA_ERR) 619 620 call schwarz_tidy() 621 call intdd_terminate() 622 623 endif 624c 625c J CPKS RHS test calculation done by setting bgj:j_cpks_rhs 626c to true 627c 628 if (.not. rtdb_get(rtdb, 'bgj:j_cpks_rhs', mt_log, 629 & 1, status)) status = .false. 630 if (status) then 631 632 call schwarz_tidy() 633 call intd_terminate() 634 635 call int_init(rtdb, 1, ao_bas_han) 636 call schwarz_init (geom, ao_bas_han) 637 call int_terminate() 638 if (CDFIT) then 639 nmo(1) = ao_bas_han 640 nmo(2) = cd_bas_han 641 call intd_init(rtdb,2,nmo) 642 else 643 call intd_init(rtdb,1,ao_bas_han) 644 endif 645c !!! Do this to be consistent with DFT gradient 646c!!! call int_app_set_no_texas(rtdb) 647c 648c Allocate and initialize temp GA's for RHS 649c 650 if (bgj_print() .gt. 0) 651 & write(*,*)'*** j cpks rhs test: nactive =',nactive 652 if (nat.gt.100) 653 & call errquit('dft_gradients: dimension error in test',0, 654 & UNKNOWN_ERR) 655 do i = 1, nat 656 if (log_mb(k_act+i-1)) then 657 g_rhs(1,i) = ga_create_atom_blocked 658 & (geom, ao_bas_han, 'CPKS RHS test x') 659 g_rhs(2,i) = ga_create_atom_blocked 660 & (geom, ao_bas_han, 'CPKS RHS test y') 661 g_rhs(3,i) = ga_create_atom_blocked 662 & (geom, ao_bas_han, 'CPKS RHS test z') 663 call ga_zero(g_rhs(1,i)) 664 call ga_zero(g_rhs(2,i)) 665 call ga_zero(g_rhs(3,i)) 666 endif 667 enddo 668 669 if (bgj_print() .gt. 0) 670 & write(LuOut,*)'*** In dft_gradients: calling j_cpks_rhs' 671 call j_cpks_rhs(iga_dens, log_mb(k_act), nactive, g_rhs) 672 673 do i = 1, nat 674 if (log_mb(k_act+i-1)) then 675 do j = 1, 3 676 if (.not.ga_destroy(g_rhs(j,i))) then 677 call errquit('j_cpks_rhs: problem destroying ga',1, 678 & GA_ERR) 679 endif 680 enddo 681 endif 682 enddo 683 endif 684c!!! BGJ test !!! 685c 686 status = ma_free_heap (l_act) 687 status = ma_free_heap (l_frc_xc) 688 status = ma_free_heap (l_frc_2el) 689 status = ma_free_heap (l_force) 690 if (ipol .eq. 2) status = ga_destroy (iga_dens(2)) 691 status = ga_destroy (iga_dens(1)) 692 status = ga_destroy (g_force) 693c 694 return 695 end 696c 697 subroutine dftg_fant(geom,natoms,a,d,forces) 698C Bonacic-Kouteck, V.; Cespiva, L.; Fantucci, P.; Pittner, J.; Kouteck, J. J. 699C Chem. Phys. 1994, 98, 490. 700C Mitric, R.; Hartmann, M.; Stanca, B.; Bonacic-Koutecky, V.; Fantucci, P.; 701C J. Phys. Chem. A.; (Article); 2001; 105(39); 8892-8905 702C The core-core repulsion has been corrected according to 703C (CC(rij) ) l/rij + D exp(-arij)). 704C Constants D and a obtained for 1e-RECP's from the fitting procedure 705C described in Appendix A necessery for corection of core-core potential: 706C DAg-Ag ) 1619.887392, aAg-Ag ) 2.49630105, 707C DAu-Au ) 1911.131090, aAu-Au ) 2.46590129, 708C DAu-Ag ) 1765.509532, aAu-Ag ) 2.48110117. 709 implicit none 710#include "errquit.fh" 711#include "geom.fh" 712#include "global.fh" 713#include "stdio.fh" 714 integer geom,natoms 715 double precision a,d,forces(3,*) 716c 717 integer i,j,B 718 double precision qi,qj 719 double precision xB(3),xj(3),rBj,ffant,drbjdxb 720 character*16 tag 721c 722 do B=1,natoms 723 if (.not. geom_cent_get(geom, B, tag, 724 & xB, qi))call errquit 725 & ('grid_acc_def: geom_cent_get failed', 0, GEOM_ERR) 726 do j=1,natoms 727 if(B.ne.j) then 728 if (.not. geom_cent_get(geom, j, tag, 729 & xj, qj))call errquit 730 & ('grid_acc_def: geom_cent_get failed', 0, GEOM_ERR) 731 rBj=sqrt((xB(1)-xj(1))**2+(xB(2)-xj(2))**2+ 732 + (xB(3)-xj(3))**2) 733 ffant=d*exp(-a*rBj)*(-a) 734 do i=1,3 735 drBjdxb=1d0/rBj*(xB(i)-xj(i)) 736 forces(i,B)=forces(i,B) + ffant*drBjdxB 737 enddo 738 endif 739 740 enddo 741 enddo 742! if(ga_nodeid().eq.0) 743! W write(luout,*) 744! W ' Bonacic-Koutecky Fantucci-Repulsive Term ',dft_fant 745 return 746 end 747