1 logical function int_ecp_init(ecpidin,soidin,indx_grad) 2* $Id$ 3 implicit none 4*::cr::7 5*--------------------------------------------------* 6* COPYRIGHT (C) 1994, 1995, 1996, 1997, 1998, 1999 * 7* Pacific Northwest National Laboratory, * 8* Battelle Memorial Institute. * 9*--------------------------------------------------* 10*------------> All Rights Reserved <---------------* 11*--------------------------------------------------* 12* 13#include "errquit.fh" 14#include "mafdecls.fh" 15#include "bas.fh" 16#include "nwc_const.fh" 17#include "basP.fh" 18#include "basdeclsP.fh" 19#include "geobasmapP.fh" 20#include "bas_exndcf_dec.fh" 21#include "bas_ibs_dec.fh" 22#include "geom.fh" 23#include "geomP.fh" 24#include "apiP.fh" 25#include "ecp_nwc.fh" 26#include "int_nbf.fh" 27#include "stdio.fh" 28#include "sym.fh" 29*::passed 30 integer ecpidin ! [input] ecp basis handle 31 integer soidin ! [input] so potential handle 32 integer indx_grad ! [input] gradient level 0=energy 1=grad 2=hess 33*::local 34c 35 integer ecpid ! lexical index for ecp basis handle 36 integer soid ! lexical index for ecp basis handle 37 integer geom ! geometry handle 38 integer geomecp ! geometry handle from ecp 39 integer geomso ! geometry handle from so pot. 40* integer num_ecp ! number of ecp centers (counted) 41 integer ncenters ! number of centers 42 integer n_zeta_c_e ! length of ecp exponent array 43 integer n_zeta_c_c ! length of ecp coefficient array 44 integer nz_add ! increment for exp/coef counters 45* integer icent ! counter for center loop 46* integer iucent ! unique center of icent 47 integer l_ecp_sz ! size of ang info for ecp pointer arrays 48 integer l_so ! maximum angular momentum of so pot. 49 integer lmax_both 50 integer idum 51 integer basis_handle 52 integer l_bas 53 integer bas_id 54 integer basis_ncont 55 integer i_cont, ic, iuc, ie, icf, inp, ing, il, ig 56 integer j_cont, jc, juc, je, jcf, jnp, jng, jl, jg 57 integer lecp_mem, mem_ecp, lecp_mem1 58* integer h_tlist, k_tlist 59 integer nV 60 logical o_both, o_ecp, o_so 61 double precision dpdum(2),qfact 62 logical oskel 63 integer ictr,ic1,ic2 64 integer jctr,jc1,jc2 65* 66 external ecp_init_bd 67* 68#include "bas_exndcf_sfn.fh" 69#include "bas_ibs_sfn.fh" 70*<<<------------------------------------------------------------> 71*<<<-First pass is to explode all information into arrays 72*<<<-The memory used can be reduced by keeping track of the 73*<<< tag to basis-unique tag information so that the pointers 74*<<< for exponent and coefficient arrays is similar to the 75*<<< storage in the basis set object. It may be possible to 76*<<< pass the entire exndcf array for the ecp basis to the 77*<<< integral routines with the proper pointer mechanism in 78*<<< place. 79*<<< 80*<<< RAK Apr 1996 81*<<<------------------------------------------------------------> 82c 83c 84c.. check initialization 85 if (init_ecp_init) then 86 write(6,*)' already called int_ecp_init' 87 call errquit('int_ecp_init error',911, INT_ERR) 88 endif 89 init_ecp_init = .true. 90c.. determine if ecp/so or both are needed 91 o_ecp = ecpidin.ne.0 92 o_so = soidin.ne.0 93 o_both = o_ecp.and.o_so 94c.. lexical basis index 95 ecpid = -1 96 soid = -1 97 if (o_ecp) ecpid = ecpidin + BASIS_HANDLE_OFFSET 98 if (o_so) soid = soidin + BASIS_HANDLE_OFFSET 99 100c.. determine geom 101 if (o_ecp) then 102 geomecp = ibs_geom(ecpid) 103 geom = geomecp 104 endif 105 if (o_so) then 106 geomso = ibs_geom(soid) 107 geom = geomso 108 endif 109 if (o_both) then 110 if (geomso.ne.geomecp) then 111 write(luout,*)' geometry used to load ecp ',geomecp 112 write(luout,*)' geometry used to load so ',geomso 113 call errquit('int_ecp_init: ecp/so geometry mismatch',911, 114 & GEOM_ERR) 115 endif 116 geom = geomecp 117 endif 118* 119* build up ecp/so information for heap 120c.... 121c.. get total number of centers 122 if (.not.geom_ncent(geom,ncenters)) call errquit 123 & ('int_ecp_init: geom_ncent failed?',911, GEOM_ERR) 124c.. get number of ecp/so/both centers 125 n_ecp = 0 126 if (o_both) then 127 if (.not.ecpso_ncent 128 & (geom,soidin,ecpidin,n_ecp)) call errquit 129 & ('int_ecp_init: ecpso_ncent failed',911, INT_ERR) 130 else if (o_ecp) then 131 if (.not.ecp_ncent 132 & (geom,ecpidin,n_ecp)) call errquit 133 & ('int_ecp_init: ecp_ncent failed?',911, INT_ERR) 134 else if (o_so) then 135 if (.not.so_ncent 136 & (geom,soidin,n_ecp)) call errquit 137 & ('int_ecp_init: so_ncent failed',911, INT_ERR) 138 else 139 call errquit('int_ecp_init: no so or ecp basis',911, 140 & BASIS_ERR) 141 endif 142 if (n_ecp.eq.0) call errquit 143 & ('int_ecp_init: fatal error no ecp/so centers',911, 144 & BASIS_ERR) 145 146c.. allocate space for coordintates of ecp/so centers 147 if (.not.ma_alloc_get(mt_dbl, 148 & (3*n_ecp), 149 & 'ecp center coords', 150 & h_xyzecp, k_xyzecp)) call errquit 151 & (' int_ecp_init: ecp center coords ma failed ',911, 152 & BASIS_ERR) 153 call dfill((3*n_ecp),0.0d00,dbl_mb(k_xyzecp),1) 154 155c.. get coordinates for ecp centers. 156 if (o_both) then 157 if (.not.ecpso_get_coords(geom,soidin,ecpidin, 158 & n_ecp,dbl_mb(k_xyzecp))) 159 & call errquit('int_ecp_init: ecpso_get_coords failed',911, 160 & INT_ERR) 161 else if (o_ecp) then 162 if (.not.ecp_get_coords(geom,ecpidin, 163 & n_ecp,dbl_mb(k_xyzecp))) 164 & call errquit('int_ecp_init: ecp_get_coords failed',911, 165 & INT_ERR) 166 else if (o_so) then 167 if (.not.so_get_coords(geom,soidin, 168 & n_ecp,dbl_mb(k_xyzecp))) 169 & call errquit('int_ecp_init: so_get_coords failed',911, 170 & INT_ERR) 171 else 172 call errquit('int_ecp_init: no so or ecp basis',911, 173 & INT_ERR) 174 endif 175* write(6,*)' coordinates after read ' 176* call output(dbl_mb(k_xyzecp),1,3,1,n_ecp,3,n_ecp,1) 177c 178c... now comes the tricky part 179c 180c allocate and fill an exponent and coeff array for the ecp basis 181c 182* write(6,*)'inside ecp init' 183* if (.not.ecp_print(ecpidin)) stop ' ecp_print error' 184 if (o_both) then 185 if (.not.bas_num_prim(ecpidin,nz_add)) call errquit 186 & ('int_ecp_init:bas_num_prim failed?',911, BASIS_ERR) 187 n_zeta_c_e = nz_add 188 if (.not.bas_num_prim(soidin,nz_add)) call errquit 189 & ('int_ecp_init:bas_num_prim failed?',911, 190 & BASIS_ERR) 191 n_zeta_c_e = n_zeta_c_e + nz_add 192 if (.not.bas_num_coef(ecpidin,nz_add)) call errquit 193 & ('int_ecp_init:bas_num_coef failed?',911, BASIS_ERR) 194 n_zeta_c_c = nz_add 195 if (.not.bas_num_coef(soidin,nz_add)) call errquit 196 & ('int_ecp_init:bas_num_prim failed?',911, BASIS_ERR) 197 n_zeta_c_c = n_zeta_c_c + nz_add 198 else if (o_ecp) then 199 if (.not.bas_num_prim(ecpidin,nz_add)) call errquit 200 & ('int_ecp_init:bas_num_prim failed?',911, BASIS_ERR) 201 n_zeta_c_e = nz_add 202 if (.not.bas_num_coef(ecpidin,nz_add)) call errquit 203 & ('int_ecp_init:bas_num_coef failed?',911, BASIS_ERR) 204 n_zeta_c_c = nz_add 205 else if (o_so) then 206 if (.not.bas_num_prim(soidin,nz_add)) call errquit 207 & ('int_ecp_init:bas_num_prim failed?',911, BASIS_ERR) 208 n_zeta_c_e = nz_add 209 if (.not.bas_num_coef(soidin,nz_add)) call errquit 210 & ('int_ecp_init:bas_num_coef failed?',911, BASIS_ERR) 211 n_zeta_c_c = nz_add 212 else 213 call errquit('int_ecp_init: no so or ecp basis',911, BASIS_ERR) 214 endif 215*.. error check for general contraction 216 if (n_zeta_c_e.ne.n_zeta_c_c) then 217 write(6,*)' possible general contraction on ecp/so basis' 218 call errquit ('int_ecp_init: n_zeta_c_e .ne. n_zeta_c_c', 219 & (n_zeta_c_e- n_zeta_c_c), BASIS_ERR) 220 endif 221 n_zeta_c = n_zeta_c_e ! use e one 222c 223cc AJL/Begin/SPIN ECPs 224c 225c... determine the number of channels needed (Spin-polarised ECPs) 226 ecp_channels = 1 227 if (o_ecp) then 228 if (.not.ecp_get_high_chan(ecpidin,ecp_channels)) 229 & call errquit('int_ecp_init error',911, INT_ERR) 230 end if 231cc AJL/Debug 232c write(6,*) 'Number of ECP Channels: ', ecp_channels 233c write(6,*) '(1 = Spin Independent, Default; 2 = Spin Polarised)' 234cc AJL/End 235c 236c.. allocate space for ecp/so exponents and ecp/so coefficients 237 if (.not.ma_alloc_get(mt_dbl, 238 & n_zeta_c, 239 & 'ecp exponents', 240 & h_ecp_e, k_ecp_e)) call errquit 241 & (' int_ecp_init: ecp exponent ma failed ',911, BASIS_ERR) 242 call dfill(n_zeta_c,0.0d00,dbl_mb(k_ecp_e),1) 243 if (.not.ma_alloc_get(mt_dbl, 244 & n_zeta_c, 245 & 'ecp coefficients', 246 & h_ecp_c, k_ecp_c)) call errquit 247 & (' int_ecp_init: ecp coefficients ma failed ',911, BASIS_ERR) 248 call dfill(n_zeta_c,0.0d00,dbl_mb(k_ecp_c),1) 249c 250cc AJL/Begin/SPIN ECPs 251 if (ecp_channels.gt.1) then 252 if (.not.ma_alloc_get(mt_dbl, 253 & n_zeta_c, 254 & 'beta ecp exponents', 255 & h_ecp_e_beta, k_ecp_e_beta)) call errquit 256 & (' int_ecp_init: beta ecp exponent ma failed ', 257 & 911, BASIS_ERR) 258 call dfill(n_zeta_c,0.0d00,dbl_mb(k_ecp_e_beta),1) 259 if (.not.ma_alloc_get(mt_dbl, 260 & n_zeta_c, 261 & 'beta ecp coefficients', 262 & h_ecp_c_beta, k_ecp_c_beta)) call errquit 263 & (' int_ecp_init: beta ecp coefficients ma failed ', 264 & 911, BASIS_ERR) 265 call dfill(n_zeta_c,0.0d00,dbl_mb(k_ecp_c_beta),1) 266 end if 267cc AJL/End 268c 269c... determine maximum angular momentum of ecp basis 270 l_ecp = 0 271 if (o_ecp) then 272 if (.not.bas_high_angular(ecpidin,l_ecp)) call errquit 273 & ('int_ecp_init:ecp: bas_high_angular failed',911, BASIS_ERR) 274 endif 275 if (o_so) then 276 if (.not.bas_high_angular(soidin,l_so)) call errquit 277 & ('int_ecp_init:so: bas_high_angular failed',911, BASIS_ERR) 278 l_ecp = max(l_ecp,l_so) 279 endif 280c 281 l_ecp_sz = l_ecp + 2 ! (-1->Lval) 282c 283c... allocate space for n_prim_C(0:4,-1:l_ecp_max,n_C*2), 284 if (.not.ma_alloc_get(mt_int, 285 & (5*l_ecp_sz*n_ecp*2), 286 & 'ecp n_prim_C', 287 & h_ecp_nprim_c, k_ecp_nprim_c)) call errquit 288 & (' int_ecp: ecp nprim_c ma failed ',911, BASIS_ERR) 289 call ifill((5*l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_nprim_c),1) 290c... allocate space for n_coef_C(-1:l_ecp_max,n_C*2) 291 if (.not.ma_alloc_get(mt_int, 292 & (l_ecp_sz*n_ecp*2), 293 & 'ecp n_coef_C', 294 & h_ecp_ncoef_c, k_ecp_ncoef_c)) call errquit 295 & (' int_ecp: ecp ncoef_c ma failed ',911, BASIS_ERR) 296 call ifill((l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_ncoef_c),1) 297c... allocate space for ind_C 298 if (.not.ma_alloc_get(mt_int, 299 & (l_ecp_sz*n_ecp*2), 300 & 'ecp ind_C', 301 & h_ecp_ind_c, k_ecp_ind_c)) call errquit 302 & (' int_ecp: ecp ind_c ma failed ',911, BASIS_ERR) 303 call ifill((l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_ind_c),1) 304c... allocate space for ind_z 305 if (.not.ma_alloc_get(mt_int, 306 & (l_ecp_sz*n_ecp*2), 307 & 'ecp ind_z', 308 & h_ecp_ind_z, k_ecp_ind_z)) call errquit 309 & (' int_ecp: ecp ind_z ma failed ',911, BASIS_ERR) 310 call ifill((l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_ind_z),1) 311c... allocate space for l_C 312 if (.not.ma_alloc_get(mt_int, 313 & n_ecp, 314 & 'ecp l_C', 315 & h_ecp_l_c, k_ecp_l_c)) call errquit 316 & (' int_ecp: ecp l_c ma failed ',911, BASIS_ERR) 317 call ifill(n_ecp,0,int_mb(k_ecp_l_c),1) 318c... allocate space for ecp center pointer list 319 if (.not.ma_alloc_get(mt_int, 320 & n_ecp, 321 & 'ecp lexical indeces for ecp centers', 322 & h_ecp_lip, k_ecp_lip)) call errquit 323 & (' int_ecp: ecp_lip ma failed',911, BASIS_ERR) 324 call ifill(n_ecp,0,int_mb(k_ecp_lip),1) 325c 326cc AJL/Begin/SPIN ECPs 327c... allocate space for n_prim_C(0:4,-1:l_ecp_max,n_C*2), 328 if (ecp_channels.gt.1) then 329 if (.not.ma_alloc_get(mt_int, 330 & (5*l_ecp_sz*n_ecp*2), 331 & 'beta ecp n_prim_C', 332 & h_ecp_nprim_c_beta, k_ecp_nprim_c_beta)) call errquit 333 & (' int_ecp: beta ecp nprim_c ma failed ',911, BASIS_ERR) 334 call ifill((5*l_ecp_sz*n_ecp*2), 335 & 0,int_mb(k_ecp_nprim_c_beta),1) 336c... allocate space for n_coef_C(-1:l_ecp_max,n_C*2) 337 if (.not.ma_alloc_get(mt_int, 338 & (l_ecp_sz*n_ecp*2), 339 & 'beta ecp n_coef_C', 340 & h_ecp_ncoef_c_beta, k_ecp_ncoef_c_beta)) call errquit 341 & (' int_ecp: beta ecp ncoef_c ma failed ',911, BASIS_ERR) 342 call ifill((l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_ncoef_c_beta),1) 343c... allocate space for ind_C 344 if (.not.ma_alloc_get(mt_int, 345 & (l_ecp_sz*n_ecp*2), 346 & 'beta ecp ind_C', 347 & h_ecp_ind_c_beta, k_ecp_ind_c_beta)) call errquit 348 & (' int_ecp: beta ecp ind_c ma failed ',911, BASIS_ERR) 349 call ifill((l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_ind_c_beta),1) 350c... allocate space for ind_z 351 if (.not.ma_alloc_get(mt_int, 352 & (l_ecp_sz*n_ecp*2), 353 & 'beta ecp ind_z', 354 & h_ecp_ind_z_beta, k_ecp_ind_z_beta)) call errquit 355 & (' int_ecp: beta ecp ind_z ma failed ',911, BASIS_ERR) 356 call ifill((l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_ind_z_beta),1) 357c... allocate space for l_C 358 if (.not.ma_alloc_get(mt_int, 359 & n_ecp, 360 & 'beta ecp l_C', 361 & h_ecp_l_c_beta, k_ecp_l_c_beta)) call errquit 362 & (' int_ecp: beta ecp l_c ma failed ',911, BASIS_ERR) 363 call ifill(n_ecp,0,int_mb(k_ecp_l_c_beta),1) 364c... allocate space for ecp center pointer list 365 if (.not.ma_alloc_get(mt_int, 366 & n_ecp, 367 & 'beta ecp lexical indices for ecp centers', 368 & h_ecp_lip_beta, k_ecp_lip_beta)) call errquit 369 & (' int_ecp: beta ecp_lip ma failed',911, BASIS_ERR) 370 call ifill(n_ecp,0,int_mb(k_ecp_lip_beta),1) 371 endif 372cc AJL/End 373c 374c 375 call int_ecp_build_ecp_ptrs(ecpidin,soidin,geom, 376 & o_both,o_ecp,o_so, 377 & ncenters, 378 & n_ecp, 379 & l_ecp, 380 & n_zeta_c, 381cc AJL/Begin/SPIN ECPs 382 & 1, ! Channel to be filled. Would be 2 for beta channel 383 & ecp_channels, ! Total no. of channels, for print_ptrs 384cc AJL/End 385 & int_mb(k_ecp_nprim_c), 386 & int_mb(k_ecp_ncoef_c), 387 & int_mb(k_ecp_ind_c), 388 & int_mb(k_ecp_ind_z), 389 & int_mb(k_ecp_l_c), 390 & int_mb(k_ecp_lip), 391 & dbl_mb(k_ecp_e), 392 & dbl_mb(k_ecp_c) ) 393c 394cc AJL/Begin/SPIN ECPs 395 if (ecp_channels.gt.1) ! Setup beta options 396 & call int_ecp_build_ecp_ptrs(ecpidin,soidin,geom, 397 & o_both,o_ecp,o_so, 398 & ncenters, 399 & n_ecp, 400 & l_ecp, 401 & n_zeta_c, 402 & 2, ! Beta channel 403 & ecp_channels, ! Total no. of channels, for print_ptrs 404 & int_mb(k_ecp_nprim_c_beta), 405 & int_mb(k_ecp_ncoef_c_beta), 406 & int_mb(k_ecp_ind_c_beta), 407 & int_mb(k_ecp_ind_z_beta), 408 & int_mb(k_ecp_l_c_beta), 409 & int_mb(k_ecp_lip_beta), 410 & dbl_mb(k_ecp_e_beta), 411 & dbl_mb(k_ecp_c_beta)) 412cc AJL/End 413c 414*... allocate space for c2s and s2c internal ecp transformation routines 415c 416c determine lmax among ao basis and ecp basis 417c l_ecp currently has Lmax for ecp basis 418 if (.not.ecp_get_parent_handle(ecpidin,basis_handle)) 419 & call errquit 420 & ('int_ecp_init: ecp_get_parent_handle failed',911, 421 & BASIS_ERR) 422 if (.not.bas_high_angular(basis_handle,l_bas)) call errquit 423 & ('int_ecp_init: bas_high_angular failed for ao handle', 424 & 911, BASIS_ERR) 425 lmax_both = max(l_ecp,l_bas) + l_bas + 2 426 call ecp_init_c2s(lmax_both,dpdum,dpdum,idum,1,1,.true.,mem_c2s) 427 if (.not.ma_alloc_get(mt_dbl, 428 & mem_c2s, 429 & 'ecp c2s routines', 430 & h_ecp_c2s, k_ecp_c2s)) call errquit 431 & ('int_ecp_init: ma failed for c2s',911, MA_ERR) 432 call dfill(mem_c2s,0.0d00,dbl_mb(k_ecp_c2s),1) 433 call ecp_init_c2s(lmax_both, 434 & dbl_mb(k_ecp_c2s),dbl_mb(k_ecp_c2s),mem_c2s,1,1,.false., 435 & idum) 436c 437c initialize constants for ecp integral code 438c 439 call ecp_init_con() 440c 441c determine maximum memory for ecp integrals 442c 443 bas_id = basis_handle + basis_handle_offset 444 basis_ncont = ncont_tot_gb(bas_id) 445 ig = ibs_geom(bas_id) 446 jg = ig 447 mem_ecp = 0 448* if (.not.bas_print(basis_handle)) stop ' error' 449* if (.not.gbs_map_print(basis_handle)) stop ' error ' 450* if (.not.bas_print(ecpidin)) stop 'error' 451* if (.not.gbs_map_print(ecpidin)) stop ' error' 452 oskel=.true. 453 do ictr=1,ncenters 454 if (sym_atom(ig, ictr, qfact)) then 455 if(.not.(bas_ce2cnr(basis_handle,ictr,ic1,ic2))) 456 . call errquit('intecp:basce2cnr failed',0, BASIS_ERR) 457 do i_cont = ic1,ic2 458 ic = sf_ibs_cn2ce(i_cont,bas_id) 459 iuc = sf_ibs_cn2ucn(i_cont,bas_id) 460 ie = infbs_cont(cont_iexp,iuc,bas_id) 461 icf = infbs_cont(cont_icfp,iuc,bas_id) 462 inp = infbs_cont(cont_nprim,iuc,bas_id) 463 ing = infbs_cont(cont_ngen,iuc,bas_id) 464 il = infbs_cont(cont_type,iuc,bas_id) 465 do jctr=ictr,ncenters 466 if (sym_atom(ig, jctr, qfact)) then 467 if(.not.(bas_ce2cnr(basis_handle,jctr,jc1,jc2))) 468 . call errquit('intecp:basce2cnr failed',0, 469 & BASIS_ERR) 470 do j_cont = jc1,jc2 471 jc = sf_ibs_cn2ce(j_cont,bas_id) 472 juc = sf_ibs_cn2ucn(j_cont,bas_id) 473 je = infbs_cont(cont_iexp,juc,bas_id) 474 jcf = infbs_cont(cont_icfp,juc,bas_id) 475 jnp = infbs_cont(cont_nprim,juc,bas_id) 476 jng = infbs_cont(cont_ngen,juc,bas_id) 477 jl = infbs_cont(cont_type,juc,bas_id) 478 nV = int_nbf_x(il)*int_nbf_x(jl)*ing*jng 479 lecp_mem = 90000 480 if (indx_grad.eq.0) then 481 if (o_ecp) then 482 call int_ecp_hf1( 483 & coords(1,ic,ig), 484 & dbl_mb(mb_exndcf(ie,bas_id)), 485 & dbl_mb(mb_exndcf(icf,bas_id)), 486 & inp,ing,il,ic, 487 & coords(1,jc,jg), 488 & dbl_mb(mb_exndcf(je,bas_id)), 489 & dbl_mb(mb_exndcf(jcf,bas_id)), 490 & jnp,jng,jl,jc, 491 & dpdum,nV,dpdum,lecp_mem1,.true.) 492 lecp_mem = max(lecp_mem,lecp_mem1) 493 endif 494 if (o_so) then 495 call intso_hf1( 496 & coords(1,ic,ig), 497 & dbl_mb(mb_exndcf(ie,bas_id)), 498 & dbl_mb(mb_exndcf(icf,bas_id)), 499 & inp,ing,il,ic, 500 & coords(1,jc,jg), 501 & dbl_mb(mb_exndcf(je,bas_id)), 502 & dbl_mb(mb_exndcf(jcf,bas_id)), 503 & jnp,jng,jl,jc, 504 & dpdum,nV,dpdum,lecp_mem1,.true.) 505 lecp_mem = max(lecp_mem,lecp_mem1) 506 endif 507 elseif (indx_grad.eq.1) then 508 if (o_ecp) then 509 call intd_ecp_hf1( 510 & coords(1,ic,ig), 511 & dbl_mb(mb_exndcf(ie,bas_id)), 512 & dbl_mb(mb_exndcf(icf,bas_id)), 513 & inp,ing,il,ic, 514 & coords(1,jc,jg), 515 & dbl_mb(mb_exndcf(je,bas_id)), 516 & dbl_mb(mb_exndcf(jcf,bas_id)), 517 & jnp,jng,jl,jc, 518 & dpdum,nV,ncenters,dpdum,lecp_mem1, 519 & .true.) 520 lecp_mem = max(lecp_mem,lecp_mem1) 521 endif 522 if (o_so) then 523 call intd_so_hf1( 524 & coords(1,ic,ig), 525 & dbl_mb(mb_exndcf(ie,bas_id)), 526 & dbl_mb(mb_exndcf(icf,bas_id)), 527 & inp,ing,il,ic, 528 & coords(1,jc,jg), 529 & dbl_mb(mb_exndcf(je,bas_id)), 530 & dbl_mb(mb_exndcf(jcf,bas_id)), 531 & jnp,jng,jl,jc, 532 & dpdum,nV,ncenters,dpdum,lecp_mem1, 533 & .true.) 534 lecp_mem = max(lecp_mem,lecp_mem1) 535 endif 536 elseif (indx_grad.eq.2) then 537 if (o_ecp) then 538 call intdd_ecp_hf1( 539 & coords(1,ic,ig), 540 & dbl_mb(mb_exndcf(ie,bas_id)), 541 & dbl_mb(mb_exndcf(icf,bas_id)), 542 & inp,ing,il,ic, 543 & coords(1,jc,jg), 544 & dbl_mb(mb_exndcf(je,bas_id)), 545 & dbl_mb(mb_exndcf(jcf,bas_id)), 546 & jnp,jng,jl,jc, 547 & dpdum,nV,ncenters,dpdum,lecp_mem1, 548 & .true.) 549 lecp_mem = max(lecp_mem,lecp_mem1) 550 endif 551 if (o_so) then 552 call intdd_so_hf1( 553 & coords(1,ic,ig), 554 & dbl_mb(mb_exndcf(ie,bas_id)), 555 & dbl_mb(mb_exndcf(icf,bas_id)), 556 & inp,ing,il,ic, 557 & coords(1,jc,jg), 558 & dbl_mb(mb_exndcf(je,bas_id)), 559 & dbl_mb(mb_exndcf(jcf,bas_id)), 560 & jnp,jng,jl,jc, 561 & dpdum,nV,ncenters,dpdum,lecp_mem1, 562 & .true.) 563 lecp_mem = max(lecp_mem,lecp_mem1) 564 endif 565 else 566 call errquit( 567 . 'int_ecp_init:unknown initializ',911, 568 & INT_ERR) 569 endif 570c ! nint block used in api routine int_hf1sp 571 lecp_mem = lecp_mem + nV 572c 573cc AJL/Begin/SPIN ECPs 574cc Need to account for scenario where we have two channels so 575cc double the number of ECPs could be generated 576 if (o_ecp.and.ecp_channels.gt.1) then 577 if (indx_grad.eq.0) then 578 call int_ecp_hf1_beta( 579 & coords(1,ic,ig), 580 & dbl_mb(mb_exndcf(ie,bas_id)), 581 & dbl_mb(mb_exndcf(icf,bas_id)), 582 & inp,ing,il,ic, 583 & coords(1,jc,jg), 584 & dbl_mb(mb_exndcf(je,bas_id)), 585 & dbl_mb(mb_exndcf(jcf,bas_id)), 586 & jnp,jng,jl,jc, 587 & dpdum,nV,dpdum,lecp_mem,.true.) 588 elseif (indx_grad.eq.1) then 589 call intd_ecp_hf1_beta( 590 & coords(1,ic,ig), 591 & dbl_mb(mb_exndcf(ie,bas_id)), 592 & dbl_mb(mb_exndcf(icf,bas_id)), 593 & inp,ing,il,ic, 594 & coords(1,jc,jg), 595 & dbl_mb(mb_exndcf(je,bas_id)), 596 & dbl_mb(mb_exndcf(jcf,bas_id)), 597 & jnp,jng,jl,jc, 598 & dpdum,nV,ncenters,dpdum,lecp_mem,.true.) 599 elseif (indx_grad.eq.2) then 600 call intdd_ecp_hf1_beta( 601 & coords(1,ic,ig), 602 & dbl_mb(mb_exndcf(ie,bas_id)), 603 & dbl_mb(mb_exndcf(icf,bas_id)), 604 & inp,ing,il,ic, 605 & coords(1,jc,jg), 606 & dbl_mb(mb_exndcf(je,bas_id)), 607 & dbl_mb(mb_exndcf(jcf,bas_id)), 608 & jnp,jng,jl,jc, 609 & dpdum,nV,ncenters,dpdum,lecp_mem,.true.) 610 endif 611c Add data storage space to lecp_mem 612 lecp_mem = lecp_mem + nV 613 endif 614cc AJL/End 615c 616* write(6,*)' i_cont ',i_cont 617* write(6,*)' j_cont ',j_cont 618 mem_ecp = max(mem_ecp,lecp_mem) 619* write(6,*)' lecp_mem from int_ecp_init is ', 620* & lecp_mem,mem_ecp 621 enddo 622 endif 623 enddo 624 enddo 625 endif 626 enddo 627* write(6,*)' scr 1e memory without ecp',mem_1e 628 mem_1e = max(mem_1e,mem_ecp) 629* write(6,*)' scr 1e memory with ecp',mem_1e 630* 631 call int_nbf_max(basis_handle,mem_ecp) 632c 633 if (indx_grad.eq.0) then 634 mem_ecp=mem_ecp*mem_ecp 635 elseif (indx_grad.eq.1) then 636 mem_ecp=mem_ecp*mem_ecp*3*ncenters 637 else 638 mem_ecp=mem_ecp*mem_ecp*3*3*(ncenters*(ncenters-1)/2+ncenters) 639 mem_1e = max(mem_1e,mem_ecp) 640 endif 641 isz_1e = max (isz_1e,mem_ecp) 642c 643 int_ecp_init = .true. 644* 645 end 646 subroutine int_ecp_terminate() 647 implicit none 648#include "errquit.fh" 649#include "mafdecls.fh" 650#include "ecp_nwc.fh" 651 652 logical status 653 654 status = .true. 655 status = status .and. MA_free_heap(h_xyzecp) 656 status = status .and. MA_free_heap(h_ecp_e) 657 status = status .and. MA_free_heap(h_ecp_c) 658 status = status .and. MA_free_heap(h_ecp_nprim_c) 659 status = status .and. MA_free_heap(h_ecp_ncoef_c) 660 status = status .and. MA_free_heap(h_ecp_ind_c) 661 status = status .and. MA_free_heap(h_ecp_ind_z) 662 status = status .and. MA_free_heap(h_ecp_l_c) 663 status = status .and. MA_free_heap(h_ecp_c2s) 664 status = status .and. MA_free_heap(h_ecp_lip) 665 h_xyzecp = 0 666 k_xyzecp = 0 667 h_ecp_e = 0 668 k_ecp_e = 0 669 h_ecp_c = 0 670 k_ecp_c = 0 671 h_ecp_nprim_c = 0 672 k_ecp_nprim_c = 0 673 h_ecp_ncoef_c = 0 674 k_ecp_ncoef_c = 0 675 h_ecp_ind_c = 0 676 k_ecp_ind_c = 0 677 h_ecp_ind_z = 0 678 k_ecp_ind_z = 0 679 h_ecp_l_c = 0 680 k_ecp_l_c = 0 681 h_ecp_c2s = 0 682 k_ecp_c2s = 0 683 h_ecp_lip = 0 684 k_ecp_lip = 0 685 n_zeta_c = 0 686 l_ecp = 0 687 n_ecp = 0 688 init_ecp_init = .false. 689c 690cc AJL/Begin 691 if (ecp_channels.gt.1) then 692 status = status .and. MA_free_heap(h_ecp_e_beta) 693 status = status .and. MA_free_heap(h_ecp_c_beta) 694 status = status .and. MA_free_heap(h_ecp_nprim_c_beta) 695 status = status .and. MA_free_heap(h_ecp_ncoef_c_beta) 696 status = status .and. MA_free_heap(h_ecp_ind_c_beta) 697 status = status .and. MA_free_heap(h_ecp_ind_z_beta) 698 status = status .and. MA_free_heap(h_ecp_l_c_beta) 699 status = status .and. MA_free_heap(h_ecp_lip_beta) 700 701 h_ecp_e_beta = 0 702 k_ecp_e_beta = 0 703 h_ecp_c_beta = 0 704 k_ecp_c_beta = 0 705 h_ecp_nprim_c_beta = 0 706 k_ecp_nprim_c_beta = 0 707 h_ecp_ncoef_c_beta = 0 708 k_ecp_ncoef_c_beta = 0 709 h_ecp_ind_c_beta = 0 710 k_ecp_ind_c_beta = 0 711 h_ecp_ind_z_beta = 0 712 k_ecp_ind_z_beta = 0 713 h_ecp_l_c_beta = 0 714 k_ecp_l_c_beta = 0 715 h_ecp_lip_beta = 0 716 k_ecp_lip_beta = 0 717 ecp_channels = 0 ! Clear out ECP channels and we are done 718 end if 719cc AJL/End 720c 721 if (status) return 722 call errquit 723 & ('int_ecp_terminate: error freeing heap',911, MEM_ERR) 724 end 725 subroutine int_ecp_build_ecp_ptrs(ecpidin,soidin,geom, 726 & o_both, o_ecp, o_so, 727 & ncenters, 728 & n_ecp, 729 & l_ecp, 730 & nz_ecp, 731cc AJL/Begin/SPIN ECPs 732 & channel, ! spin channel 733 & channels, ! total number of channels 734cc AJL/End 735 & n_prim_C, 736 & n_coef_C, 737 & ind_C, 738 & ind_Z, 739 & l_C, 740 & i_cent_C, 741 & c_exp, 742 & c_coef) 743 implicit none 744#include "errquit.fh" 745#include "bas.fh" 746#include "nwc_const.fh" 747#include "basP.fh" 748#include "basdeclsP.fh" 749#include "geobasmapP.fh" 750#include "geomP.fh" 751#include "mafdecls.fh" 752#include "bas_exndcf_dec.fh" 753#include "bas_ibs_dec.fh" 754c 755 integer ecpidin ! [input] ecp basis set handle 756 integer soidin ! [input] so pot. handle 757 logical o_both, o_ecp, o_so ![input] which ones to set up 758 integer ncenters ! [input] number of centers 759 integer n_ecp ! [input] number of ecp centers 760 integer l_ecp ! [input] maximum angular momentum in ecp basis 761 integer nz_ecp ! [input] number of prims/coeffs in stored data structure. 762cc AJL/Begin/SPIN ECPs 763 integer channel ! [input] Current spin channel used for ECPs 764 integer channels ! [input] Total channels, used for print 765cc AJL/End 766 integer n_prim_C(0:4,-1:l_ecp,n_ecp,2) ! [output] 767 integer n_coef_C(-1:l_ecp,n_ecp,2) ! [output] 768 integer ind_C(-1:l_ecp,n_ecp,2) ! [output] 769 integer ind_z(-1:l_ecp,n_ecp,2) ! [output] 770 integer l_C(n_ecp) ! [output] 771 integer i_cent_C(n_ecp) ! [output] 772 double precision c_exp(nz_ecp) ! [output] 773 double precision c_coef(nz_ecp) ! [output] 774c 775 logical on_ecp, on_so 776 integer geom ! geometry handle 777 integer ecpid ! lexical basis set handle 778 integer soid ! lexical so pot handle 779 integer icent ! counter for centers 780 integer iucent ! unique map of icent 781 integer p_nprim ! counter/pointer into exp/coeff array 782 integer num_ecp ! ecp center as counted 783 integer F_cont ! first contraction on center iucent 784 integer L_cont ! last contraction on center iucent 785 integer iucont ! contraction counter 786 integer type ! function type (-1 = local, 0-lval is non-local) 787 integer nprim ! number of prims in a contraction 788 integer ncoef ! number of coefficients in a contraction 789 integer iexp ! pointer into exndcf for exponents 790 integer icfp ! pointer into exndcf for coefficients 791 integer irexp ! pointer into exndcf for r-exponents 792 integer n0,n1,n2 ! r-exponent count 793 integer n3,n4 ! r-exponent count 794 integer ie ! ecp/so counter 795cc AJL/Begin/SPIN ECPs 796 integer mychan ! channel of contracted ECP 797cc AJL/End 798c 799#include "bas_exndcf_sfn.fh" 800#include "bas_ibs_sfn.fh" 801c 802* write(6,*)' inside build' 803 ecpid = ecpidin + BASIS_HANDLE_OFFSET 804 soid = soidin + BASIS_HANDLE_OFFSET 805c... zero arrays 806 call ifill((5*(l_ecp+2)*n_ecp*2),0,n_prim_C,1) 807 call ifill(((l_ecp+2)*n_ecp*2),0,n_coef_C,1) 808 call ifill(((l_ecp+2)*n_ecp*2),0,ind_C,1) 809 call ifill(((l_ecp+2)*n_ecp*2),0,ind_z,1) 810 call ifill(n_ecp,0,l_C,1) 811 call ifill(n_ecp,0,i_cent_C,1) 812 call dfill(nz_ecp,0.0d00,c_exp,1) 813 call dfill(nz_ecp,0.0d00,c_coef,1) 814c 815 if (o_both) then 816 if (.not. 817 & ecpso_list_ncent(geom,soidin,ecpidin,num_ecp,i_cent_C)) 818 & call errquit('ecp_ptrs: ecpso_list failed',911, BASIS_ERR) 819 else if (o_ecp) then 820 if (.not.ecp_list_ncent(geom,ecpidin,num_ecp,i_cent_C)) 821 & call errquit('ecp_ptrs: ecp_list failed',911, BASIS_ERR) 822 else if (o_so) then 823 if (.not.so_list_ncent(geom,soidin,num_ecp,i_cent_C)) 824 & call errquit('ecp_ptrs: ecpso_list failed',911, BASIS_ERR) 825 else 826 call errquit('ecp_ptrs: should not get here',911, BASIS_ERR) 827 endif 828 if (num_ecp.ne.n_ecp) call errquit 829 & ('ecp_ptrs: mismatch in num_ecp and n_ecp', 830 & (num_ecp-n_ecp), BASIS_ERR) 831c 832 p_nprim = 0 833 do ie = 1,n_ecp 834 icent = i_cent_C(ie) 835 if (o_ecp) then 836 on_ecp = bas_isce(geom,ecpidin,icent) 837 else 838 on_ecp = .false. 839 endif 840 if (o_so) then 841 on_so = bas_isce(geom,soidin,icent) 842 else 843 on_so = .false. 844 endif 845 if (on_ecp) then 846 iucent = sf_ibs_ce2uce(icent,ecpid) 847 F_cont = infbs_tags(TAG_FCONT,iucent,ecpid) 848 L_cont = infbs_tags(TAG_LCONT,iucent,ecpid) 849 do iucont = F_cont, L_cont 850 type = infbs_cont(CONT_TYPE,iucont,ecpid) 851 nprim = infbs_cont(CONT_NPRIM,iucont,ecpid) 852 ncoef = nprim*infbs_cont(CONT_NGEN,iucont,ecpid) 853 if (nprim.ne.ncoef) then 854 write(6,*)'general contraction ecp basis are invalid now' 855 call errquit('int_ecp_build_ecp_ptrs: error',911, 856 & BASIS_ERR) 857 endif 858 iexp = infbs_cont(CONT_IEXP,iucont,ecpid) 859 icfp = infbs_cont(CONT_ICFP,iucont,ecpid) 860 irexp = infbs_cont(CONT_IREXP,iucont,ecpid) 861cc AJL/Begin/SPIN ECPs 862 mychan = infbs_cont(CONT_CHANNEL,iucont,ecpid) 863 if (mychan.eq.channel.or.mychan.eq.0) then 864cc AJL: If mychan = 0, applies the ECP to all channels 865 if ((p_nprim+nprim).gt.nz_ecp) call errquit 866 & ('int_ecp_build_ecp_ptrs:too many coefficients',911, 867 & BASIS_ERR) 868 call ecp_get_n3( 869 & c_exp((p_nprim+1)), 870 & dbl_mb(mb_exndcf(iexp,ecpid)), 871 & c_coef((p_nprim+1)), 872 & dbl_mb(mb_exndcf(icfp,ecpid)), 873 & dbl_mb(mb_exndcf(irexp,ecpid)), 874 & nprim,n0,n1,n2,n3,n4) 875 n_prim_C(0,type,ie,1) = n0 876 n_prim_C(1,type,ie,1) = n1 877 n_prim_C(2,type,ie,1) = n2 878 n_prim_C(3,type,ie,1) = n3 879 n_prim_C(4,type,ie,1) = n4 880 ind_C(type,ie,1) = p_nprim+1 881 ind_Z(type,ie,1) = p_nprim+1 882 n_coef_c(type,ie,1) = nprim 883 l_C(ie) = max(type,l_C(ie)) 884 p_nprim = p_nprim+nprim 885 endif 886cc AJL/End 887 enddo 888 endif 889 if (on_so) then 890 iucent = sf_ibs_ce2uce(icent,soid) 891 F_cont = infbs_tags(TAG_FCONT,iucent,soid) 892 L_cont = infbs_tags(TAG_LCONT,iucent,soid) 893 do iucont = F_cont, L_cont 894 type = infbs_cont(CONT_TYPE,iucont,soid) 895 nprim = infbs_cont(CONT_NPRIM,iucont,soid) 896 ncoef = nprim*infbs_cont(CONT_NGEN,iucont,soid) 897 if (nprim.ne.ncoef) then 898 write(6,*)'general contraction so basis are invalid now' 899 call errquit('int_ecp_build_ecp_ptrs: error',911, 900 & BASIS_ERR) 901 endif 902 iexp = infbs_cont(CONT_IEXP,iucont,soid) 903 icfp = infbs_cont(CONT_ICFP,iucont,soid) 904 irexp = infbs_cont(CONT_IREXP,iucont,soid) 905 if ((p_nprim+nprim).gt.nz_ecp) call errquit 906 & ('int_ecp_build_ecp_ptrs:too many coefficients',911, 907 & BASIS_ERR) 908 call ecp_get_n3( 909 & c_exp((p_nprim+1)), 910 & dbl_mb(mb_exndcf(iexp,soid)), 911 & c_coef((p_nprim+1)), 912 & dbl_mb(mb_exndcf(icfp,soid)), 913 & dbl_mb(mb_exndcf(irexp,soid)), 914 & nprim,n0,n1,n2,n3,n4) 915 n_prim_C(0,type,ie,2) = n0 916 n_prim_C(1,type,ie,2) = n1 917 n_prim_C(2,type,ie,2) = n2 918 n_prim_C(3,type,ie,2) = n3 919 n_prim_C(4,type,ie,2) = n4 920 ind_C(type,ie,2) = p_nprim+1 921 ind_Z(type,ie,2) = p_nprim+1 922 n_coef_c(type,ie,2) = nprim 923 l_C(ie) = max(type,l_C(ie)) 924 p_nprim = p_nprim+nprim 925 enddo 926 endif 927 enddo 928c 929cc AJL/Begin/SPIN ECPS (Debug) 930c call print_ecp_ptrs(n_ecp, 931c & l_ecp,nz_ecp, 932c & channel, ! current channel [alpha/beta] 933c & channels, ! total number of channels 934c & n_prim_C, 935c & n_coef_C, 936c & ind_C, 937c & ind_Z, 938c & l_C, 939c & c_exp, 940c & c_coef) 941cc AJL/End 942c 943 end 944 subroutine print_ecp_ptrs(n_ecp, 945 & l_ecp,nz_ecp, 946cc AJL/Begin/SPIN ECPs 947 & channel, ! current channel [alpha/beta] 948 & channels, ! total number of channels 949cc AJL/End 950 & n_prim_C, 951 & n_coef_C, 952 & ind_C, 953 & ind_Z, 954 & l_C, 955 & c_exp, 956 & c_coef) 957 implicit none 958c 959 integer n_ecp ! [input] number of ecp centers 960 integer l_ecp ! [input] maximum angular momentum in ecp basis 961 integer nz_ecp ! [input] number of prims/coeffs in stored data structure. 962cc AJL/Begin/SPIN ECPs 963 integer channel ! [input] Current spin channel used for ECPs 964 integer channels ! [input] Total number of spin channels used 965cc AJL/End 966 integer n_prim_C(0:4,-1:l_ecp,n_ecp,2) ! [output] 967 integer n_coef_C(-1:l_ecp,n_ecp,2) ! [output] 968 integer ind_C(-1:l_ecp,n_ecp,2) ! [output] 969 integer ind_Z(-1:l_ecp,n_ecp,2) ! [output] 970 integer l_C(n_ecp) ! [output] 971 double precision c_exp(nz_ecp) ! [output] 972 double precision c_coef(nz_ecp) ! [output] 973* 974 integer i, j, k, l 975* integer low0, high0, low1, high1, low2, high2 976* integer ir 977 integer pe, pek 978* 979 write(6,*)' print_ecp_ptrs: start' 980c 981cc AJL/Begin/SPIN ECPs 982 if (channels.eq.1) then ! ECP variable for number of channels 983 write(6,*) 'Both Channels' 984 else if (channel.eq.1) then ! Spin-polarised check? 985 write(6,*) 'Alpha Channel' 986 else if (channel.eq.2) then 987 write(6,*) 'Beta Channel' 988 endif 989cc AJL/End 990c 991 write(6,*)' exponents and coefficients array' 992 do i = 1,nz_ecp 993 write(6,10000)i,c_exp(i),i,c_coef(i) 994 enddo 99510000 format(1x,'exp(',i5,') =',f12.6,2x,'coeff(',i5,') =',f12.6) 996 do i = 1,n_ecp 997 write(6,*)' l_c(',i,')',l_c(i) 998 enddo 999 do l = 1,2 1000 do i=1,n_ecp 1001 do j=-1,l_ecp 1002 write(6,*)' n_coef_C(',j,',',i,',',l,') = ', 1003 & n_coef_C(j,i,l) 1004 write(6,*)' ind_C(',j,',',i,',',l,') = ', 1005 & ind_C(j,i,l) 1006 write(6,*)' ind_z(',j,',',i,',',l,') = ', 1007 & ind_z(j,i,l) 1008 do k = 0,4 1009 write(6,*)' n_prim_C(', 1010 & k,',',j,',',i,',',l,') =', 1011 & n_prim_C(k,j,i,l) 1012 enddo 1013 enddo 1014 enddo 1015 enddo 1016 write(6,*)' --- ' 1017 write(6,*)' i = 1:',n_ecp 1018 do l = 1,2 1019 if (l.eq.1) 1020 & write(6,*)' ecp' 1021 if (l.eq.2) 1022 & write(6,*)' sop' 1023 do i=1,n_ecp 1024 write(6,*)' angular momentum max on center ', 1025 & i,'is',l_c(i) 1026 write(6,*)' j = -1:',l_ecp 1027 do j=-1,l_ecp 1028 write(6,*)' lval = ',j 1029 pe = ind_c(j,i,l)-1 1030 do k=1,n_prim_C(0,j,i,l) 1031 pek = pe + k 1032 write(6,*)' 0 ',c_exp(pek),c_coef(pek) 1033 enddo 1034 pe = pe + n_prim_C(0,j,i,l) 1035 do k=1,n_prim_C(1,j,i,l) 1036 pek = pe + k 1037 write(6,*)' 1 ',c_exp(pek),c_coef(pek) 1038 enddo 1039 pe = pe + n_prim_C(1,j,i,l) 1040 do k=1,n_prim_C(2,j,i,l) 1041 pek = pe + k 1042 write(6,*)' 2 ',c_exp(pek),c_coef(pek) 1043 enddo 1044 pe = pe + n_prim_C(2,j,i,l) 1045 do k=1,n_prim_C(3,j,i,l) 1046 pek = pe + k 1047 write(6,*)' 3 ',c_exp(pek),c_coef(pek) 1048 enddo 1049 pe = pe + n_prim_C(3,j,i,l) 1050 do k=1,n_prim_C(4,j,i,l) 1051 pek = pe + k 1052 write(6,*)' 4 ',c_exp(pek),c_coef(pek) 1053 enddo 1054 write(6,*)' ' 1055 enddo 1056 write(6,*)' ' 1057 enddo 1058 enddo 1059 write(6,*)' print_ecp_ptrs: finish' 1060 end 1061 subroutine ecp_get_n3( 1062 & c_exp, 1063 & ecp_exp, 1064 & c_coef, 1065 & ecp_coef, 1066 & grexp, 1067 & nprim,new0,new1,new2,new3,new4) 1068 implicit none 1069#include "errquit.fh" 1070c 1071c data structures in the ecp code assume all 0 exponents, all 1 exponents 1072c and all 2 exponents are contiguous. Not enforced in the input. 1073c 1074c 1075 integer nprim, new0, new1, new2, new3, new4 1076 double precision c_exp(nprim) 1077 double precision ecp_exp(nprim) 1078 double precision c_coef(nprim) 1079 double precision ecp_coef(nprim) 1080 double precision grexp(nprim) 1081c 1082 integer i, j, ival, iptr 1083 integer new(0:4) 1084c 1085c assumes no general contraction in the ecp basis 1086c 1087*:debugn3: write(6,*)' ' 1088*:debugn3: write(6,*)'-----------------------------------------------------' 1089*:debugn3: write(6,*)' ecp exp ' 1090*:debugn3: call output (ecp_exp,1,nprim,1,1,nprim,1,1) 1091*:debugn3: write(6,*)' ecp coef ' 1092*:debugn3: call output (ecp_coef,1,nprim,1,1,nprim,1,1) 1093*:debugn3: write(6,*)' grexp ' 1094*:debugn3: call output (grexp,1,nprim,1,1,nprim,1,1) 1095 call ifill(5,0,new,1) 1096 iptr = 0 1097 do i = 0,4 1098 do j = 1,nprim 1099 ival = int(grexp(j) + 0.00002d00) 1100*:debugn3: write(6,*)' i, j, ival ',i, j, ival 1101 if (ival.lt.0.or.ival.gt.4) then 1102 write(6,*)' ival = ',ival 1103 write(6,*)' grexp(j) = ',grexp(j) 1104 call errquit 1105 & ('ecp_get_n3: r-exponent not equal to 0,1,2,3,4',911, 1106 & BASIS_ERR) 1107 endif 1108 if (ival.eq.i) then 1109 new(i) = new(i) + 1 1110 iptr = iptr + 1 1111*:debugn3: write(6,*)' i, j, ival,iptr,new ',i, j, ival,iptr,new 1112 c_exp(iptr) = ecp_exp(j) 1113 c_coef(iptr) = ecp_coef(j) 1114 endif 1115 enddo 1116 enddo 1117*:debugn3: write(6,*) ' iprt after loop,nprim ',iptr,nprim 1118 new0 = new(0) 1119 new1 = new(1) 1120 new2 = new(2) 1121 new3 = new(3) 1122 new4 = new(4) 1123*:debugn3: write(6,*)' c_exp ' 1124*:debugn3: call output (c_exp,1,nprim,1,1,nprim,1,1) 1125*:debugn3: write(6,*)' c_coef ' 1126*:debugn3: call output (c_coef,1,nprim,1,1,nprim,1,1) 1127*:debugn3: write(6,*)'-----------------------------------------------------' 1128*:debugn3: write(6,*)' ' 1129c 1130 end 1131 block data ecp_init_bd 1132c 1133c Block data structure to initialize common block for ecp 1134c 1135#include "ecp_nwc.fh" 1136 data h_xyzecp /0/ ! MA handle for ecp center coordinates 1137 data k_xyzecp /0/ ! MA index for ecp center coordinates 1138 data h_ecp_e /0/ ! MA handle for ecp exponents 1139 data k_ecp_e /0/ ! MA index for ecp exponents 1140 data h_ecp_c /0/ ! MA handle for ecp coefficients 1141 data k_ecp_c /0/ ! MA index for ecp coefficients 1142 data h_ecp_nprim_c /0/ ! MA handle for n_prim_C (see ecp_integral) 1143 data k_ecp_nprim_c /0/ ! MA index for n_prim_C (see ecp_integral) 1144 data h_ecp_ncoef_c /0/ ! MA handle for n_coef_C (see ecp_integral) 1145 data k_ecp_ncoef_c /0/ ! MA index for n_coef_C (see ecp_integral) 1146 data h_ecp_ind_c /0/ ! MA handle for int_C (see ecp_integral) 1147 data k_ecp_ind_c /0/ ! MA index for int_C (see ecp_integral) 1148 data h_ecp_l_c /0/ ! MA handle for l_C (see ecp_integral) 1149 data k_ecp_l_c /0/ ! MA index for l_C (see ecp_integral) 1150 data h_ecp_c2s /0/ ! MA handle for c2s routines 1151 data k_ecp_c2s /0/ ! MA index for c2s routines 1152 data h_ecp_lip /0/ ! MA handle for ecp center index 1153 data k_ecp_lip /0/ ! MA index for ecp center index 1154cc AJL/Begin/SPIN ECPs 1155 data h_ecp_e_beta /0/ ! MA handle for ecp exponents 1156 data k_ecp_e_beta /0/ ! MA index for ecp exponents 1157 data h_ecp_c_beta /0/ ! MA handle for ecp coefficients 1158 data k_ecp_c_beta /0/ ! MA index for ecp coefficients 1159 data h_ecp_nprim_c_beta /0/ ! MA handle for n_prim_C (see ecp_integral) 1160 data k_ecp_nprim_c_beta /0/ ! MA index for n_prim_C (see ecp_integral) 1161 data h_ecp_ncoef_c_beta /0/ ! MA handle for n_coef_C (see ecp_integral) 1162 data k_ecp_ncoef_c_beta /0/ ! MA index for n_coef_C (see ecp_integral) 1163 data h_ecp_ind_c_beta /0/ ! MA handle for ind_C (see ecp_integral) 1164 data k_ecp_ind_c_beta /0/ ! MA index for ind_C (see ecp_integral) 1165 data h_ecp_ind_z_beta /0/ ! MA handle for ind_Z (see ecp_integral) 1166 data k_ecp_ind_z_beta /0/ ! MA index for ind_Z (see ecp_integral) 1167 data h_ecp_l_c_beta /0/ ! MA handle for l_C (see ecp_integral) 1168 data k_ecp_l_c_beta /0/ ! MA index for l_C (see ecp_integral) 1169 data h_ecp_lip_beta /0/ ! MA handle for ecp center index 1170 data k_ecp_lip_beta /0/ ! MA index for ecp center index 1171cc AJL/End 1172 data n_zeta_c /0/ ! length of ecp exp/coef array 1173 data l_ecp /0/ ! high ang for ecp basis 1174 data n_ecp /0/ ! number of ecp centers (from API) 1175cc AJL/Begin/SPIN ECPs 1176 data ecp_channels /0/ ! number of ecp channels [Alpha/Beta] 1177cc AJL/End 1178 data init_ecp_init /.false./ ! logical saying if ecp is init-ed 1179 end 1180 subroutine int_ecp_hf1( 1181 & xyza,expa,coefa,a_nprim,a_ngen,La,ictra, 1182 & xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb, 1183 & ecp_ints,sz_ints,scr,lscr,dryrun) 1184 implicit none 1185#include "mafdecls.fh" 1186#include "ecp_nwc.fh" 1187* 1188 integer a_nprim, a_ngen, La, ictra 1189 integer b_nprim, b_ngen, Lb, ictrb 1190 double precision expa(a_nprim), expb(b_nprim) 1191 double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen) 1192 integer sz_ints ! [input] buffer size for ecp_ints 1193 integer lscr ! [input] length of scratch array 1194 double precision xyza(3), xyzb(3) ! [input] a and b center coords. 1195 double precision ecp_ints(sz_ints) ! [output] ecp integrals 1196 double precision scr(lscr) ! [scratch] array 1197 logical dryrun ! [input] compute vs calculate memory requirements. 1198* 1199* write(6,*)' lscr IN ecp_hf1:',lscr 1200* if (.not.dryrun) then 1201* write(6,*)' int_ecp_hf1: coords a ',xyza 1202* write(6,*)' int_ecp_hf1: coords b ',xyzb 1203* endif 1204c 1205 call ecp_integral( 1206 & xyza, 1207 & expa, 1208 & coefa, 1209 & a_nprim,a_ngen,La,ictra, 1210 & xyzb, 1211 & expb, 1212 & coefb, 1213 & b_nprim,b_ngen,Lb,ictrb, 1214 & dbl_mb(k_xyzecp), 1215 & dbl_mb(k_ecp_e),dbl_mb(k_ecp_c), 1216 & int_mb(k_ecp_nprim_c), 1217 & int_mb(k_ecp_ncoef_c), ! new name is n_colc_C 1218 & int_mb(k_ecp_ind_z), 1219 & int_mb(k_ecp_ind_c), 1220 & n_zeta_c, 1221 & n_zeta_c, 1222 & int_mb(k_ecp_l_c), 1223 & int_mb(k_ecp_lip), 1224 & n_ecp,l_ecp, 1225 & 0, 1226 & dbl_mb(k_ecp_c2s),mem_c2s, 1227 & ecp_ints,sz_ints,1, ! nblk 1 for ecp integrals only 1228 & dryrun,scr,lscr, 1229 & 0) ! ibug 1230* 1231 end 1232 subroutine intdd_ecp_hf1( 1233 & xyza,expa,coefa,a_nprim,a_ngen,La,ictra, 1234 & xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb, 1235 & ecp_grad,sz_grad,nat,scr,lscr,dryrun) 1236 implicit none 1237#include "mafdecls.fh" 1238#include "ecp_nwc.fh" 1239* 1240 integer a_nprim, a_ngen, La, ictra 1241 integer b_nprim, b_ngen, Lb, ictrb 1242 integer nat 1243 double precision expa(a_nprim), expb(b_nprim) 1244 double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen) 1245 integer sz_grad ! [input] buffer size for ecp_grad 1246 integer lscr ! [input] length of scratch array 1247 double precision xyza(3), xyzb(3) ! [input] a and b center coords. 1248 double precision ecp_grad(sz_grad,3,3,(nat*(nat-1)/2+nat)) ! [output] ecp integrals 1249 double precision scr(lscr) ! [scratch] array 1250 logical dryrun ! [input] compute vs calculate memory requirements. 1251* 1252* write(6,*)' lscr IN d_ecp_hf1:',lscr 1253* if (.not.dryrun) then 1254* write(6,*)' intd_ecp_hf1: coords a ',xyza 1255* write(6,*)' intd_ecp_hf1: coords b ',xyzb 1256* endif 1257c 1258 call ecp_hessian( 1259 & xyza, 1260 & expa, 1261 & coefa, 1262 & a_nprim,a_ngen,La,ictra, 1263 & xyzb, 1264 & expb, 1265 & coefb, 1266 & b_nprim,b_ngen,Lb,ictrb, 1267 & dbl_mb(k_xyzecp), 1268 & dbl_mb(k_ecp_e),dbl_mb(k_ecp_c), 1269 & int_mb(k_ecp_nprim_c), 1270 & int_mb(k_ecp_ncoef_c), 1271 & int_mb(k_ecp_ind_z), 1272 & int_mb(k_ecp_ind_c), 1273 & n_zeta_c, 1274 & n_zeta_c, 1275 & int_mb(k_ecp_l_c), 1276 & int_mb(k_ecp_lip), 1277 & n_ecp,l_ecp, 1278 & 0, 1279 & dbl_mb(k_ecp_c2s),mem_c2s, 1280 & ecp_grad,sz_grad,1,nat, ! nblock = 1 for ECP only 1281 & dryrun,scr,lscr, 1282 & 0) ! ibug 1283* 1284 end 1285 subroutine intd_ecp_hf1( 1286 & xyza,expa,coefa,a_nprim,a_ngen,La,ictra, 1287 & xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb, 1288 & ecp_grad,sz_grad,nat,scr,lscr,dryrun) 1289 implicit none 1290#include "mafdecls.fh" 1291#include "ecp_nwc.fh" 1292* 1293 integer a_nprim, a_ngen, La, ictra 1294 integer b_nprim, b_ngen, Lb, ictrb 1295 integer nat 1296 double precision expa(a_nprim), expb(b_nprim) 1297 double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen) 1298 integer sz_grad ! [input] buffer size for ecp_grad 1299 integer lscr ! [input] length of scratch array 1300 double precision xyza(3), xyzb(3) ! [input] a and b center coords. 1301 double precision ecp_grad(sz_grad,3,nat) ! [output] ecp integrals 1302 double precision scr(lscr) ! [scratch] array 1303 logical dryrun ! [input] compute vs calculate memory requirements. 1304* 1305* write(6,*)' lscr IN d_ecp_hf1:',lscr 1306* if (.not.dryrun) then 1307* write(6,*)' intd_ecp_hf1: coords a ',xyza 1308* write(6,*)' intd_ecp_hf1: coords b ',xyzb 1309* endif 1310c 1311 call ecp_gradient( 1312 & xyza, 1313 & expa, 1314 & coefa, 1315 & a_nprim,a_ngen,La,ictra, 1316 & xyzb, 1317 & expb, 1318 & coefb, 1319 & b_nprim,b_ngen,Lb,ictrb, 1320 & dbl_mb(k_xyzecp), 1321 & dbl_mb(k_ecp_e),dbl_mb(k_ecp_c), 1322 & int_mb(k_ecp_nprim_c), 1323 & int_mb(k_ecp_ncoef_c), 1324 & int_mb(k_ecp_ind_z), 1325 & int_mb(k_ecp_ind_c), 1326 & n_zeta_c, 1327 & n_zeta_c, 1328 & int_mb(k_ecp_l_c), 1329 & int_mb(k_ecp_lip), 1330 & n_ecp,l_ecp, 1331 & 0, 1332 & dbl_mb(k_ecp_c2s),mem_c2s, 1333 & ecp_grad,sz_grad,1,nat, ! nblock = 1 for ECP only 1334 & dryrun,scr,lscr, 1335 & 0) ! ibug 1336* 1337 end 1338 subroutine intso_hf1( 1339 & xyza,expa,coefa,a_nprim,a_ngen,La,ictra, 1340 & xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb, 1341 & ecp_ints,sz_ints,scr,lscr,dryrun) 1342 implicit none 1343#include "mafdecls.fh" 1344#include "ecp_nwc.fh" 1345* 1346 integer a_nprim, a_ngen, La, ictra 1347 integer b_nprim, b_ngen, Lb, ictrb 1348 double precision expa(a_nprim), expb(b_nprim) 1349 double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen) 1350 integer sz_ints ! [input] buffer size for ecp_ints 1351 integer lscr ! [input] length of scratch array 1352 double precision xyza(3), xyzb(3) ! [input] a and b center coords. 1353 double precision ecp_ints(sz_ints) ! [output] ecp integrals 1354 double precision scr(lscr) ! [scratch] array 1355 logical dryrun ! [input] compute vs calculate memory requirements. 1356* 1357* write(6,*)' lscr IN so_hf1:',lscr 1358c 1359 call ecp_integral( 1360 & xyza, 1361 & expa, 1362 & coefa, 1363 & a_nprim,a_ngen,La,ictra, 1364 & xyzb, 1365 & expb, 1366 & coefb, 1367 & b_nprim,b_ngen,Lb,ictrb, 1368 & dbl_mb(k_xyzecp), 1369 & dbl_mb(k_ecp_e),dbl_mb(k_ecp_c), 1370 & int_mb(k_ecp_nprim_c), 1371 & int_mb(k_ecp_ncoef_c), ! new name is n_colc_C 1372 & int_mb(k_ecp_ind_z), 1373 & int_mb(k_ecp_ind_c), 1374 & n_zeta_c, 1375 & n_zeta_c, 1376 & int_mb(k_ecp_l_c), 1377 & int_mb(k_ecp_lip), 1378 & n_ecp,l_ecp, 1379 & 0, 1380 & dbl_mb(k_ecp_c2s),mem_c2s, 1381 & ecp_ints,sz_ints,3, ! nblk 1 for ecp integrals only 1382 & dryrun,scr,lscr, 1383 & 0) ! ibug 1384* 1385 end 1386 subroutine intd_so_hf1( 1387 & xyza,expa,coefa,a_nprim,a_ngen,La,ictra, 1388 & xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb, 1389 & ecp_grad,sz_grad,nat,scr,lscr,dryrun) 1390 implicit none 1391#include "mafdecls.fh" 1392#include "ecp_nwc.fh" 1393* 1394 integer a_nprim, a_ngen, La, ictra 1395 integer b_nprim, b_ngen, Lb, ictrb 1396 integer nat 1397 double precision expa(a_nprim), expb(b_nprim) 1398 double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen) 1399 integer sz_grad ! [input] buffer size for ecp_grad 1400 integer lscr ! [input] length of scratch array 1401 double precision xyza(3), xyzb(3) ! [input] a and b center coords. 1402 double precision ecp_grad(sz_grad,3,nat) ! [output] ecp integrals 1403 double precision scr(lscr) ! [scratch] array 1404 logical dryrun ! [input] compute vs calculate memory requirements. 1405* 1406* write(6,*)' lscr IN d_ecp_hf1:',lscr 1407* if (.not.dryrun) then 1408* write(6,*)' intd_ecp_hf1: coords a ',xyza 1409* write(6,*)' intd_ecp_hf1: coords b ',xyzb 1410* endif 1411c 1412 call ecp_gradient( 1413 & xyza, 1414 & expa, 1415 & coefa, 1416 & a_nprim,a_ngen,La,ictra, 1417 & xyzb, 1418 & expb, 1419 & coefb, 1420 & b_nprim,b_ngen,Lb,ictrb, 1421 & dbl_mb(k_xyzecp), 1422 & dbl_mb(k_ecp_e),dbl_mb(k_ecp_c), 1423 & int_mb(k_ecp_nprim_c), 1424 & int_mb(k_ecp_ncoef_c), 1425 & int_mb(k_ecp_ind_z), 1426 & int_mb(k_ecp_ind_c), 1427 & n_zeta_c, 1428 & n_zeta_c, 1429 & int_mb(k_ecp_l_c), 1430 & int_mb(k_ecp_lip), 1431 & n_ecp,l_ecp, 1432 & 0, 1433 & dbl_mb(k_ecp_c2s),mem_c2s, 1434 & ecp_grad,sz_grad,3,nat, ! nblock = 3 for SO only 1435 & dryrun,scr,lscr, 1436 & 0) ! ibug 1437* 1438 end 1439 subroutine intdd_so_hf1( 1440 & xyza,expa,coefa,a_nprim,a_ngen,La,ictra, 1441 & xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb, 1442 & ecp_grad,sz_grad,nat,scr,lscr,dryrun) 1443 implicit none 1444#include "mafdecls.fh" 1445#include "ecp_nwc.fh" 1446* 1447 integer a_nprim, a_ngen, La, ictra 1448 integer b_nprim, b_ngen, Lb, ictrb 1449 integer nat 1450 double precision expa(a_nprim), expb(b_nprim) 1451 double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen) 1452 integer sz_grad ! [input] buffer size for ecp_grad 1453 integer lscr ! [input] length of scratch array 1454 double precision xyza(3), xyzb(3) ! [input] a and b center coords. 1455 double precision ecp_grad(sz_grad,3,3,(nat*(nat-1)/2+nat)) ! [output] ecp integrals 1456 double precision scr(lscr) ! [scratch] array 1457 logical dryrun ! [input] compute vs calculate memory requirements. 1458* 1459* write(6,*)' lscr IN d_ecp_hf1:',lscr 1460* if (.not.dryrun) then 1461* write(6,*)' intd_ecp_hf1: coords a ',xyza 1462* write(6,*)' intd_ecp_hf1: coords b ',xyzb 1463* endif 1464c 1465 call ecp_hessian( 1466 & xyza, 1467 & expa, 1468 & coefa, 1469 & a_nprim,a_ngen,La,ictra, 1470 & xyzb, 1471 & expb, 1472 & coefb, 1473 & b_nprim,b_ngen,Lb,ictrb, 1474 & dbl_mb(k_xyzecp), 1475 & dbl_mb(k_ecp_e),dbl_mb(k_ecp_c), 1476 & int_mb(k_ecp_nprim_c), 1477 & int_mb(k_ecp_ncoef_c), 1478 & int_mb(k_ecp_ind_z), 1479 & int_mb(k_ecp_ind_c), 1480 & n_zeta_c, 1481 & n_zeta_c, 1482 & int_mb(k_ecp_l_c), 1483 & int_mb(k_ecp_lip), 1484 & n_ecp,l_ecp, 1485 & 0, 1486 & dbl_mb(k_ecp_c2s),mem_c2s, 1487 & ecp_grad,sz_grad,3,nat, ! nblock = 1 for ECP only 1488 & dryrun,scr,lscr, 1489 & 0) ! ibug 1490* 1491 end 1492c 1493cc AJL/Begin/SPIN ECPs 1494cc All the subroutines designed for use with the beta ECPs 1495 1496 subroutine int_ecp_hf1_beta( 1497 & xyza,expa,coefa,a_nprim,a_ngen,La,ictra, 1498 & xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb, 1499 & ecp_ints,sz_ints,scr,lscr,dryrun) 1500 implicit none 1501#include "mafdecls.fh" 1502#include "ecp_nwc.fh" 1503* 1504 integer a_nprim, a_ngen, La, ictra 1505 integer b_nprim, b_ngen, Lb, ictrb 1506 double precision expa(a_nprim), expb(b_nprim) 1507 double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen) 1508 integer sz_ints ! [input] buffer size for ecp_ints 1509 integer lscr ! [input] length of scratch array 1510 double precision xyza(3), xyzb(3) ! [input] a and b center coords. 1511 double precision ecp_ints(sz_ints) ! [output] ecp integrals 1512 double precision scr(lscr) ! [scratch] array 1513 logical dryrun ! [input] compute vs calculate memory requirements. 1514* 1515* write(6,*)' lscr IN ecp_hf1:',lscr 1516* if (.not.dryrun) then 1517* write(6,*)' int_ecp_hf1: coords a ',xyza 1518* write(6,*)' int_ecp_hf1: coords b ',xyzb 1519* endif 1520c 1521 call ecp_integral( 1522 & xyza, 1523 & expa, 1524 & coefa, 1525 & a_nprim,a_ngen,La,ictra, 1526 & xyzb, 1527 & expb, 1528 & coefb, 1529 & b_nprim,b_ngen,Lb,ictrb, 1530 & dbl_mb(k_xyzecp), 1531 & dbl_mb(k_ecp_e_beta), 1532 & dbl_mb(k_ecp_c_beta), 1533 & int_mb(k_ecp_nprim_c_beta), 1534 & int_mb(k_ecp_ncoef_c_beta), ! new name is n_colc_C 1535 & int_mb(k_ecp_ind_z_beta), 1536 & int_mb(k_ecp_ind_c_beta), 1537 & n_zeta_c, 1538 & n_zeta_c, 1539 & int_mb(k_ecp_l_c_beta), 1540 & int_mb(k_ecp_lip_beta), 1541 & n_ecp,l_ecp, 1542 & 0, 1543 & dbl_mb(k_ecp_c2s),mem_c2s, 1544 & ecp_ints,sz_ints,1, ! nblk 1 for ecp integrals only 1545 & dryrun,scr,lscr, 1546 & 0) ! ibug 1547* 1548 end 1549 subroutine intdd_ecp_hf1_beta( 1550 & xyza,expa,coefa,a_nprim,a_ngen,La,ictra, 1551 & xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb, 1552 & ecp_grad,sz_grad,nat,scr,lscr,dryrun) 1553 implicit none 1554#include "mafdecls.fh" 1555#include "ecp_nwc.fh" 1556* 1557 integer a_nprim, a_ngen, La, ictra 1558 integer b_nprim, b_ngen, Lb, ictrb 1559 integer nat 1560 double precision expa(a_nprim), expb(b_nprim) 1561 double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen) 1562 integer sz_grad ! [input] buffer size for ecp_grad 1563 integer lscr ! [input] length of scratch array 1564 double precision xyza(3), xyzb(3) ! [input] a and b center coords. 1565 double precision ecp_grad(sz_grad,3,3,(nat*(nat-1)/2+nat)) 1566 ! [output] ecp integrals 1567 double precision scr(lscr) ! [scratch] array 1568 logical dryrun ! [input] compute vs calculate memory requirements. 1569* 1570* write(6,*)' lscr IN d_ecp_hf1:',lscr 1571* if (.not.dryrun) then 1572* write(6,*)' intd_ecp_hf1: coords a ',xyza 1573* write(6,*)' intd_ecp_hf1: coords b ',xyzb 1574* endif 1575c 1576 call ecp_hessian( 1577 & xyza, 1578 & expa, 1579 & coefa, 1580 & a_nprim,a_ngen,La,ictra, 1581 & xyzb, 1582 & expb, 1583 & coefb, 1584 & b_nprim,b_ngen,Lb,ictrb, 1585 & dbl_mb(k_xyzecp), 1586 & dbl_mb(k_ecp_e_beta), 1587 & dbl_mb(k_ecp_c_beta), 1588 & int_mb(k_ecp_nprim_c_beta), 1589 & int_mb(k_ecp_ncoef_c_beta), 1590 & int_mb(k_ecp_ind_z_beta), 1591 & int_mb(k_ecp_ind_c_beta), 1592 & n_zeta_c, 1593 & n_zeta_c, 1594 & int_mb(k_ecp_l_c_beta), 1595 & int_mb(k_ecp_lip_beta), 1596 & n_ecp,l_ecp, 1597 & 0, 1598 & dbl_mb(k_ecp_c2s),mem_c2s, 1599 & ecp_grad,sz_grad,1,nat, ! nblock = 1 for ECP only 1600 & dryrun,scr,lscr, 1601 & 0) ! ibug 1602* 1603 end 1604 subroutine intd_ecp_hf1_beta( 1605 & xyza,expa,coefa,a_nprim,a_ngen,La,ictra, 1606 & xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb, 1607 & ecp_grad,sz_grad,nat,scr,lscr,dryrun) 1608 implicit none 1609#include "mafdecls.fh" 1610#include "ecp_nwc.fh" 1611* 1612 integer a_nprim, a_ngen, La, ictra 1613 integer b_nprim, b_ngen, Lb, ictrb 1614 integer nat 1615 double precision expa(a_nprim), expb(b_nprim) 1616 double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen) 1617 integer sz_grad ! [input] buffer size for ecp_grad 1618 integer lscr ! [input] length of scratch array 1619 double precision xyza(3), xyzb(3) ! [input] a and b center coords. 1620 double precision ecp_grad(sz_grad,3,nat) ! [output] ecp integrals 1621 double precision scr(lscr) ! [scratch] array 1622 logical dryrun ! [input] compute vs calculate memory requirements. 1623* 1624* write(6,*)' lscr IN d_ecp_hf1:',lscr 1625* if (.not.dryrun) then 1626* write(6,*)' intd_ecp_hf1: coords a ',xyza 1627* write(6,*)' intd_ecp_hf1: coords b ',xyzb 1628* endif 1629c 1630 call ecp_gradient( 1631 & xyza, 1632 & expa, 1633 & coefa, 1634 & a_nprim,a_ngen,La,ictra, 1635 & xyzb, 1636 & expb, 1637 & coefb, 1638 & b_nprim,b_ngen,Lb,ictrb, 1639 & dbl_mb(k_xyzecp), 1640 & dbl_mb(k_ecp_e_beta), 1641 & dbl_mb(k_ecp_c_beta), 1642 & int_mb(k_ecp_nprim_c_beta), 1643 & int_mb(k_ecp_ncoef_c_beta), 1644 & int_mb(k_ecp_ind_z_beta), 1645 & int_mb(k_ecp_ind_c_beta), 1646 & n_zeta_c, 1647 & n_zeta_c, 1648 & int_mb(k_ecp_l_c_beta), 1649 & int_mb(k_ecp_lip_beta), 1650 & n_ecp,l_ecp, 1651 & 0, 1652 & dbl_mb(k_ecp_c2s),mem_c2s, 1653 & ecp_grad,sz_grad,1,nat, ! nblock = 1 for ECP only 1654 & dryrun,scr,lscr, 1655 & 0) ! ibug 1656* 1657 end 1658cc AJL/End 1659