1* $Id$ 2 subroutine atscf_ecp(geom,basis,tag,hatom,nbas, 3 & zeta,cont,ucont,nz) 4* 5*.... computes the ecp contribution to the atomic scf potential 6* 7 implicit none 8#include "errquit.fh" 9#include "mafdecls.fh" 10#include "bas.fh" 11#include "geom.fh" 12#include "stdio.fh" 13c 14c::-passed 15 integer geom ! [input] geometry handle 16 integer basis ! [input] basis set handle 17 character*16 tag ! [input] tag for atom 18 double precision hatom(*) ! [output] ecp contribution to the V 19 integer nbas(4) ! [input] atomic scf internal structure of 20* number of basis functions per sym 21* or shell type s, p, d, f 22 integer nz ! [input] number of exponents and coefs 23 double precision zeta(nz) ! [input] exponents 24 double precision cont(nz) ! [input] scaled atomic scf coefs 25 double precision ucont(nz)! [input] unscaled coefs (from bas) 26c::-local 27 integer bgeom ! geometry handle that basis was loaded with 28 integer nat ! number of atoms 29 integer icent ! center index 30 character*16 gtag ! dummy geometry tag 31 character*16 btag ! dummy basis tag 32 integer clo ! first non-uniuqe contraction on atom 33 integer chi ! last non-uniuqe contraction on atom 34 integer i,j ! dummy index 35 integer sz_hatom ! size of hatom array 36 integer ecpis ! ecp basis handle 37c 38cMG_start pseudopotentials do not work for periodic systems yet 39c 40 integer isys 41c 42 if (.not. geom_systype_get(geom, isys)) call errquit 43 $ ('atscf_ecp: systype?', 0, UNKNOWN_ERR) 44 if (isys.gt.0) return 45cMG_end 46c 47 gtag = ' ' 48 btag = ' ' 49c 50c compare geometry handles input vs. what bas used to read basis info 51c 52 if (.not.bas_geom(basis,bgeom)) call errquit 53 & ('atscf_ecp: bas_geom failed ',911, BASIS_ERR) 54 if (geom.ne.bgeom) call errquit 55 & ('atscf_ecp: geom .ne. bgeom Fatal mismatch error',911, 56 & GEOM_ERR) 57c 58c determine the center index in the global picture 59c 60 if (.not.geom_ncent(geom,nat)) call errquit 61 & ('atscf_ecp: geom_ncent failed ',911, GEOM_ERR) 62c 63 do icent = 1,nat 64 if (.not.geom_cent_tag(geom, icent, gtag)) call errquit 65 & ('atscf_ecp: geom_cent_tag failed ',911, GEOM_ERR) 66 if (gtag.eq.tag) goto 00001 67 enddo 6800001 continue 69*debug: write(6,*)'tags 1', tag, gtag 70c 71c determine size of hatom (primitive triangles in s,p,d,f) 72c 73 sz_hatom = 0 74 do i = 1,4 75 j = nbas(i) 76 sz_hatom = sz_hatom + j*(j+1)/2 77 enddo 78*debug: write(6,*)' atscf_ecp: size of hatom is ',sz_hatom 79 call dfill(sz_hatom,0.0d00,hatom,1) 80c 81c check if tag is an ecp center ? If not do not compute anything 82c 83 if (.not.geom_ecp_get(geom,icent)) return 84c 85c get ecpis .. ecp basis set handle 86c 87 if (.not.bas_get_ecp_handle(basis,ecpis)) call errquit 88 & ('atscf_ecp: bas_get_ecp_handle failed ',911, BASIS_ERR) 89c 90c gather info for call to ecp_integral 91c 92 if (.not.bas_ce2cnr(ecpis,icent,clo,chi)) 93 & call errquit('atscf_ecp:bas_ce2cnr failed',911, BASIS_ERR) 94c Check if center is in the given ecp basis set 95 if (.not.(bas_isce(geom,ecpis,icent))) then 96 write(luout,*)' the lexical center ',icent, ' with tag ', 97 & gtag,' is not in the ecp basis ' 98 call errquit('atscf_ecp: tag mismatch fatal error',911, 99 & BASIS_ERR) 100 endif 101 call atscf_ecp_setup_hatom(geom,basis,ecpis, 102 & hatom,sz_hatom,nbas,zeta,cont,ucont,nz,icent) 103 end 104 subroutine atscf_ecp_setup_hatom(geom,basis,ecpis, 105 & hatom,sz_hatom,nbas,zeta,cont,ucont,nz,icent) 106 implicit none 107#include "errquit.fh" 108* 109* routine that sets up pointers for atomic ecp computation 110* 111#include "mafdecls.fh" 112#include "bas.fh" 113c::-passed 114 integer geom ! [input] geometry handle 115 integer basis ! [input] basis set handle 116 integer ecpis ! [input] ecp basis set handle 117 integer nbas(4) ! [input] atomic scf internal structure of 118* number of basis functions per sym 119* or shell type s, p, d, f 120 integer nz ! [input] number of exponents and coefs 121 double precision zeta(nz) ! [input] exponents 122 double precision cont(nz) ! [input] scaled atomic scf coefs 123 double precision ucont(nz)! [input] unscaled coefs (from bas) 124 integer sz_hatom ! [input] size of hatom array 125 double precision hatom(sz_hatom) ! [output] ecp part of V 126 integer icent ! [input] lexical index of center 127c::-local 128 integer sz_zeta ! size of exponent and coef arrays for ecp info 129 integer lmaxbs ! high angular momentum in basis set 130 integer lmax_all_ecp ! high angular momentum in ecp basis 131 integer lscr, lbuf ! size of scratch and buffer arrays 132 integer h_zetac, k_zetac ! MA handle and pointer ecp exponents 133 integer h_coefc, k_coefc ! MA handle and pointer ecp coefs 134 integer h_nprimc, k_nprimc ! MA handle and pointer ecp n_prim_c 135 integer h_ncoefc, k_ncoefc ! MA handle and pointer ecp n_coef_c 136 integer h_indc, k_indc ! MA handle and pointer ecp ind_c 137 integer h_scr, k_scr ! MA handle and pointer scratch 138 integer h_buf, k_buf ! MA handle and pointer buffer 139c 140c... allocate space for info on ecp functions 141c.. n_zeta_c number of exponents, need a max 142c.1. zeta_c = exponents of ecp (flat for all ecp info) 143c.2. coef_c = contraction coeefs (flat for all ecp info) 144c.3. n_prim_c = info about primitives rval, type and center 145* (0:4,-1:lmax,necp=1) 146c.4. n_coef_c = info about prims summed over rval just type 147* and center (-1:lmax,necp=1) 148c.5. ind_c = first coeff/exponent for each type (-1:lmax,necp=1) 149c.. n_zeta_C = total number ecp exponents summed over centers 150c.. l_C = max projector on center C (necp=1) 151c.. n_c = number of ecp centers necp=1=n_c 152c.. l_ecp_max = max l val of projector on any center (lmax) 153c 154c.. only need 5 arrays use arbitrary length for zeta_c and coef_c 155* with max for test 156c.. get lmax in ecp basis and use this. 157c.. 158c get ecp primitive exponents in ma 159c 160 sz_zeta = 500 ! should be okay for most atoms 161 162c.. allocate memory for exponents 163 if (.not. 164 & ma_alloc_get(mt_dbl,sz_zeta, 165 & 'atscf ecp exponents', 166 & h_zetac,k_zetac)) call errquit 167 & ('atscf_ecp_setup_hatom: ma_alloc get failed',911, MA_ERR) 168 call dfill(sz_zeta,0.0d00,dbl_mb(k_zetac),1) 169 170c.. allocate memory for coeffs 171 if (.not. 172 & ma_alloc_get(mt_dbl,sz_zeta, 173 & 'atscf ecp coeffs', 174 & h_coefc,k_coefc)) call errquit 175 & ('atscf_ecp_setup_hatom: ma_alloc get failed',911, MA_ERR) 176 call dfill(sz_zeta,0.0d00,dbl_mb(k_coefc),1) 177 178c.. get lmax for all ecps 179 if (.not.bas_high_angular(ecpis,lmax_all_ecp)) call errquit 180 & ('atscf_ecp_setup_hatom: bas_high_angular failed',911, 181 & BASIS_ERR) 182 183c.. get lmax for all basis set 184 if (.not.bas_high_angular(basis,lmaxbs)) call errquit 185 & ('atscf_ecp_setup_hatom: bas_high_angular failed',911, 186 & BASIS_ERR) 187 188c.. allocate memory for n_prim_c array (see ecp_integral for def) 189 if (.not. 190 & ma_alloc_get(mt_int,(5*(lmax_all_ecp+2)), 191 & 'atscf ecp nprimc', 192 & h_nprimc,k_nprimc)) call errquit 193 & ('atscf_ecp_setup_hatom: ma_alloc get failed',911, 194 & MA_ERR) 195 call ifill((5*(lmax_all_ecp+2)),0,int_mb(k_nprimc),1) 196 197c.. allocate memory for n_coef_c array (see ecp_integral for def) 198 if (.not. 199 & ma_alloc_get(mt_int,(lmax_all_ecp+2), 200 & 'atscf ecp ncoefc', 201 & h_ncoefc,k_ncoefc)) call errquit 202 & ('atscf_ecp_setup_hatom: ma_alloc get failed',911, MA_ERR) 203 call ifill((lmax_all_ecp+2),0,int_mb(k_ncoefc),1) 204 205c.. allocate memory for ind_c array (see ecp_integral for def) 206 if (.not. 207 & ma_alloc_get(mt_int,(lmax_all_ecp+2), 208 & 'atscf ecp nindc', 209 & h_indc,k_indc)) call errquit 210 & ('atscf_ecp_setup_hatom: ma_alloc get failed',911, MA_ERR) 211 call ifill((lmax_all_ecp+2),0,int_mb(k_indc),1) 212 213c.. allocate memory for ecp integral buffer array 214 lbuf = (lmaxbs+1)*(lmaxbs+2)/2 215 lbuf = lbuf*lbuf 216 if (.not. 217 & ma_alloc_get(mt_dbl,lbuf, 218 & 'atscf ecp integrals', 219 & h_buf,k_buf)) call errquit 220 & ('atscf_ecp_setup_hatom: ma_alloc get failed',911, MA_ERR) 221* does not need to be zeroed 222 223c.. allocate memory for scratch array used by ecp code 224 lscr = max(lbuf*10,20000) 225 if (.not. 226 & ma_alloc_get(mt_dbl,lscr, 227 & 'atscf ecp scr', 228 & h_scr,k_scr)) call errquit 229 & ('atscf_ecp_setup_hatom: ma_alloc get failed',911, MA_ERR) 230* does not need to be zeroed 231 232c.. call routine that does computation 233 call atscf_ecp_build_hatom(geom,basis,ecpis, 234 & hatom,sz_hatom,nbas,zeta,cont,ucont,nz,icent, 235 & sz_zeta,dbl_mb(k_zetac),dbl_mb(k_coefc), 236 & lmax_all_ecp, 237 & int_mb(k_nprimc),int_mb(k_ncoefc),int_mb(k_indc), 238 & dbl_mb(k_scr),lscr,dbl_mb(k_buf),lbuf) 239 240c.. free allocated memory 241 if (.not.ma_free_heap(h_zetac)) call errquit 242 & ('atscf_ecp: ma_free_heap failed for h_zetac',911, MEM_ERR) 243 if (.not.ma_free_heap(h_coefc)) call errquit 244 & ('atscf_ecp: ma_free_heap failed for h_coefc',911, MEM_ERR) 245 if (.not.ma_free_heap(h_nprimc)) call errquit 246 & ('atscf_ecp: ma_free_heap failed for h_nprimc',911, MEM_ERR) 247 if (.not.ma_free_heap(h_ncoefc)) call errquit 248 & ('atscf_ecp: ma_free_heap failed for h_ncoefc',911, MEM_ERR) 249 if (.not.ma_free_heap(h_indc)) call errquit 250 & ('atscf_ecp: ma_free_heap failed for h_indc',911, MEM_ERR) 251 if (.not.ma_free_heap(h_buf)) call errquit 252 & ('atscf_ecp: ma_free_heap failed for h_buf',911, MEM_ERR) 253 if (.not.ma_free_heap(h_scr)) call errquit 254 & ('atscf_ecp: ma_free_heap failed for h_scr',911, MEM_ERR) 255 end 256 subroutine atscf_ecp_build_hatom(geom,basis,ecpis, 257 & hatom,sz_hatom,nbas,zeta,cont,ucont,nz,icent, 258 & sz_zetac,zeta_c,coef_c,lmax_all,n_prim_c, 259 & n_coef_c, ind_c, scr,lscr, ecp_ints, l_ecp_ints) 260 implicit none 261#include "errquit.fh" 262* 263* this routine fills pointer and coeff arrays with appropriate 264* ecp information and then generates the ecp integrals 265* 266#include "stdio.fh" 267#include "mafdecls.fh" 268#include "bas.fh" 269#include "nwc_const.fh" 270#include "basP.fh" 271#include "bas_ibs_dec.fh" 272#include "bas_exndcf_dec.fh" 273#include "basdeclsP.fh" 274#include "geobasmapP.fh" 275#include "ecp_nwc.fh" 276c 277c::-passed 278 integer geom ! [input] geometry handle 279 integer basis ! [input] basis set handle 280 integer ecpis ! [input] ecp basis set handle 281 integer nbas(4) ! [input] atomic scf internal structure of 282* number of basis functions per sym 283* or shell type s, p, d, f 284 integer nz ! [input] number of exponents and coefs 285 double precision zeta(nz) ! [input] exponents 286 double precision cont(nz) ! [input] scaled atomic scf coefs 287 double precision ucont(nz)! [input] unscaled coefs (from bas) 288 integer sz_hatom ! [input] size of hatom array 289 double precision hatom(sz_hatom) ! [output] ecp part of V 290 integer icent ! [input] lexical index of center 291 integer sz_zetac ! [input] size of ecp exp/coef arrays 292* [s] == [scratch] used/filled locally and down tree but not 293* important above 294* [i] == [input] 295 double precision zeta_c(sz_zetac) ! [s] ecp exponent array 296 double precision coef_c(sz_zetac) ! [s] ecp coefs array 297 integer lmax_all ! [i] ecp max angular momentum 298 integer n_prim_c(0:4,-1:lmax_all) ! [s] ecp code pointer array 299 integer n_coef_c(-1:lmax_all) ! [s] ecp code pointer array 300 integer ind_c(-1:lmax_all) ! [s] ecp code pointer array 301 integer lscr ! [input] length of scratch buffer 302 double precision scr(lscr) ! [s] scratch buffer for ecp code 303 integer l_ecp_ints ! [input] size of ecp int buffer 304 double precision ecp_ints(l_ecp_ints) ! [s] ecp int buffer 305c::-local 306 double precision C(3) ! coordinates 307 integer type ! function type 308 integer ncoef ! number of coefs in contraction 309 integer nprim ! number of prims in contraction 310 integer ucent ! unique basis center 311 integer nn ! dummy contraction index 312 integer nn_off ! pointer offset for building ecp arrays 313 integer f_cont ! first contraction on center 314 integer l_cont ! last contraction on center 315 integer ecp ! lexical basis set index for ecpis 316 integer atscf_n_zetac ! size of n_zeta_c in reality 317 integer n0, n1, n2 ! dummy r exponent indecies 318 integer n3, n4 ! dummy r exponent indecies 319 integer iexp,icfp,irexp ! private basis set pointer info 320 integer l_c ! max ang on ecp center 321* for atomic code the next three are the same = 1 322 integer i_cent_C ! lexical atom index of ecp center 323 integer i_c_A ! lexical atom index of basis function center A 324 integer i_c_B ! lexical atom index of basis function center B 325 integer i, j, k, ! loop indices 326 & lj, lk, ljkoff ! pointer offsets 327 double precision 328 & expa, coefa, expb, coefb ! prim exp and coefs 329 integer sz_ints ! size of ints block 330 integer nbf_s ! nbf of shell in spherical d=5 331 integer nbf_x ! nbf of shell in cart. d=6 332 integer lscr_guess ! dummy arg to check memory in ecp code 333 integer cnt_hatom ! counter in hatom array 334c::-statement functions 335#include "bas_exndcf_sfn.fh" 336#include "bas_ibs_sfn.fh" 337c 338c.. get lexical basis set index for ecp basis set handle 339 ecp = ecpis + Basis_handle_offset 340 341c.. get unique center index. 342 ucent = sf_ibs_ce2uce(icent,ecp) 343*debug: write(6,*)' icent/ucent ',icent,'/',ucent 344c 345c.. get number of contractions/primitives on ecp tag 346 atscf_n_zetac = infbs_tags(Tag_Nprim,ucent,ecp) 347 348c.. check to make sure ecp arrays are big enough 349 if (atscf_n_zetac.gt.sz_zetac) then 350 write(luout,*)' nprimc/ncoefc array too small ' 351 write(luout,*)'atscf_n_zetac = ',atscf_n_zetac 352 write(luout,*)'sz_zetac = ',sz_zetac 353 write(luout,*)'contact nwchem-support@emsl.pnl.gov' 354 call errquit 355 & ('atscf_ecp_build_hatom: fatal error: ',911, BASIS_ERR) 356 endif 357 358c.. get first/last contraction on ecp tag 359 f_cont = infbs_tags(Tag_Fcont,ucent,ecp) 360 l_cont = infbs_tags(Tag_Lcont,ucent,ecp) 361*debug: write(6,*)' f/l cont',f_cont,'/',l_cont 362c 363c.. intialize some variables 364 l_c = -1000 365 nn_off = 1 366 367c.. loop over contractions 368 do nn = f_cont,l_cont 369c 370 type = infbs_cont(Cont_Type,nn,ecp) 371 nprim = infbs_cont(Cont_Nprim,nn,ecp) 372 ncoef = nprim*infbs_cont(Cont_Ngen,nn,ecp) 373 if (nprim.ne.ncoef) then 374 write(luout,*) 375 & 'general contraction ecp basis are invalid now' 376 call errquit('atscf_ecp_build_hatom: error',911, BASIS_ERR) 377 endif 378 iexp = infbs_cont(Cont_Iexp,nn,ecp) 379 icfp = infbs_cont(Cont_Icfp,nn,ecp) 380 irexp = infbs_cont(Cont_Irexp,nn,ecp) 381 if ((nn_off+nprim-1).gt.sz_zetac) call errquit 382 & ('atscf_ecp_build_hatom: too many exponents/coefficents',911, 383 & BASIS_ERR) 384*rak: call dcopy(nprim,dbl_mb(mb_exndcf(iexp,ecp)),1, 385*rak: & zeta_c(nn_off),1) 386*rak: call dcopy(nprim,dbl_mb(mb_exndcf(icfp,ecp)),1, 387*rak: & coef_c(nn_off),1) 388 call ecp_get_n3( 389 & zeta_c(nn_off), 390 & dbl_mb(mb_exndcf(iexp,ecp)), 391 & coef_c(nn_off), 392 & dbl_mb(mb_exndcf(icfp,ecp)), 393 & dbl_mb(mb_exndcf(irexp,ecp)),nprim,n0,n1,n2,n3,n4) 394 n_prim_c(0,type) = n0 395 n_prim_c(1,type) = n1 396 n_prim_c(2,type) = n2 397 n_prim_c(3,type) = n3 398 n_prim_c(4,type) = n4 399 ind_c(type) = nn_off 400 n_coef_c(type) = nprim 401 l_c = max(type,l_c) 402 nn_off = nn_off + nprim 403 enddo 404 i_cent_c = 1 405 i_c_A = 1 406 i_c_B = 1 407c 408*debug: write(6,*)' exponent list for ecp ' 409*debug: call output(zeta_c,1,atscf_n_zetac,1,1,atscf_n_zetac,1,1) 410*debug: write(6,*)' contraction list for ecp ' 411*debug: call output(coef_c,1,atscf_n_zetac,1,1,atscf_n_zetac,1,1) 412c 413c atom centered at origin 414 call dfill(3,0.0d00,C,1) 415*debug: write(6,*)'nbas',nbas 416*debug: write(6,*)'zeta',zeta 417*debug: write(6,*)'ucont',ucont 418 ljkoff = 0 419 cnt_hatom = 0 420c.. loop over shell types or atomic scf symmetries 421 do i = 1,4 422 type = i-1 ! type in NWChem terms 423 sz_ints = (type+1)*(type+2)/2 424 sz_ints = sz_ints*sz_ints ! size of integral block 425 lj = ljkoff 426 427c.. loop over primitives for <bra| 428 do j = 1,nbas(i) 429 lj = lj + 1 430 expa = zeta(lj) 431 coefa = ucont(lj) 432 lk = ljkoff 433c.. loop over primitives for |ket> 434 do k = 1,j 435 lk=lk+1 436 expb = zeta(lk) 437 coefb = ucont(lk) 438 lscr_guess = lscr 439*debug: write(6,*)' lscr ',lscr 440 441 call ecp_integral( 442 & C,expa,coefa,1,1,type,i_c_A, 443 & C,expb,coefb,1,1,type,i_c_B, 444 & C,zeta_c,coef_c,n_prim_C,n_coef_c, 445 & ind_c, ind_c, atscf_n_zetac, atscf_n_zetac, 446 & l_c,i_cent_C,1,lmax_all,0, 447 & dbl_mb(k_ecp_c2s),mem_c2s, 448 & ecp_ints,sz_ints,1, ! nblk = 1 for ecp integrals 449 & .true., 450 & scr,lscr_guess, 451 & 0) 452c.. make sure scratch array is big enough 453 if (lscr_guess.gt.lscr) then 454 write(luout,*)' lscr_guess =',lscr_guess 455 write(luout,*)' lscr =',lscr 456 write(luout,*)' contact nwchem-support@emsl.pnl.gov' 457 call errquit('atscf_ecp_build_hatom: fatal error',911, 458 & BASIS_ERR) 459 else 460c.. compute ecp integrals <bra|ecp|ket> 461 call ecp_integral( 462 & C,expa,coefa,1,1,type,i_c_A, 463 & C,expb,coefb,1,1,type,i_c_B, 464 & C,zeta_c,coef_c,n_prim_C,n_coef_c, 465 & ind_c, ind_c, atscf_n_zetac, atscf_n_zetac, 466 & l_c,i_cent_C,1,lmax_all,0, 467 & dbl_mb(k_ecp_c2s),mem_c2s, 468 & ecp_ints,sz_ints,1, ! nblk = 1 for ecp integrals 469 & .false., 470 & scr,lscr, 471 & 0) 472* write(6,*)' ecp integrals, cart' 473* call output(ecp_ints,1,sz_ints,1,1,sz_ints,1,1) 474 nbf_x = (type+1)*(type+2)/2 475 nbf_s = 2*type+1 476 477c.. transform cartesian block to spherical since the atomic scf 478* is a spherical code. This makes the integral block 479* diagonal with the same values for all 2*l+1 components. 480 481 call spcart_tran1e(ecp_ints,scr, 482 & nbf_x,nbf_x,type,1, 483 & nbf_s,nbf_s,type,1, 484 & .false.) 485* write(6,*)' ecp integrals, spherical' 486* call output(ecp_ints,1,(nbf_s*nbf_s), 487* & 1,1,(nbf_s*nbf_s),1,1) 488 endif 489 490c.. compute index into hatom array 491 cnt_hatom = cnt_hatom + 1 492c... use the zero component of the (-l:l) square spherical integral 493* block. nbf_x is now the index into odd rank square matrix 494* (e.g., find the center element) 495 nbf_x = (nbf_s*nbf_s - 1)/2 + 1 496 hatom(cnt_hatom) = ecp_ints(nbf_x) 497 enddo 498 enddo 499 500c.. compute offset into scalar atomic scf zeta and ucont arrays 501 ljkoff = ljkoff + nbas(i) 502 enddo 503 end 504