1* 2* $Id$ 3* 4 logical function raktask_intdd_3c(rtdb) 5 implicit none 6#include "errquit.fh" 7#include "stdio.fh" 8#include "bas.fh" 9#include "geom.fh" 10#include "mafdecls.fh" 11#include "global.fh" 12* 13*::functions 14 logical int_normalize 15 external int_normalize 16*::passed 17 integer rtdb 18*::local 19 integer basis, geom, nbf, cn_nbf_max, nshell 20 integer nat, nat3 21 integer size, sizesq, sizeg, scr_size 22 integer maxg1, maxg2, maxs1, maxs2 23 logical status 24 integer hbuf, kbuf, hscr, kscr 25 integer hfd, kfd, hfdsq, kfdsq, hxyz, kxyz 26 integer hgradp, kgradp 27 integer hgradm, kgradm 28 integer hbufsum, kbufsum 29* 30 raktask_intdd_3c = .false. 31* 32 if (.not.geom_create(geom,'geometry')) call errquit 33 & ('geom create failed',911, GEOM_ERR) 34 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) call errquit 35 & ('geom_rtdb_load failed',911, RTDB_ERR) 36c 37 if (.not.bas_create(basis,'ao basis')) call errquit 38 & ('bas_create failed',911, BASIS_ERR) 39 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) call errquit 40 & ('bas_rtdb_load failed',911, RTDB_ERR) 41c 42 write(luout,*)' geom/basis loaded' 43 call util_flush(luout) 44c 45 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 46c 47 if (.not. bas_print(basis)) 48 $ call errquit(' basis print failed', 0, BASIS_ERR) 49 if (.not. gbs_map_print(basis)) 50 $ call errquit(' gbs_map_print failed', 0, UNKNOWN_ERR) 51c 52 if (.not.bas_numbf(basis,nbf)) call errquit 53 & ('numbf failed',911, BASIS_ERR) 54c 55 if (.not.bas_numcont(basis,nshell)) call errquit 56 & ('numbf failed',911, BASIS_ERR) 57c 58 if (.not.geom_ncent(geom,nat)) stop 'geom_ncent fe' 59 write(luout,*) 'number of atoms ', nat 60 nat3 = 3*nat 61c 62 if (.not.bas_nbf_cn_max(basis,cn_nbf_max)) 63 & stop 'bas_nbf_cn_max' 64c 65 size = (cn_nbf_max**4)*78 66 sizesq = (cn_nbf_max**4)*12*12 67 sizeg = (cn_nbf_max**4)*12 68 size = size + size /10 69 sizesq = sizesq + sizesq/10 70 sizeg = sizeg + sizeg /10 71c 72 call intdd_init(rtdb,1,basis) 73 call int_mem_print() 74 call int_mem_1e(maxg1,maxs1) 75 call int_mem_2e4c(maxg2,maxs2) 76* 77 write(luout,*)' maxg1 :',maxg1 78 write(luout,*)' maxs1 :',maxs1 79 write(luout,*)' maxg2 :',maxg2 80 write(luout,*)' size :',size 81 write(luout,*)' sizesq:',sizesq 82 write(luout,*)' sizeg :',sizeg 83 write(luout,*)' maxs2 :',maxs2 84 call util_flush(luout) 85* 86 scr_size = 2*maxs2 87 status = 88 & ma_alloc_get(mt_dbl,size,'intdd buffer',hbuf,kbuf) 89 status = status.and. 90 & ma_alloc_get(mt_dbl,sizesq,'intdd buffer summed', 91 & hbufsum,kbufsum) 92 status = status.and. 93 & ma_alloc_get(mt_dbl,scr_size,'scr buffer',hscr,kscr) 94 status = status.and. 95 & ma_alloc_get(mt_dbl,size,'intdd fd buffer',hfd,kfd) 96 status = status.and. 97 & ma_alloc_get(mt_dbl,sizesq,'intdd fd sq buffer',hfdsq,kfdsq) 98 status = status.and. 99 & ma_alloc_get(mt_dbl,3*nat,'coords',hxyz,kxyz) 100 status = status.and. 101 & ma_alloc_get(mt_dbl,sizeg,'grad +',hgradp,kgradp) 102 status = status.and. 103 & ma_alloc_get(mt_dbl,sizeg,'grad -',hgradm,kgradm) 104 if (.not.status) stop ' memory alloc failed rak23 (1)' 105* 106 call raktask_intdd_3ca(geom,basis,nbf,nshell,cn_nbf_max, 107 & size,dbl_mb(kbuf), 108 & scr_size,dbl_mb(kscr), 109 & dbl_mb(kgradp), 110 & dbl_mb(kgradm), 111 & dbl_mb(kfd), 112 & dbl_mb(kfdsq), 113 & nat,dbl_mb(kxyz), 114 & size,sizesq,sizeg, 115 & dbl_mb(kbufsum)) 116* 117 call intdd_terminate() 118 raktask_intdd_3c = bas_destroy(basis) 119 raktask_intdd_3c = raktask_intdd_3c.and. 120 & geom_destroy(geom) 121 raktask_intdd_3c = raktask_intdd_3c.and. 122 & ma_free_heap(hscr) 123 raktask_intdd_3c = raktask_intdd_3c.and. 124 & ma_free_heap(hbuf) 125 raktask_intdd_3c = raktask_intdd_3c.and. 126 & ma_free_heap(hfd) 127 raktask_intdd_3c = raktask_intdd_3c.and. 128 & ma_free_heap(hfdsq) 129 raktask_intdd_3c = raktask_intdd_3c.and. 130 & ma_free_heap(hxyz) 131 raktask_intdd_3c = raktask_intdd_3c.and. 132 & ma_free_heap(hgradp) 133 raktask_intdd_3c = raktask_intdd_3c.and. 134 & ma_free_heap(hgradm) 135 raktask_intdd_3c = raktask_intdd_3c.and. 136 & ma_free_heap(hbufsum) 137 end 138 subroutine raktask_intdd_3ca(geom,basis,nbf,nshell,cn_nbf_max, 139 & lbuf,buf,lscr,scr,gradp, gradm,buffd,buffdsq,nat,xyz, 140 & lbuffd, lbuffdsq, lgrad, bufsum) 141 implicit none 142#include "mafdecls.fh" 143#include "stdio.fh" 144#include "bas.fh" 145#include "nwc_const.fh" 146#include "geomP.fh" 147#include "inp.fh" 148c::functions 149 logical rakdd_checktrans 150 external rakdd_checktrans 151c::passed 152 integer geom 153 integer basis 154 integer nbf 155 integer nshell 156 integer cn_nbf_max 157 integer lbuf 158 integer lscr 159 integer nat 160 integer lbuffd, lbuffdsq, lgrad 161 double precision buffd(lbuffd), buffdsq(lbuffdsq) 162 double precision gradp(lgrad), gradm(lgrad) 163 double precision buf(lbuf) 164 double precision bufsum(lbuffdsq) 165 double precision scr(lscr) 166 double precision xyz(3,nat) 167* 168 integer hbufcp, kbufcp 169 integer nzero,ncount,count 170 integer iatom, jatom, katom, latom 171 integer ish, jsh, ksh 172 integer ilo, ihi, inbf 173 integer jlo, jhi, jnbf 174 integer klo, khi, knbf 175 integer nint, ninth, nintg 176 integer idatom(4) 177 integer idatoms(4) 178 integer idatomp(4) 179 integer idatomm(4) 180 integer atoms2move(4) 181 integer num_atoms2move, atom1, atom2 182 integer nat3 183 integer ixyz, zatom 184 double precision thresh, delta, scale, normmax, normmin, norm 185 character*99 string,strings, stringe 186* 187 normmax = 0.0d00 188 normmin = 1.2345678d05 189 call int_acc_high() 190* store original coordintates 191 nat3 = 3*nat 192 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 193* 194 thresh = 1.0d-07 195 delta = 0.001d00 196* 197 write(luout,*)' ',nshell,' total shells ' 198 call util_flush(luout) 199* do ish = 1, nshell 200 do ish = 1, 1 201 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) 202 & stop 'cn2bfr error i' 203 inbf = ihi - ilo + 1 204 strings = ' ' 205 call util_date(strings) 206 if (.not.bas_cn2ce(basis,ish,iatom)) 207 & stop 'bas_cn2ce error i' 208* do jsh = 1, nshell 209 do jsh = 1, 1 210 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) 211 & stop 'cn2bfr error j' 212 jnbf = jhi - jlo + 1 213 if (.not.bas_cn2ce(basis,jsh,jatom)) 214 & stop 'bas_cn2ce error j' 215* do ksh = 1, nshell 216 do ksh = 1, 6 217 if (.not.bas_cn2bfr(basis,ksh,klo,khi)) 218 & stop 'cn2bfr error k' 219 knbf = khi - klo + 1 220 if (.not.bas_cn2ce(basis,ksh,katom)) 221 & stop 'bas_cn2ce error k' 222 if (.not.(iatom.eq.jatom.and. 223 & jatom.eq.katom)) then 224 nint = inbf*jnbf*knbf 225 nintg = nint*12 226 ninth = nint*78 227 call dfill (lbuf,0.0d00,buf,1) 228 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 229 call intdd_2e3c(basis,ish,basis,jsh,ksh, 230 & lscr,scr,lbuf,buf,idatom) 231 call rakdd_checkt_78(nint,buf) 232 write(luout,*)idatom 233 idatoms(1) = idatom(1) 234 idatoms(2) = idatom(2) 235 idatoms(3) = idatom(3) 236 idatoms(4) = idatom(4) 237 write(6,*)' idatom 1:', idatom 238 write(6,*)' idatoms 1:', idatoms 239 if (.not.ma_alloc_get(mt_dbl,lbuf, 240 & 'copy of buf',hbufcp,kbufcp)) stop ' ma alloc fail' 241 call dcopy(lbuf,buf,1,dbl_mb(kbufcp),1) 242 call rakdd_fill(12,nint,dbl_mb(kbufcp),bufsum,idatoms) 243 write(6,*)' idatom 2:', idatom 244 write(6,*)' idatoms 2:', idatoms 245 if (.not.ma_free_heap(hbufcp)) stop ' ma_free fail' 246 call rakdd_print_dd(nint,buf,idatom) 247 if (.not.rakdd_checktrans(nint,bufsum)) then 248 call rakdd_print_dd(nint,buf,idatom) 249 call rakdd_printsum(nint,bufsum,idatoms) 250 endif 251 nzero = 0 252 ncount = 0 253 do count = 1,ninth 254 if (abs(buf(count)).gt.thresh) then 255 ncount = ncount + 1 256 else 257 nzero = nzero + 1 258 endif 259 enddo 260 write(luout,*)nzero,'+',ncount,' != ',ninth 261 if ((nzero+ncount).ne.ninth) 262 & write(luout,*)nzero,'+',ncount,' != ',ninth 263 atoms2move(1) = iatom 264 atoms2move(2) = iatom 265 atoms2move(3) = jatom 266 atoms2move(4) = katom 267 num_atoms2move = 0 268 do atom1 = 1,4 269 do atom2 = 2,4 270 if (atom1.ne.atom2) then 271 if (atoms2move(atom1).eq.atoms2move(atom2)) 272 & atoms2move(atom2) = 0 273 endif 274 enddo 275 enddo 276 num_atoms2move = 0 277 do atom1 = 1,4 278 if (atoms2move(atom1).gt.0) 279 & num_atoms2move = num_atoms2move + 1 280 enddo 281* 282 call dfill (lbuffdsq,0.0d00,buffdsq,1) 283 write(6,*)' atoms2move ',atoms2move 284 do zatom = 1,4 285 if (atoms2move(zatom).gt.0) then 286 do ixyz = 1,3 287 nintg = nint*12 288 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 289 coords(ixyz,atoms2move(zatom),geom) = 290 & coords(ixyz,atoms2move(zatom),geom) + delta 291 call dfill(lgrad,0.0d00,gradp,1) 292 call intd_2e3c(basis,ish,basis,jsh,ksh, 293 & lscr,scr,lgrad,gradp,idatomp) 294 write(string,*)' grad +',ixyz,zatom 295* call rakdd_printgrad(nint,gradp,string) 296 write(6,*)'idatomp:',idatomp 297 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 298 coords(ixyz,atoms2move(zatom),geom) = 299 & coords(ixyz,atoms2move(zatom),geom) - delta 300 call dfill(lgrad,0.0d00,gradm,1) 301 call intd_2e3c(basis,ish,basis,jsh,ksh, 302 & lscr,scr,lgrad,gradm,idatomm) 303 write(string,*)' grad -',ixyz,zatom 304* call rakdd_printgrad(nint,gradm,string) 305 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 306 write(6,*)'idatomm:',idatomm 307 call daxpy(nintg,-1.0d00,gradm,1,gradp,1) 308 write(string,*)' grad diff',ixyz,zatom 309* call rakdd_printgrad(nint,gradp,string) 310 scale = 1.0d00/(2.0d00*delta) 311 call dscal(nintg,scale,gradp,1) 312 write(string,*)' grad diff scaled',ixyz,zatom 313* call rakdd_printgrad(nint,gradp,string) 314 call rakdd_grad_fill(gradp,idatomp, 315 & idatomm,ixyz, 316 & zatom,idatoms, 317 & atoms2move(zatom),buffdsq,nint) 318 enddo 319 endif 320 enddo ! zatom 321 call rakdd_fill_b(buffdsq,nint) ! check if symmetric 322 call rakdd_printsq(buffdsq,nint) 323* call rakdd_print_both(80,buffdsq,bufsum,nint) 324 call rakdd_compare3c( 325 & ish,jsh,ksh, 326 & (nint*12*12),buffdsq,bufsum,norm) 327 normmax = max(normmax,norm) 328 normmin = min(normmin,norm) 329 endif ! 4 atoms the same 330 enddo ! ksh 331 enddo ! jsh 332 string = ' ' 333 stringe = ' ' 334 call util_date(stringe) 335 write(string,'(1x,i5,a3,a27,a2,a27)') 336 & ish,' ::', 337 & strings(1:25),'::', 338 & stringe(1:25) 339 do jsh = 1,inp_strlen(string) 340 ksh = ichar(string(jsh:jsh)) 341 if (ksh.eq.10) string(jsh:jsh) = ' ' 342 enddo 343 write(luout,'(1x,a)')string(1:inp_strlen(string)) 344 enddo ! ish 345 write(luout,'(1x,a,1pd20.10)') 346 & ' maximum difference norm over all quartets:',normmax 347 write(luout,'(1x,a,1pd20.10)') 348 & ' minimum difference norm over all quartets:',normmin 349* 35010000 format(1x,a,4(i4),1x,a,4(i3),1x,a,4(1x,i5),1x,a,i6,/, 351 & 47x,a,4(1x,i5),1x,a,i2) 352 end 353 subroutine rakdd_compare3c( 354 & ish,jsh,ksh, 355 & len,buffd,buf,norm) 356 implicit none 357#include "stdio.fh" 358#include "mafdecls.fh" 359 integer ish,jsh,ksh 360 integer len 361 double precision buffd(len), buf(len) 362* 363 double precision ddot 364 external ddot 365* 366 double precision norm, thresh 367 integer h_diff, k_diff 368* 369 thresh = 1.0d-07 370 if (.not.ma_alloc_get(mt_dbl,len, 371 & 'diff buf',h_diff,k_diff)) stop ' ma alloc fail' 372 373 call dcopy(len,buf,1,dbl_mb(k_diff),1) 374 call daxpy(len,-1.0d00,buffd,1,dbl_mb(k_diff),1) 375 norm = ddot(len,dbl_mb(k_diff),1,dbl_mb(k_diff),1) 376* if (norm.gt.thresh) 377* & 378 write(luout,*)' ',ish,jsh,ksh, 379 & ' norm = ',norm 380 if (.not.ma_free_heap(h_diff)) stop ' ma free fail' 381 end 382 subroutine rakdd_printsq(buffsq,nint) 383 implicit none 384 integer nint 385 double precision buffsq (3,4,3,4,nint) 386c 387 integer ia, ja, ix, jx, int 388 logical doit 389c 390 do int = 1,nint 391 do ia = 1,4 392 do ja = 1,4 393 if (ia.le.ja) then 394 do ix = 1,3 395 do jx = 1,3 396 doit = .true. 397 if (ia.eq.ja) doit = ix.le.jx 398 if (doit) then 399 write(6,10000)ix,ia,jx,ja,int, 400 & buffsq(ix,ia,jx,ja,int) 401 endif 402 enddo 403 enddo 404 endif 405 enddo 406 enddo 407 enddo 40810000 format(1x,'dd(',4(i4,','),i4,') =',1pd20.10) 409 end 410