1 program testbasis 2c $Id$ 3 implicit none 4#include "errquit.fh" 5c 6#include "nwc_const.fh" 7#include "bas.fh" 8#include "rtdb.fh" 9#include "mafdecls.fh" 10#include "geom.fh" 11#include "inp.fh" 12#include "stdio.fh" 13c 14c::friends functions 15 logical bas_ucontinfo 16 logical bas_getu_exponent 17 logical bas_getu_coeff 18 logical bas_setu_exponent 19 logical bas_setu_coeff 20 external bas_ucontinfo 21 external bas_getu_exponent 22 external bas_getu_coeff 23 external bas_setu_exponent 24 external bas_setu_coeff 25c 26 integer rtdb, geom, basis, ngen,nprim, iang 27* integer ncenters 28 integer sphcart, i,j 29 integer nbf,ibf,iat,icont, type 30 character*256 namename, nametrans 31* character*16 drivtags(20) 32* double precision coords(3,20), charge(20) 33 integer lexcf 34 parameter (lexcf=400) 35 double precision exp(lexcf), coeff(lexcf) 36 logical status 37 double precision expnt_new(3), coeff_new(4,3) 38 logical bas_add_ucnt, bas_geobas_build 39 external bas_add_ucnt, bas_geobas_build 40 data expnt_new/1., 2., 3./ 41 data coeff_new/-1., -2., -3., 0.0, -4., -5., -6., 0.0, -7., -8., 42 $ -9., 0.0/ 43c 44 if (.not. ma_init(MT_DBL, 400 000, 400 000)) 45 & call errquit('testbasis: ma_init failed',911, MA_ERR) 46 47 status = rtdb_open('testbasis.rtdb','empty', rtdb) 48 if (.not.status) call errquit 49 & ('testbasis rtdb open failed',911, RTDB_ERR) 50c 51c write(luout,*)' rtdb handle ', rtdb 52*old status = bas_321g_load(rtdb) 53c 54 open(unit=luin,file='testbasis.nw', 55 & form='formatted', 56 & access='sequential', 57 & status='old', 58 & err=99565) 59 call input_parse(rtdb) 60c 61 if(.not.geom_create(geom,'geometry')) then 62 write(LuOut,*)' error getting geometry handle ' 63 stop ' error ' 64 endif 65c 66 if(.not.geom_rtdb_load(rtdb,geom,'geometry')) 67 & call errquit('error loading geometry',911, GEOM_ERR) 68 status=geom_print(geom) 69c 70 basis = 0 71 if(.not.bas_create(basis,'xy basis')) then 72 write(LuOut,*)' error getting basis handle ' 73 stop ' error ' 74 endif 75 basis = 0 76 if(.not.bas_create(basis,'wz basis')) then 77 write(LuOut,*)' error getting second basis handle ' 78 stop ' error ' 79 endif 80 basis = 0 81 if(.not.bas_create(basis,'ao basis')) then 82 write(LuOut,*)' error getting third basis handle ' 83 stop ' error ' 84 endif 85c 86c 87 status = bas_rtdb_load(rtdb,geom,basis,'ao basis') 88 status = bas_print(basis) 89 status = gbs_map_print(basis) 90c 91 status = bas_continfo(basis,1,iang,nprim,ngen,sphcart) 92 write(LuOut,*)' user:query: status cont 1 ',status 93 write(LuOut,*)' user:query: type cont 1 ',iang 94 write(LuOut,*)' user:query: nprim cont 1 ',nprim 95 write(LuOut,*)' user:query: ngen cont 1 ',ngen 96 write(LuOut,*)' user:query: sphcart cont 1 ',sphcart 97 if (nprim*ngen.gt.lexcf) call errquit 98 & (' lexcf too small ',(nprim*ngen), BASIS_ERR) 99 call dfill(lexcf,0.0d00,exp ,1) 100 call dfill(lexcf,0.0d00,coeff,1) 101 status = bas_get_exponent(basis,1,exp) 102 status = bas_get_coeff(basis,1,coeff) 103 write(LuOut,*)' user exponenents and coefficients ' 104 do 00100 i=1,nprim 105 write(LuOut,*)exp(i),(coeff(i+(j-1)*nprim),j=1,ngen) 10600100 continue 107c 108 status = bas_ucontinfo(basis,1,type,nprim,ngen,sphcart) 109 write(LuOut,*)' unique:query: status cont 1 ',status 110 write(LuOut,*)' unique:query: type cont 1 ',type 111 write(LuOut,*)' unique:query: nprim cont 1 ',nprim 112 if (nprim*ngen.gt.lexcf) call errquit 113 & (' lexcf too small ',(nprim*ngen), BASIS_ERR) 114 write(LuOut,*)' unique:query: ngen cont 1 ',ngen 115 write(LuOut,*)' unique:query: sphcart cont 1 ',sphcart 116 call dfill(lexcf,0.0d00,exp ,1) 117 call dfill(lexcf,0.0d00,coeff,1) 118 status = bas_getu_exponent(basis,1,exp) 119 status = bas_getu_coeff(basis,1,coeff) 120 write(LuOut,*)' unique exponenents and coefficients ' 121 do 00200 i=1,nprim 122 write(LuOut,*)exp(i),(coeff(i+(j-1)*nprim),j=1,ngen) 12300200 continue 124c 125 exp(1) = 565.6589 126 coeff(1) = 6.021023 127 status = bas_setu_exponent(basis,1,exp,(nprim+1)) 128 status = bas_setu_coeff(basis,1,coeff,(nprim*ngen)+1) 129 status = bas_setu_exponent(basis,1,exp,nprim) 130 status = bas_setu_coeff(basis,1,coeff,nprim*ngen) 131c 132c 133 status = bas_ucontinfo(basis,1,type,nprim,ngen,sphcart) 134 write(LuOut,*)' modified ' 135 write(LuOut,*)' unique:query: status cont 1 ',status 136 write(LuOut,*)' unique:query: type cont 1 ',type 137 write(LuOut,*)' unique:query: nprim cont 1 ',nprim 138 write(LuOut,*)' unique:query: ngen cont 1 ',ngen 139 write(LuOut,*)' unique:query: sphcart cont 1 ',sphcart 140 if (nprim*ngen.gt.lexcf) call errquit 141 & (' lexcf too small ',(nprim*ngen), BASIS_ERR) 142 status = bas_getu_exponent(basis,1,exp) 143 status = bas_getu_coeff(basis,1,coeff) 144 write(LuOut,*)' exponenents and coefficients ' 145 do 00300 i=1,nprim 146 write(LuOut,*)exp(i),(coeff(i+(j-1)*nprim),j=1,ngen) 14700300 continue 148c 149c 150c Try adding new contractions on an existing center 151c 152 write(LuOut,*) ' adding 3*3 d function on H' 153 if (.not. bas_add_ucnt(basis, 'H', 2, 3, 3, expnt_new, 154 & expnt_new, coeff_new, 4, 'none', .false.)) 155 & write(LuOut,*) ' basis_add_ucnt failed' 156 if (.not. bas_print(basis)) write(LuOut,*) ' print ?' 157c 158c Try adding new contractions on a new center 159c 160 write(LuOut,*) ' adding 2*3 g function on Cl' 161 if (.not. bas_add_ucnt(basis, 'Cl', 4, 2, 3, expnt_new, 162 $ expnt_new,coeff_new, 4, 'none', .false.)) 163 & write(LuOut,*) ' basis_add_ucnt failed' 164 if (.not. bas_print(basis)) write(LuOut,*) ' print ?' 165c 166c 167 write(LuOut,'(///,1x,a)')' bas_print_all ' 168 status = bas_print_all() 169 call bas_err_info(' testbasis bas_err_info info ') 170c 171 status = bas_high_angular(basis,iang) 172 write(LuOut,*)' high angular momentum ', iang 173 status = bas_version() 174 namename = ' ' 175 nametrans = ' ' 176 status = bas_name(basis,namename, nametrans) 177 write(LuOut,*)' handle : ',basis, status 178 i = inp_strlen(namename) 179 write(LuOut,'(a,a,a)')' name : <', 180 & namename(1:i),'>' 181 i = inp_strlen(nametrans) 182 write(LuOut,*)' translated name: <', 183 & nametrans(1:i),'>' 184c 185 write(LuOut,*)' basis function to ce/cn check ' 186 status = bas_numbf(basis,nbf) 187 do 00400 ibf=1,nbf 188 status = bas_bf2ce(basis,ibf,iat) 189 status = bas_bf2cn(basis,ibf,icont) 190 write(LuOut,*)' basis function ',ibf, 191 & ' is on center ',iat, 192 & ' and in contraction',icont 19300400 continue 194c 195 status = bas_geobas_build(basis) 196 write(LuOut,*)' basis handle ',basis 197 call bas_geomap_check(basis) 198c 199 status = bas_nbf_cn_max(basis,nbf) 200 write(LuOut,*)' max block of nbf on a shell for basis is:',nbf 201c 202 status = bas_nbf_ce_max(basis,nbf) 203 write(LuOut,*)' max block of nbf on an atom for basis is:',nbf 204c 205 status = bas_ncoef_cn_max(basis,nbf) 206 write(LuOut,*)' max num coefs in any contraction ',nbf 207c 208 status = bas_nprim_cn_max(basis,nbf) 209 write(LuOut,*)' max num prims in any contraction ',nbf 210c 211c 212 stop ' testbasis done ' 21399565 continue 214 stop ' error opening testbasis.nw' 215 end 216 subroutine bas_geomap_check(basisin) 217 implicit none 218#include "nwc_const.fh" 219#include "bas.fh" 220#include "basP.fh" 221#include "geobasmapP.fh" 222#include "basdeclsP.fh" 223#include "mafdecls.fh" 224#include "bas_ibs_dec.fh" 225#include "geom.fh" 226#include "stdio.fh" 227c 228 integer basisin 229c 230 logical status 231 integer basis, geom 232 integer ncont, icont, ucont, ucont2 233 integer cent_api, cent_sf, cent_mb 234 integer inat, nat, atom_sf, atom_mb 235 integer bflo_api, bflo_sf, bflo_mb 236 integer bfhi_api, bfhi_sf, bfhi_mb 237 integer cnlo_api, cnlo_sf, cnlo_mb 238 integer cnhi_api, cnhi_sf, cnhi_mb 239c 240#include "bas_ibs_sfn.fh" 241c 242 status = bas_geom(basisin,geom) 243 basis = basisin + BASIS_HANDLE_OFFSET 244 write(LuOut,*)' ' 245 write(LuOut,*)' ' 246 write(LuOut,*)' ' 247 write(LuOut,*)' basis geometry mapping checks ' 248 write(LuOut,*)' basis handle is :',basisin 249 write(LuOut,*)' geometry handle is :',geom 250 write(LuOut,*)' ' 251 write(LuOut,*)' current memory heap pointers' 252 write(LuOut,*)' ibs_cn2ucn ma handle :',ibs_cn2ucn(H_ibs,basis) 253 write(LuOut,*)' ibs_cn2ucn ma pointer:',ibs_cn2ucn(K_ibs,basis) 254 write(LuOut,*)' ibs_cn2ucn size :',ibs_cn2ucn(SZ_ibs,basis) 255 write(LuOut,*)' ' 256 write(LuOut,*)' ibs_cn2ce ma handle :',ibs_cn2ce(H_ibs,basis) 257 write(LuOut,*)' ibs_cn2ce ma pointer :',ibs_cn2ce(K_ibs,basis) 258 write(LuOut,*)' ibs_cn2ce size :',ibs_cn2ce(SZ_ibs,basis) 259 write(LuOut,*)' ' 260 write(LuOut,*)' ibs_ce2uce ma handle :',ibs_ce2uce(H_ibs,basis) 261 write(LuOut,*)' ibs_ce2uce ma pointer:',ibs_ce2uce(K_ibs,basis) 262 write(LuOut,*)' ibs_ce2uce size :',ibs_ce2uce(SZ_ibs,basis) 263 write(LuOut,*)' ' 264 write(LuOut,*)' ibs_cn2bfr ma handle :',ibs_cn2bfr(H_ibs,basis) 265 write(LuOut,*)' ibs_cn2bfr ma pointer:',ibs_cn2bfr(K_ibs,basis) 266 write(LuOut,*)' ibs_cn2bfr size :',ibs_cn2bfr(SZ_ibs,basis) 267 write(LuOut,*)' ' 268 write(LuOut,*)' ibs_ce2cnr ma handle :',ibs_ce2cnr(H_ibs,basis) 269 write(LuOut,*)' ibs_ce2cnr ma pointer:',ibs_ce2cnr(K_ibs,basis) 270 write(LuOut,*)' ibs_ce2cnr size :',ibs_ce2cnr(SZ_ibs,basis) 271 write(LuOut,*)' ' 272 write(LuOut,*)' ' 273 write(LuOut,*)' ' 274c 275 status = bas_numcont(basisin,ncont) 276 status = geom_ncent(geom,nat) 277 write(LuOut,*)' total number of contractions is ',ncont 278 write(LuOut,*)' total number of centers is ',nat 279 write(LuOut,*)' ' 280 write(LuOut,*)' ' 281 write(LuOut,*)' ' 282 write(LuOut,*)' cont 2 uniq cont check ' 283 do icont = 0,ncont 284 ucont = 23234 285 ucont2 = ucont 286 ucont = sf_ibs_cn2ucn(icont,basis) 287 ucont2 = int_mb(mb_ibs_cn2ucn(icont,basis)) 288 write(LuOut,'(a,i5,a,i5,i5)') 289 & ' contraction ',icont, 290 & ' maps to unique contraction ',ucont,ucont2 291 enddo 292 293 write(LuOut,*)' ' 294 write(LuOut,*)' ' 295 write(LuOut,*)' cont 2 center check ' 296 do icont = 0,ncont 297 status = bas_cn2ce(basisin,icont,cent_api) 298 cent_sf = sf_ibs_cn2ce(icont,basis) 299 cent_mb = int_mb(mb_ibs_cn2ce(icont,basis)) 300 write(LuOut,'(a,i5,a,3i5,a)') 301 & ' contraction ',icont,' is on center ', 302 & cent_api,cent_sf,cent_mb, 303 & ' ... (sould be 3 identical numbers)' 304 enddo 305c 306 write(LuOut,*)' ' 307 write(LuOut,*)' ' 308 write(LuOut,*)' center 2 unique center check ' 309 do inat = 1,nat 310 atom_sf = sf_ibs_ce2uce(inat,basis) 311 atom_mb = int_mb(mb_ibs_ce2uce(inat,basis)) 312 write(LuOut,'(a,i5,a,2i5,a)') 313 & ' center ',inat, 314 & ' maps to unique center/tag ',atom_sf,atom_mb, 315 & ' ... (sould be 2 identical numbers)' 316 enddo 317c 318 write(LuOut,*)' ' 319 write(LuOut,*)' ' 320 write(LuOut,*)' contraction 2 basis function range ' 321 do icont = 1,ncont 322 status = bas_cn2bfr(basisin,icont,bflo_api,bfhi_api) 323 bflo_sf = sf_ibs_cn2bfr(1,icont,basis) 324 bfhi_sf = sf_ibs_cn2bfr(2,icont,basis) 325 bflo_mb = int_mb(mb_ibs_cn2bfr(1,icont,basis)) 326 bfhi_mb = int_mb(mb_ibs_cn2bfr(2,icont,basis)) 327 write(LuOut,'(a,i5,a,6i5,a)') 328 & ' contraction ',icont, 329 & ' maps to range of basis functions ', 330 & bflo_api,bfhi_api,bflo_sf,bfhi_sf,bflo_mb,bfhi_mb, 331 & ' ... (sould be 3 sets identical numbers)' 332 enddo 333c 334 write(LuOut,*)' ' 335 write(LuOut,*)' ' 336 write(LuOut,*)' center 2 contraction range ' 337 do inat = 1,nat 338 status = bas_ce2cnr(basisin,inat,cnlo_api,cnhi_api) 339 cnlo_sf = sf_ibs_ce2cnr(1,inat,basis) 340 cnhi_sf = sf_ibs_ce2cnr(2,inat,basis) 341 cnlo_mb = int_mb(mb_ibs_ce2cnr(1,inat,basis)) 342 cnhi_mb = int_mb(mb_ibs_ce2cnr(2,inat,basis)) 343 write(LuOut,'(a,i5,a,6i5,a)') 344 & ' center ',inat, 345 & ' maps to range of contractons', 346 & cnlo_api,cnhi_api,cnlo_sf,cnhi_sf,cnlo_mb,cnhi_mb, 347 & ' ... (sould be 3 sets identical numbers)' 348 enddo 349c 350 end 351