1module gdma 2 3! Distributed Multipole Analysis for Gaussian Wavefunctions 4! 5! Copyright (C) 2005-08 Anthony J. Stone 6! 7! This program is free software; you can redistribute it and/or 8! modify it under the terms of the GNU General Public License 9! as published by the Free Software Foundation; either version 2 10! of the License, or (at your option) any later version. 11! 12! This program is distributed in the hope that it will be useful, 13! but WITHOUT ANY WARRANTY; without even the implied warranty of 14! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! GNU General Public License for more details. 16! 17! You should have received a copy of the GNU General Public License 18! along with this program; if not, write to 19! the Free Software Foundation, Inc., 51 Franklin Street, 20! Fifth Floor, Boston, MA 02110-1301, USA. 21! 22! This version has been modified by Andy Simmonett (03/16) to link 23! into Psi4, rather than serve as a standalone executable. 24 25USE input 26use iso_c_binding 27 28USE dma 29USE atom_grids, ONLY: debug_grid => debug 30USE timing, ONLY: start_timer, timer, time_and_date 31IMPLICIT NONE 32 33INTEGER, PARAMETER :: dp=kind(1d0) 34 35CHARACTER(LEN=100) :: file 36CHARACTER(LEN=80) :: buffer 37CHARACTER(LEN=20) :: key 38CHARACTER(LEN=8) :: whichd="SCF" 39CHARACTER(LEN=24) :: datestring 40 41! Maximum number of sites is number of atoms + nextra 42INTEGER :: nextra=16 43INTEGER :: ncoorb, maxl, cmax, nprim, nx, num, ich, mul 44INTEGER, ALLOCATABLE :: shell_type(:) 45INTEGER :: i, j, k, kp=0 46LOGICAL :: eof, fchk, first, ok=.false. 47INTEGER open_status, infile 48 49REAL(dp), ALLOCATABLE :: densty(:,:), dtri(:) 50INTEGER :: ir=5 ! Input stream 51 52LOGICAL :: verbose=.false., debug(0:2)=.false. 53 54CONTAINS 55 56 57subroutine run_gdma(c_outfilename, c_datfilename) bind(c, name='run_gdma') 58!character(kind=c_char,len=1), intent(in) :: c_outfilename 59CHARACTER(kind=C_CHAR) :: c_outfilename(*), c_datfilename(*) 60character(len=:), allocatable :: outfilename, datfilename 61integer i, nchars 62integer outfile 63 64i = 1 65do 66 if (c_outfilename(i) == c_null_char) exit 67 i = i + 1 68end do 69nchars = i - 1 ! Exclude null character from Fortran string 70!allocate(character(len=nchars) :: outfilename) 71outfilename = (repeat(' ', nchars)) 72outfilename = transfer(c_outfilename(1:nchars), outfilename) 73i = 1 74do 75 if (c_datfilename(i) == c_null_char) exit 76 i = i + 1 77end do 78nchars = i - 1 ! Exclude null character from Fortran string 79!allocate(character(len=nchars) :: datfilename) 80datfilename = (repeat(' ', nchars)) 81datfilename = transfer(c_datfilename(1:nchars), datfilename) 82 83! 84! Added file IO (ACS 03/16) 85! 86infile = 51 87outfile = 52 88open (unit=infile, file=datfilename, status='old', & 89 iostat=open_status, action='read', position='rewind') 90if ( open_status /= 0 ) then 91 write(outfile, *) 'Could not open GDMA input for reading.', & 92 'unit = ', infile 93 stop 94endif 95open (unit=outfile, file=outfilename, status='old', & 96 iostat=open_status, action='write', position='append') 97if ( open_status /= 0 ) then 98 write(outfile, *) 'Could not open psi4 output for writing.', & 99 'unit = ', infile 100 stop 101endif 102!!! 103write(outfile, "(15x,a/)") & 104 " G D M A", & 105 " by Anthony Stone", & 106 " version 2.2.06, 22 June 2011", & 107 "Distributed Multipoles from Gaussian wavefunctions" 108 109call time_and_date(datestring) 110write(outfile, "(/2A)") "Starting at ", datestring 111 112call start_timer 113 114punchfile="dma.punch" 115nat=0 116fchk=.false. 117first=.true. 118do 119 call read_line(eof, infile) 120 if (eof) exit 121 call readu(key) 122 select case(key) 123 case("","NOTE","!") 124 cycle 125 case("VERBOSE") 126 verbose=.true. 127 case("QUIET") 128 debug=.false. 129 verbose=.false. 130 case("DEBUG") 131 debug(0)=.true. 132 do while (item<nitems) 133 call readi(k) 134 if (k>0) then 135 debug(k)=.true. 136 else 137 debug(-k)=.false. 138 end if 139 end do 140 debug_grid=.true. 141 verbose=.true. 142 case("ANGSTROM") 143 rfact=bohr 144 case("BOHR") 145 rfact=1d0 146 case("SI") 147 Qfactor(0)=echarge 148 do k=1,20 149 Qfactor(k)=Qfactor(k-1)*bohr 150 end do 151 case("AU") 152 Qfactor=1d0 153 case("COMMENT","TITLE") 154 call reada(buffer) 155 write(outfile, "(/a/)") trim(buffer) 156 case("DENSITY") 157 if (fchk) call die & 158 ("Specify density to use before reading data file",.true.) 159 call readu(whichd) 160 case("FILE","READ") 161 nat=0 162 fchk=.false. 163 ok=.false. 164 first=.true. 165 if (allocated(dtri)) deallocate(dtri) 166 do while (item<nitems) 167 call readu(key) 168 select case(key) 169 case("DENSITY") 170 call readu(whichd) 171 case default 172 call reread(-1) 173 call reada(file) 174 open(unit=9,file=file,status="old",iostat=k) 175 if (k .ne. 0) then 176 call die("Can't open file "//file,.true.) 177 endif 178 end select 179 end do 180 ir=9 181 call get_data(whichd,ok,outfile) 182 close(9) 183 ir=5 184 fchk=.true. 185 case ("HERE") 186 call get_data(whichd,ok,outfile) 187 fchk=.true. 188 case("NAMES") 189 if (.not. fchk) call die & 190 ("Read data file before specifying atom names",.false.) 191 call read_line(eof, infile) 192 do i=1,nat 193 call geta(name(i)) 194 end do 195 case("GO","START","MULTIPOLES") 196 if (.not. ok) then 197 call die (trim(whichd)//" density not found",.false.) 198 endif 199 if (first) then 200 ! convert density matrix to triangular form 201 allocate(dtri(nx)) 202 k=0 203 do i=1,num 204 do j=1,i 205 k=k+1 206 dtri(k)=densty(i,j) 207 end do 208 end do 209 deallocate(densty) 210 first=.false. 211 endif 212 write(outfile, "(//2A/)") "Using "//trim(whichd)//" density matrix", & 213 " from file "//trim(file) 214 call dma_main(dtri,kp,infile, outfile) 215 !call timer 216 case("RESET") 217 nat=0 218 fchk=.false. 219 ok=.false. 220 first=.true. 221 deallocate(dtri) 222 whichd="SCF" 223 case("FINISH") 224 exit 225 case default 226 call die("Keyword "//trim(key)//" not recognized",.true.) 227 end select 228end do 229 230call time_and_date(datestring) 231write(outfile, "(/2A)") "Finished at ", datestring 232close(outfile) 233close(infile) 234end subroutine run_gdma 235 236 237!----------------------------------------------------------------------- 238 239SUBROUTINE get_data(whichd,ok,outfile) 240 241IMPLICIT NONE 242 243CHARACTER(LEN=*), INTENT(IN) :: whichd 244INTEGER, INTENT(IN) :: outfile 245LOGICAL, INTENT(OUT) :: ok 246 247INTEGER :: atom, i, j, k, n, nn, aok 248REAL(dp) :: e, rt3v2, td(5,6), tf(7,10), tg(9,15) 249REAL(dp), ALLOCATABLE :: temp(:,:) 250LOGICAL eof 251CHARACTER :: text*40, buffer*80, ww*2, density_header*24, type*1 252CHARACTER(LEN=8) :: dummy 253REAL(dp), PARAMETER :: PI=3.14159265358979d0 254CHARACTER(LEN=2), DIMENSION(54) :: element=(/"H ", "He", & 255 "Li", "Be", "B ", "C ", "N ", "O ", "F ", "Ne", & 256 "Na", "Mg", "Al", "Si", "P ", "S ", "Cl", "Ar", & 257 "K ", "Ca", "Sc", "Ti", "V ", "Cr", "Mn", "Fe", "Co", "Ni", & 258 "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", & 259 "Rb", "Sr", "Y ", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", & 260 "Ag", "Cd", "In", "Sn", "Sb", "Te", "I ", "Xe"/) 261REAL(dp), PARAMETER :: rt2=1.4142135623731d0, & 262 rt3=1.73205080756888d0, rt5=2.23606797749979d0, rt7=2.64575131106459d0 263REAL(dp), PARAMETER :: rt10=rt2*rt5, rt14=rt2*rt7, rt21=rt3*rt7, & 264 rt35=rt5*rt7 265INTEGER, PARAMETER :: v400=1, v040=2, v004=3, v310=4, v301=5, & 266 v130=6, v031=7, v103=8, v013=9, v220=10, v202=11, v022=12, & 267 v211=13, v121=14, v112=15 268CHARACTER(LEN=5) :: label(-5:5)=(/"h(s)","g(s)","f(s)","d(s)","sp ", & 269 "s ","p ","d ","f ","g ","h "/) 270 271! Conversion from normalised spherical form to normalised Cartesian 272! Schlegel & Frisch, IJQC (1995) 54, 83-87. 273rt3v2=rt3/2d0 274! d functions 275! 1 2 3 4 5 6 276! xx yy zz xy xz yz 277! 200 020 002 110 101 011 278td(1,:)=(/-0.5d0, -0.5d0, 1d0, 0d0, 0d0, 0d0/) 279td(2,:)=(/0d0, 0d0, 0d0, 0d0, 1d0, 0d0/) 280td(3,:)=(/0d0, 0d0, 0d0, 0d0, 0d0, 1d0/) 281td(4,:)=(/rt3v2, -rt3v2, 0d0, 0d0, 0d0, 0d0/) 282td(5,:)=(/0d0, 0d0, 0d0, 1d0, 0d0, 0d0/) 283 284! f functions 285! 1 2 3 4 5 6 7 8 9 10 286! xxx yyy zzz xxy xxz xyy yyz xzz yzz xyz 287! 300 030 003 210 201 120 021 102 012 111 288tf(:,:)=0d0 289! 30 290tf(1,3)=1d0; tf(1,5)=-1.5d0/sqrt(5d0); tf(1,7)=-1.5d0/sqrt(5d0) 291! 31c ( F+1 in Gaussian notation ) 292tf(2,1)=-sqrt(3d0/8d0); tf(2,6)=-sqrt(1.2d0)/4d0; tf(2,8)=sqrt(1.2d0) 293! 31s ( F-1 ) 294tf(3,2)=-sqrt(3d0/8d0); tf(3,4)=-sqrt(1.2d0)/4d0; tf(3,9)=sqrt(1.2d0) 295! 32c ( F+2 ) 296tf(4,5)=sqrt(0.75d0); tf(4,7)=-sqrt(0.75d0) 297! 32s 298tf(5,10)=1d0 299! 33c 300tf(6,1)=sqrt(10d0)/4d0; tf(6,6)=-0.75d0*sqrt(2d0) 301! 33s 302tf(7,2)=-sqrt(10d0)/4d0; tf(7,4)=0.75d0*sqrt(2d0) 303 304! g functions 305! 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 306! xxxx yyyy zzzz xxxy xxxz xyyy yyyz xzzz yzzz xxyy xxzz yyzz xxyz xyyz xyzz 307! 400 040 004 310 301 130 031 103 013 220 202 022 211 121 112 308! 12 23 D 23 23,12 13 23,12 12 D 23 12 309tg=0d0 310! 40 311tg(1,v400)=0.375d0; tg(1,v040)=0.375d0; tg(1,v004)=1d0 312tg(1,v220)=0.75d0*rt3/rt35; tg(1,v202)=-3d0*rt3/rt35; tg(1,v022)=-3d0*rt3/rt35 313! 41c 314tg(2,v103)=rt10/rt7; tg(2,v301)=-0.75d0*rt10/rt7; tg(2,v121)=-0.75d0*rt2/rt7 315! 41s 316tg(3,v013)=rt10/rt7; tg(3,v031)=-0.75d0*rt10/rt7; tg(3,v211)=-0.75d0*rt2/rt7 317! 42c 318tg(4,v202)=1.5d0*rt3/rt7; tg(4,v022)=-1.5d0*rt3/rt7; tg(4,v400)=-rt5/4d0; tg(4,v040)=rt5/4d0 319! 42s 320tg(5,v112)=3d0/rt7; tg(5,v310)=-rt5/(2d0*rt7); tg(5,v130)=-rt5/(2d0*rt7) 321! 43c 322tg(6,v301)=rt10/4d0; tg(6,v121)=-0.75d0*rt2 323! 43s 324tg(7,v031)=-rt10/4d0; tg(7,v211)=0.75d0*rt2 325! 44c 326tg(8,v400)=rt35/8d0; tg(8,v040)=rt35/8d0; tg(8,v220)=-0.75*rt3 327! 44s 328tg(9,v310)=rt5/2d0; tg(9,v130)=-rt5/2d0 329 330! select case(whichg) 331! case("94","G94") 332! rfact=1d0 333! case("98","G98","03","G03") 334! rfact=1d0 335! end select 336 337ok=.false. 338call stream(ir) 339read (ir,"(A80/A80)") title(1), title(2) 340density_header="Total "//trim(whichd)//" Density" 341if (verbose) then 342 write(outfile,"(a/a/a/)") "Gaussian header:", trim(title(1)), trim(title(2)) 343 write(outfile,"(3a)") 'Looking for "', trim(density_header), '"' 344end if 345do 346 read (ir,"(A)",iostat=k) buffer 347 if (k .ne. 0) then 348 call stream(5) 349 exit 350 end if 351 if (debug(0)) write(outfile, "(a)") buffer 352 text=buffer(1:40) 353 type=buffer(44:44) 354 ww=buffer(48:49) 355 if (ww .eq. "N=") then 356 read(buffer,"(55X,I6)") nn 357 else 358 nn=0 359 endif 360 if (debug(0)) write(outfile, "(a)") text 361 select case(text) 362 case("") 363 cycle 364 case("END","End","end") 365 call stream(5) 366 exit 367 case("Number of atoms") 368 read(buffer,"(55X,I6)") nat 369 if (verbose) write(outfile, "(i0,a)") nat, " atoms" 370 if (allocated(zan)) deallocate(zan,c) 371 allocate(zan(nat),c(3,nat),stat=aok) 372 if (aok>0) call die("Can't allocate atom arrays") 373 maxcen=nat 374 maxs=maxcen+nextra ! Arbitrary limit on number of sites 375 if (allocated(name)) deallocate(name) 376 allocate(name(maxs),stat=aok) 377 if (aok>0) call die("Can't allocate site-name array") 378 case("Charge") 379 read(buffer,"(55X,I6)") ich 380 if (verbose) write(outfile, "(a,i0)") "Charge ", ich 381 case("Multiplicity") 382 read(buffer,"(55X,I6)") mul 383 if (verbose) write(outfile, "(a,i0)") "Multiplicity ", mul 384 case("Number of basis functions") 385 read(buffer,"(55X,I6)") ncoorb 386 ! This number may be increased following conversion from 387 ! spherical to cartesian 388 if (verbose) write(outfile, "(i0,a)") ncoorb, " basis functions" 389 case("Number of contracted shells") 390 read(buffer,"(55X,I6)") nshell 391 if (verbose) write(outfile, "(i0,a)") nshell, " shells" 392 if (allocated(kstart)) & 393 deallocate(kstart,katom,ktype,kng,kloc,kmin,kmax,shell_type) 394 allocate (kstart(nshell), katom(nshell+1), ktype(nshell), & 395 kng(nshell), kloc(nshell), kmin(nshell), kmax(nshell), & 396 shell_type(nshell),stat=aok) 397 if (aok>0) call die("Can't allocate shell arrays") 398 shell_type=0 399 case("Highest angular momentum") 400 read(buffer,"(55X,I6)") maxl 401 if (verbose) write(outfile, "(a,i0)") "Highest angular momentum ", maxl 402 if (maxl .gt. 4) call die & 403 ("Sorry -- GDMA can only handle s, p, d, f and g basis functions",.false.) 404 case("Largest degree of contraction") 405 read(buffer,"(55X,I6)") cmax 406 if (verbose) write(outfile, "(a,i0)") "Largest contraction depth ", cmax 407 if (cmax .gt. 16) call die & 408 ("Sorry -- maximum contraction depth is 16",.false.) 409 case("Number of primitive shells") 410 read(buffer,"(55X,I6)") nprim 411 if (verbose) write(outfile, "(i0,a)") nprim, " primitive shells" 412 if (allocated(ex)) deallocate(ex,cs,cp) 413 allocate(ex(nprim), cs(nprim), cp(nprim), stat=aok) 414 if (aok>0) then 415 call die("Can't allocate arrays for primitives.") 416 end if 417 ex=0d0; cs=0d0; cp=0d0 418 case("Atomic numbers") 419 call read_line(eof) 420 do i=1,nat 421 call geti(k) 422 name(i)=element(k) 423 end do 424 if (verbose) write(outfile, "(a,20a3/(17x,20a3))") & 425 "Atoms: ", name(1:nat) 426 case("Nuclear charges") 427 call read_line(eof) 428 do i=1,nat 429 call getf(zan(i)) 430 end do 431 if (verbose) write(outfile, "(a,20i3/(16x,20i3))") & 432 "Nuclear charges:", nint(zan(1:nat)) 433 case("Current cartesian coordinates") 434 call read_line(eof) 435 if (verbose) write(outfile, "(a)") "Atom Z Position (a.u.)" 436 do i=1,nat 437 do j=1,3 438 call getf(c(j,i),rfact) 439 end do 440 if (verbose) write(outfile, "(a3,i4,3f10.5)") name(i), nint(zan(i)), c(:,i) 441 end do 442 case("Shell types") 443 call read_line(eof) 444 ! num is the original number of basis functions. n is the 445 ! increased number when conversion from spherical to 446 ! cartesian is taken into account. 447 num=0 448 n=0 449 do i=1,nshell 450 call geti(shell_type(i)) 451 select case(shell_type(i)) 452 case(0) 453 num=num+1; n=n+1 454 case(1) 455 num=num+3; n=n+3 456 case(-1) 457 num=num+4; n=n+4 458 case(2) 459 num=num+6; n=n+6 460 case(-2) 461 num=num+5; n=n+6 462 case(3) 463 num=num+10; n=n+10 464 case(-3) 465 num=num+7; n=n+10 466 case(4) 467 num=num+15; n=n+15 468 case(-4) 469 num=num+9; n=n+15 470 end select 471 end do 472 if (verbose) write(outfile, "(i0,a)") num, " basis functions" 473 if ((verbose) .and. n>num) & 474 write(outfile, "(a,i0,a)") "(", n, " after conversion to cartesian)" 475 maxbfn=n 476 if (allocated(iax)) deallocate(iax) 477 allocate(iax(n+1), stat=aok) 478 if (aok>0) then 479 call die("Can't allocate IAX array") 480 end if 481 case("Number of primitives per shell") 482 call read_line(eof) 483 k=1 484 do i=1,nshell 485 call geti(j) 486 kstart(i)=k 487 k=k+j 488 kng(i)=j 489 end do 490 if (verbose) then 491 write(outfile,"(a,20i3/(19x,20i3))") "Contraction depths:", kng(1:nshell) 492 write(outfile,"(a,i0)") "Total number of primitives required: ", k-1 493 end if 494 if (k .ne. nprim+1) call die & 495 ("Shell contractions do not match number of primitives",.false.) 496 case("Shell to atom map") 497 call read_line(eof) 498 do i=1,nshell 499 call geti(katom(i)) 500 end do 501 if (verbose) then 502 write(outfile,"(a,120i3)") "shell to atom", katom(1:nshell) 503 endif 504 case("Primitive exponents") 505 call read_line(eof) 506 do i=1,nprim 507 call getf(ex(i)) 508 end do 509! print "(a,20F10.6)", "primitive exps", ex(1:nprim) 510 case("Contraction coefficients") 511 call read_line(eof) 512 do i=1,nshell 513 do j=kstart(i),kstart(i)+kng(i)-1 514 call getf(e) 515 if (shell_type(i) .eq. 1) then 516 cp(j)=e 517 else 518 cs(j)=e 519 end if 520 ! select case(shell_type(i)) 521 ! case(-1,0) 522 ! cs(j)=e 523 ! case(1) 524 ! cp(j)=e 525 ! case(2,-2,3,-3,4,-4,5,-5) 526 ! cd(j)=e 527 ! end select 528 end do 529 end do 530 case("P(S=P) Contraction coefficients") 531 call read_line(eof) 532 do i=1,nshell 533 do j=kstart(i),kstart(i)+kng(i)-1 534 call getf(e) 535 if (shell_type(i) .eq. -1) then 536 cp(j)=e 537 end if 538 end do 539 end do 540! case("Alpha MO coefficients","Beta MO coefficients") 541! if (debug(2)) then 542! print "(/a)", text 543! allocate(temp(ncoorb,ncoorb)) 544! do i=1,ncoorb 545! do j=1,ncoorb 546! call getf(temp(j,i)) 547! end do 548! end do 549! call matwrtt(temp,1,ncoorb,1,ncoorb,format="6f12.5") 550! deallocate(temp) 551! end if 552 case default 553 if (text .eq. density_header) then 554 ncoorb=n 555 nx=n*(n+1)/2 556 allocate(densty(n,n),temp(n,n)) 557 call read_line(eof) 558 do i=1,num 559 do j=1,i 560 call getf(densty(i,j)); densty(j,i)=densty(i,j) 561 end do 562 end do 563 ok=.true. 564 else 565 ! Ignore this section 566 if (nn .gt. 0) then 567 if (type .eq. "I") then 568 do i=1,(nn+5)/6 569 read(ir,"(A8)") dummy 570 end do 571 else if (type .eq. "R") then 572 do i=1,(nn+4)/5 573 read(ir,"(A8)") dummy 574 end do 575 end if 576 endif 577 endif 578 end select 579end do 580 581if (verbose) then 582 atom=0 583 do i=1,nshell 584 if (katom(i) .ne. atom) then 585 atom=katom(i) 586 write(outfile, "(a)") name(atom) 587 end if 588 write(outfile, "(a,i0,3x,a)") "Shell ", i, label(shell_type(i)) 589 do j=kstart(i),kstart(i)+kng(i)-1 590 select case (shell_type(i)) 591 case (-1) 592 write(outfile, "(i10,f16.8, 2f14.8)") j, ex(j), cs(j), cp(j) 593 case(0) 594 write(outfile, "(i10,f16.8, 2f14.8)") j, ex(j), cs(j) 595 case(1) 596 write(outfile, "(i10,f16.8, 2f14.8)") j, ex(j), cp(j) 597 case(2,-2,3,-3,4,-4) 598 write(outfile, "(i10,f16.8, 2f14.8)") j, ex(j), cs(j) 599 end select 600 end do 601 end do 602end if 603! We use unnormalized primitive functions, so we transfer the 604! normalising factor to the contraction coefficients. This is 605! the factor for z^n exp(-e*r^2). General formula is 606! (4e)^(n/2).(2e/pi)^{3/4}/sqrt{(2n-1)!!} 607do i=1,nshell 608 do j=kstart(i),kstart(i)+kng(i)-1 609 e=ex(j) 610 select case(abs(shell_type(i))) 611 case(0,1) 612 cs(j)=cs(j)*sqrt(sqrt((2d0*e/pi)**3)) 613 cp(j)=cp(j)*sqrt(4d0*e*sqrt((2d0*e/pi)**3)) 614 case(5) ! h shell 615 cs(j)=cs(j)*(4d0*e)**2*sqrt(4d0*e*sqrt((2d0*e/pi)**3)/945d0) 616 case(4) ! g shell 617 cs(j)=cs(j)*(4d0*e)**2*sqrt(sqrt((2d0*e/pi)**3)/105d0) 618 case(3) ! f shell 619 cs(j)=cs(j)*4d0*e*sqrt(4d0*e*sqrt((2d0*e/pi)**3)/15d0) 620 case(2) ! d shell 621 cs(j)=cs(j)*4d0*e*sqrt(sqrt((2d0*e/pi)**3)/3d0) 622 end select 623 end do 624end do 625 626if (.not. ok) return 627! call matwrtt(densty,1,num,1,num,format='5F10.5', cols=5) 628 629! Deal with shell types, transforming from spherical to cartesian 630! basis if necessary 631k=0 632do i=1,nshell 633 kloc(i)=k+1 ! First basis function for shell i 634 select case(shell_type(i)) 635 case(-1) ! sp shell 636 kmin(i)=1 637 kmax(i)=4 638 ktype(i)=2 639 case(0) ! s shell 640 kmin(i)=1 641 kmax(i)=1 642 ktype(i)=1 643 case(1) ! p shell 644 kmin(i)=2 645 kmax(i)=4 646 ktype(i)=2 647 case(2,-2) ! d shell 648 kmin(i)=5 649 kmax(i)=10 650 ktype(i)=3 651 if (shell_type(i) .lt. 0) then ! Spherical d shell 652 temp(1:num,1:k)=densty(1:num,1:k) 653 temp(1:num,k+1:k+6)=matmul(densty(1:num,k+1:k+5),td) 654 temp(1:num,k+7:num+1)=densty(1:num,k+6:num) 655 num=num+1 656 densty(1:k,1:num)=temp(1:k,1:num) 657 densty(k+1:k+6,1:num)=matmul(transpose(td),temp(k+1:k+5,1:num)) 658 densty(k+7:num,1:num)=temp(k+6:num-1,1:num) 659 endif 660 case(3,-3) ! f shell 661 kmin(i)=11 662 kmax(i)=20 663 ktype(i)=4 664 if (shell_type(i) .lt. 0) then ! Spherical f shell 665 temp(1:num,1:k)=densty(1:num,1:k) 666 temp(1:num,k+1:k+10)=matmul(densty(1:num,k+1:k+7),tf) 667 if (i<nshell) temp(1:num,k+11:num+3)=densty(1:num,k+8:num) 668 num=num+3 669 densty(1:k,1:num)=temp(1:k,1:num) 670 densty(k+1:k+10,1:num)=matmul(transpose(tf),temp(k+1:k+7,1:num)) 671 if (i<nshell) densty(k+11:num,1:num)=temp(k+8:num-3,1:num) 672 endif 673 case(4,-4) ! g shell 674 kmin(i)=21 675 kmax(i)=35 676 ktype(i)=5 677 ! print "(a,i0,a,i0)", "num = ", num, " k = ", k 678 if (shell_type(i) .lt. 0) then ! Spherical g shell 679 temp(1:num,1:k)=densty(1:num,1:k) 680 temp(1:num,k+1:k+15)=matmul(densty(1:num,k+1:k+9),tg) 681 if (i<nshell) temp(1:num,k+16:num+6)=densty(1:num,k+10:num) 682 num=num+6 683 densty(1:k,1:num)=temp(1:k,1:num) 684 densty(k+1:k+15,1:num)=matmul(transpose(tg),temp(k+1:k+9,1:num)) 685 if (i<nshell) densty(k+16:num,1:num)=temp(k+10:num-6,1:num) 686 endif 687 case default 688 write (buffer,"(a,i0)") "Unrecognized or unimplemented shell type ", i 689 call die(trim(buffer),.false.) 690 end select 691 k=k+kmax(i)-kmin(i)+1 692end do 693if (k .ne. n .or. num .ne. n) call die & 694 ("Mismatch in number of basis functions",.false.) 695if (debug(1)) call matwrtt(densty,1,n,1,n,format='5F10.5', iformat="5i10", cols=5) 696 697deallocate(temp) 698 699END SUBROUTINE get_data 700 701!---------------------------------------------------------------- 702 703SUBROUTINE matwrtt(c,i1,i2,j1,j2,format,cols,iformat) 704 705! Print rows I1 to I2, columns J1 to J2, of the lower triangle of 706! the symmetric matrix C 707 708IMPLICIT NONE 709REAL(dp) :: c(:,:) 710INTEGER i1, i2, j1, j2 711CHARACTER(LEN=*), OPTIONAL :: format, iformat 712INTEGER, OPTIONAL :: cols 713 714INTEGER i, j, jstart, jfinis, ncols 715CHARACTER(LEN=20) :: fmt, ifmt 716 717fmt="1p6g12.4" 718ifmt="12i12" 719if (present(format)) fmt=format 720if (present(iformat)) ifmt=iformat 721ncols=6 722if (present(cols)) ncols=cols 723 724jfinis=j1-1 725do while (jfinis .lt. j2) 726 jstart=jfinis+1 727 jfinis=min(j2,jfinis+ncols) 728 write (6,"(/1x,"//ifmt//")") (j,j=jstart,jfinis) 729 write (6,'(1x)') 730 do i=max(i1,jstart),i2 731 write (6,"(1x,i3,1x,"//fmt//")") & 732 i,(c(i,j),j=jstart,min(i,jfinis)) 733 end do 734end do 735 736END SUBROUTINE matwrtt 737 738END module gdma 739