1 subroutine dft_inpana(rtdb) 2c 3c $Id$ 4c 5c********************************************************************* 6c 7c inpana (input analysis) 8c Analyze input to deduce nature of system, and set key flags. 9c Write pertinent information to user output. 10c 11c********************************************************************* 12 implicit none 13#include "errquit.fh" 14c 15#include "stdio.fh" 16#include "rtdb.fh" 17#include "mafdecls.fh" 18#include "global.fh" 19#include "tcgmsg.fh" 20#include "geom.fh" 21#include "bas.fh" 22#include "cdft.fh" 23#include "util.fh" 24#include "util_params.fh" 25c 26 logical even, cksetd, no_prune 27 Integer rtdb ! runtime database handle 28 integer me,ichg,nel 29 integer noc1 30 integer n, na, nb, test_sic 31 integer nshells 32 integer ictr 33 double precision bsr, radang 34 double precision anucl_charg, ckfac, t1 35 double precision ictr_coord(3), ictr_chg,smear_sigma 36 character*3 on_off_1, on_off_2 37 logical oprint_general, oprint_grid, oprint_xc, 38 & oprint_convergence, oprint_tolerances, oprint_sic, 39 , lsigma 40c 41 logical odft,rodft 42 character*4 scftype 43c 44 logical cam_exch,cam_srhf 45 double precision cam_alpha,cam_beta,cam_omega 46c 47 character*9 local_c, nonlocal_c, nineb_c 48 character*10 start_10c, NA_10c, asap_10c 49 character*10 strng1, strng2, strng3, strng4, strng5, strng6 50 character*16 tag, theory 51 character*20 rgridnames(10) 52 logical xc_gotxc 53 external xc_gotxc 54 data rgridnames /'Euler-MacLaurin','Mura-Knowles', 55 . 'Treutler-Ahlrichs','Gauss-Legendre','G-C-interv', 56 , 'Lindh','Chebyshev','Legendre','DE2','DE2D'/ 57c 58 logical status ! FA 59 logical use_nwxc 60 logical dftmp2,out1 61 double precision mp2fac 62c 63 me=ga_nodeid() 64c 65 oprint_general = util_print('general information',print_default) 66 oprint_grid = util_print('grid information',print_default) 67 oprint_xc = util_print('xc information',print_default) 68 oprint_convergence = util_print('convergence information', 69 & print_default) 70 oprint_tolerances = util_print('screening tolerance information', 71 & print_default) 72 oprint_sic = util_print('sic information',print_default) 73c 74c Figure out the number of electrons from the required total 75c charge and the sum of nuclear charges 76c 77 if (.not. rtdb_cget(rtdb, 'dft:theory', 1, theory)) 78 $ call errquit('dft_inpana: theory not specified',0, 79 & RTDB_ERR) 80 if (.not. geom_nuc_charge(geom, anucl_charg)) 81 & call errquit('dft_inpana: geom_nuc_charge failed', 0, 82 & GEOM_ERR) 83 nel = nint(anucl_charg - rcharge) 84 if (nel .le. 0) call errquit 85 $ ('dft_inpana: negative no. of electrons ?', nel, INPUT_ERR) 86 if (abs(anucl_charg - rcharge - dble(nel)) .gt. 1d-8) 87 $ call errquit('dft_inpana: non-integral # of electrons ?', 0, 88 & INPUT_ERR) 89c 90c == Coulomb Attenuation Method (CAM/LC) parameters == 91 if (.not.rtdb_get(rtdb, 'dft:cam_exch', mt_log, 1, cam_exch)) 92 & cam_exch=.false. 93 if (.not.rtdb_get(rtdb, 'dft:cam_srhf', mt_log, 1, cam_srhf)) 94 & cam_srhf=.false. 95 if (.not.rtdb_get(rtdb, 'dft:cam_omega', mt_dbl, 1,cam_omega)) 96 & cam_omega=0.d0 97 if (.not.rtdb_get(rtdb, 'dft:cam_alpha', mt_dbl, 1,cam_alpha)) 98 & cam_alpha=0.d0 99 if (.not.rtdb_get(rtdb, 'dft:cam_beta', mt_dbl, 1, cam_beta)) 100 & cam_beta=0.d0 101c 102c Check to see if calculation type is allowed. 103c 104c Even number of electrons required for RHF. 105c 106 even=mod(nel,2).eq.0 107c 108c == check for restricted open-shell dft == 109 if (.not. rtdb_get(rtdb, 'dft:rodft', mt_log, 1, 110 & rodft))rodft = .false. 111c 112 odft = .false. 113 if (.not. rtdb_cget(rtdb, 'dft:scftype', 1, scftype)) then 114c 115c If the user did not specify whether the calculation should use 116c an open-shell or closed-shell formalism then work it out from 117c the spin multiplicity and potentially from a request for rodft. 118c 119 if (mult.eq.1) then 120 scftype = 'RHF' 121 else 122 odft = .true. 123 if (rodft) then 124 scftype = 'ROHF' 125 else 126 scftype = 'UHF' 127 endif 128 endif 129c 130 else 131c 132c If the user specified what formalism to use then use what the 133c user ordered. Even if that means running a closed shell system 134c with the UHF or ROHF formalism. 135c 136 if (scftype.eq.'UHF') then 137 odft = .true. 138 else if (scftype.eq.'ROHF') then 139c fix for closed shell and rodft 140 if(mult.eq.1) then 141 odft=.false. 142 scftype='RHF' 143 ipol=1 144 else 145 odft = .true. 146 endif 147 148 elseif (scftype.eq.'RHF'.and.mult.ne.1) then 149 odft= .true. 150 scftype='UHF' 151 endif 152c 153 endif 154 if (.not. rtdb_cput(rtdb, 'dft:scftype', 1, scftype)) 155 & call errquit('dft_input: rtdb_put failed', 1500, RTDB_ERR) 156 call dft_cscf_scftype(scftype) 157c 158 if ((.not.even).and.(.not.odft).and.(.not.rodft)) then 159 if (ga_nodeid().eq.0) 160 & write(luout,"(10x,'Number of electrons: ',i10)") nel 161 call errquit( 162 & 'Odd number of electrons. Please specify a 163 &restricted open-shell or open-shell calculation',nel,INPUT_ERR) 164 end if 165c 166c odd # of electrons or not a singlet state --> LSD 167c 168 if ((.not.even).or.(mult.ne.1).or.(theory.eq.'sodft').or. 169 & (scftype.eq.'UHF').or.(scftype.eq.'ROHF')) then 170 ipol=2 171 endif 172 noc(2)=0 173c 174c Calculate number of occupied orbitals. 175c 176 if (ipol.eq.1)then 177 noc1 = nel/2 178 noc(2)= 0 179 noc(1)= noc1 180 else 181c 182c check consistency of no. elec and multiplicity 183c 184 even=mod((nel+mult-1),2).eq.0 185 if (.not.even) then 186 write(LuOut,*)' number of electrons :',nel 187 write(LuOut,*)' multiplicity :',mult 188 call errquit( 189 & ' no. of electrons and multiplicity not compatible',nel, 190 & INPUT_ERR) 191 endif 192 if(mult.gt.0) then 193 noc(2) = (nel - mult + 1)/2 194 noc(1) = nel - noc(2) 195 noc1 = noc(1) + noc(2) 196 else 197 noc(1) = (nel + mult + 1)/2 198 noc(2) = nel - noc(1) 199 noc1 = noc(1) + noc(2) 200 endif 201 endif 202 if (.not. rtdb_put(rtdb, 'dft:ipol', mt_int, 1, ipol)) 203 $ call errquit('inpana: dft:ipol put failed', 0, RTDB_ERR) 204c 205c Check to see if there are enough electrons for this 206c value of the multiplicity. 207c 208 if (noc(2).lt.0)then 209 call errquit('dft: #electrons not valid for multiplicity',mult, 210 & INPUT_ERR) 211 endif 212c 213c write noc (consistent with definition in ddscf) to rtdb 214c 215 if (.not. rtdb_put(rtdb, 'dft:noc', mt_int, 2, noc)) 216 & call errquit('inpana: rtdb_put of noc failed', 0, RTDB_ERR) 217cHvD 218 use_nwxc = util_module_avail("nwxc") 219 if (use_nwxc) then 220 call nwxc_rtdb_load(rtdb,"dft",use_nwxc) 221 endif 222 if (use_nwxc) then 223 call nwxc_getvals("nwxc_has_hfx",out1) 224 if (out1) then 225 call nwxc_getwght("nwxc_wght_hfx",xfac(1)) 226 endif 227 call nwxc_getvals("nwxc_has_mp2c",dftmp2) 228 if (dftmp2) then 229 call nwxc_getwght("nwxc_wght_mp2c",mp2fac) 230 if (.not.rtdb_put(rtdb,'dft:mp2fac', mt_dbl, 1, mp2fac)) 231 & call errquit('dft_inpana: rtdb_put failed', 2902, RTDB_ERR) 232 endif 233 call nwxc_getvals("nwxc_has_cam",cam_exch) 234 if (cam_exch) then 235 call nwxc_get_cam(cam_alpha,cam_beta,cam_omega,cam_srhf) 236 else 237 cam_alpha = 0.0d0 238 cam_beta = 0.0d0 239 cam_omega = 0.0d0 240 endif 241 call nwxc_print() 242 if (cam_exch.and.(ga_nodeid().eq.0)) then 243 write(LuOut,*) 244 write(LuOut,8202) 245 & 'Range-Separation Parameters ' 246 write(LuOut,8203) 247 write(LuOut,8201)cam_alpha,cam_beta,cam_omega,cam_srhf 248 end if ! cam_exch 249 call xc_setcamparam(rtdb,cam_exch,cam_srhf,cam_omega,cam_alpha, 250 & cam_beta) 251 endif 252cHvD 253c 254c Write new data to checkpoint file. 255c 256c Analyze any user specified XC functionals ... set if none 257c 258 cksetd = .true. 259 cksetd = cksetd .and. .not. libxcon 260c 261c Check if user has specified some type of functional, 262c if so, do not set defaults. 263c 264 do n = 1, numfunc 265 if (lcfac(n).or.nlcfac(n).or.lxfac(n).or.nlxfac(n)) 266 & cksetd = .false. 267 enddo 268 if (cksetd)then 269c 270c Set functional defaults. 271c 272 cfac(1) = 1.0d0 273 lcfac(1) = .true. 274 xfac(2) = 1.0d0 275 lxfac(2) = .true. 276c 277c Update rtdb. 278c 279 if (.not. rtdb_put(rtdb, 'dft:cfac', mt_dbl, numfunc, cfac)) 280 $ call errquit('dft_input: rtdb_put failed', 210, RTDB_ERR) 281 if (.not. rtdb_put(rtdb, 'dft:xfac', mt_dbl, numfunc, xfac)) 282 $ call errquit('dft_input: rtdb_put failed', 211, RTDB_ERR) 283 if (.not. rtdb_put(rtdb, 'dft:lcfac', mt_log, numfunc, lcfac )) 284 $ call errquit('dft_input: rtdb_put failed', 9, RTDB_ERR) 285 if (.not. rtdb_put(rtdb, 'dft:lxfac', mt_log, numfunc, lxfac )) 286 $ call errquit('dft_input: rtdb_put failed', 11, RTDB_ERR) 287 endif 288c 289cc AJL/Begin/FDE 290c Repeat for FDE XC. Is this a safe way to check we are using FDE? 291 if (frozemb_fde) then 292c 293c Analyze any user specified XC functionals ... set if none 294c 295 cksetd = .true. 296c 297c Check if user has specified some type of functional, 298c if so, do not set defaults. 299c 300 do n = 1, numfunc 301 if (lcfac_fde(n).or.nlcfac_fde(n).or. 302 & lxfac_fde(n).or.nlxfac_fde(n)) 303 & cksetd = .false. 304 enddo 305 if (cksetd)then 306c 307c Set functional defaults. 308c 309 cfac_fde(1) = 1.0d0 310 lcfac_fde(1) = .true. 311 xfac_fde(2) = 1.0d0 312 lxfac_fde(2) = .true. 313c 314c Update rtdb. 315c 316 if (.not. rtdb_put(rtdb, 317 & 'dft:frozemb:cfac', mt_dbl, numfunc, cfac_fde)) 318 $ call errquit('dft_input: rtdb_put failed', 210, RTDB_ERR) 319 if (.not. rtdb_put(rtdb, 320 & 'dft:frozemb:xfac', mt_dbl, numfunc, xfac_fde)) 321 $ call errquit('dft_input: rtdb_put failed', 211, RTDB_ERR) 322 if (.not. rtdb_put(rtdb, 323 & 'dft:frozemb:lcfac', mt_log, numfunc, lcfac_fde)) 324 $ call errquit('dft_input: rtdb_put failed', 9, RTDB_ERR) 325 if (.not. rtdb_put(rtdb, 326 & 'dft:frozemb:lxfac', mt_log, numfunc, lxfac_fde)) 327 $ call errquit('dft_input: rtdb_put failed', 11, RTDB_ERR) 328 endif 329c 330c Repeat something similar for KE functionals. 331c I'm going to set this up so it can be expanded. 332c 333 cksetd = .true. 334c 335c Check if user has specified some type of functional, 336c if so, do not set defaults. 337c 338 do n = 1, numfunc_fde 339 if (tsfac(n).gt.1.d-8) cksetd = .false. 340 enddo 341 if (cksetd)then 342c 343c Set functional defaults. 344c 345 tsfac(1) = 1.d0 346c 347c Update rtdb. 348c 349 if (.not. rtdb_put(rtdb, 350 & 'dft:frozemb:ts', mt_dbl, numfunc_fde, tsfac)) 351 $ call errquit('dft_input: rtdb_put failed', 210, RTDB_ERR) 352 end if ! tsfac default 353c 354 endif ! frozemb_fde 355c AJL/End 356c 357c Check/set defaults for convergence schemes 358c 359c Three types of convergence speedup schemes: 360c 1) based on number of cycles performed, 361c 2) based on differences of total energy less than some threshold, 362c 3) DIIS plus levelshifting if homo-lumo gap is small. 363c 364c Examples of these might be: 365c 366c 1) based on number of cycles performed, 367c ncydp = 3 368c ndamp = 40 369c ncysh = iterations 370c rlshift = 0.5 371c ncyds = iterations 372c 373c 2) based on differences of total energy being less than some threshold, 374c dampon = 1.d8 375c dampoff = 1.0d-1 376c levlon = 1.d-1 377c levloff = 1.0d-3 378c rlshift = 0.5 379c diison = 1.d-1 380c diisoff = 0.0d0 381c 382c 3) DIIS plus levelshifting if homo-lumo gap is small (default) 383c nodamping = .true. 384c ncydp = 0 385c ndamp = 0 386c ncysh = iterations 387c rlshift = 0.5 388c ncyds = iterations 389c 390 if(nodamping)then 391 damp = .false. 392 ncydp = 0 393 ndamp = 0 394 endif 395 if(nolevelshifting)then 396 levelshift = .false. 397 ncysh = 0 398 rlshift = 0.0 399 endif 400 if(nodiis)then 401 diis = .false. 402 ncyds = 0 403 endif 404 if (damp)then 405c 406c check to make sure either number of damping iterations or 407c energy criterion has been specified. 408c 409 if (ncydp.eq.0.and.dampon.eq.0.0d0)then 410 ncydp = iterations 411 endif 412 endif 413 if (levelshift)then 414c 415c check to make sure either number of levelshifting iterations or 416c energy criterion has been specified. 417c 418 if (ncysh.eq.0.and.levlon.eq.0.0d0)then 419 ncysh = iterations 420 endif 421 endif 422 if (diis)then 423c 424c check to make sure either number of diis iterations or 425c energy criterion has been specified. 426c 427 if (ncyds.eq.0.and.diison.eq.0.0d0)then 428 ncyds = iterations 429 endif 430 endif 431c 432c If convergence input based upon #cycles then turn off energy constraints. 433c 434 if (ncydp.ne.0)then 435 dampon = -999.9 436 dampoff = -999.9 437 endif 438 if (ncysh.ne.0)then 439 levlon = -999.9 440 levloff = -999.9 441 endif 442 if (ncyds.ne.0)then 443 diison = -999.9 444 diisoff = -999.9 445 endif 446c 447c check special case with damping - change 448c default of 2 to "iterations" if no other 449c convergence control specified 450c 451 if (ncysh.eq.0 .and. ncyds.eq.0 .and. ncydp.eq.2) 452 & ncydp = iterations 453c 454c check on 2-e integral and XC grid tolerances 455c 456 call dft_inpanae(rtdb) 457c 458c Check for no pruning; no_prune 459c 460 if (.not. rtdb_get(rtdb, 'dft:no_prune', mt_log, 1, 461 & no_prune))no_prune = .false. 462 if (.not. rtdb_get(rtdb, 'dft:test_sic', mt_int, 1, 463 & test_sic))test_sic = 0 464 lsigma=rtdb_get(rtdb, 'dft:smear_sigma',mt_dbl,1,smear_sigma) 465c 466 if (me.eq.0) then 467c 468c Write to output. 469c 470 if (oprint_general)then 471 write(LuOut,*) 472 call util_print_centered 473 & (LuOut,'General Information',20,.true.) 474 write(LuOut,9020) 475 if (scftype.eq."RHF")then 476 write(LuOut,9050) 477 elseif (scftype.eq."ROHF")then 478 write(LuOut,9052) 479 elseif (scftype.eq."UHF")then 480 write(LuOut,9055) 481 endif 482 write(LuOut,8150)ncenters 483 na = noc(1) 484 nb = noc(2) 485 if (ipol.eq.1)nb = noc(1) 486 ichg = nint(rcharge) 487 write(LuOut,8200)nel,na,nb,ichg,mult 488c 489cFA Writing Number of electrons on rtdb to be used 490c in create_munu4nbo() defined in hnd_efgmap_Z4.F 491 status = rtdb_parallel(.false.) 492 if (.not. rtdb_put(rtdb, 'prop:Nel',mt_int, 493 $ 1,nel)) 494 $ call errquit('prop_input-EFGZ4-nel: rtdb_put failed', 495 $ 555, RTDB_ERR) 496 if (.not. rtdb_put(rtdb, 'prop:Nocc_a',mt_int, 497 $ 1,na)) 498 $ call errquit('prop_input-EFGZ4-Nocc_a: rtdb_put failed', 499 $ 555, RTDB_ERR) 500 if (.not. rtdb_put(rtdb, 'prop:Nocc_b',mt_int, 501 $ 1,nb)) 502 $ call errquit('prop_input-EFGZ4-Nocc_b: rtdb_put failed', 503 $ 555, RTDB_ERR) 504 status = rtdb_parallel(.true.) 505cFA 506 on_off_1 = 'off' 507 if (oskel)on_off_1 = 'on' 508 on_off_2 = 'off' 509 if (oadapt)on_off_2 = 'on' 510 write(LuOut,9056)on_off_1, on_off_2 511 write(LuOut,9030)iterations 512 if (direct)write(LuOut,9110) 513 if (.not. bas_numcont(AO_bas_han, nshells)) 514 & call errquit('rdinput:rdinput:',86, BASIS_ERR) 515 write(LuOut,4001)nbf, nshells 516 if (nbf_cd.gt.0)then 517 write(LuOut,9130) 518 if (.not. bas_numcont(CD_bas_han, nshells_cd)) 519 & call errquit('rdinput:rdinput:',87, BASIS_ERR) 520 write(LuOut,4002)nbf_cd, nshells_cd 521 endif 522 if (nbf_xc.gt.0)then 523 write(LuOut,9120) 524 if (.not. bas_numcont(XC_bas_han, nshells_xc)) 525 & call errquit('rdinput:rdinput:',88, BASIS_ERR) 526 write(LuOut,4003)nbf_xc, nshells_xc 527 endif 528 write(LuOut,9035)e_conv 529 if (d_conv.gt.0)then 530 write(LuOut,9040)d_conv 531 endif 532 if (g_conv.gt.0)then 533 write(LuOut,9045)g_conv 534 endif 535 call util_flush(LuOut) 536 endif 537 if (oprint_xc)then 538 write(LuOut,*) 539 call util_print_centered 540 & (LuOut,'XC Information',20,.true.) 541 if (libxcon) call nwchem_libxc_print_header() 542c 543c Write out XC info. Combo info first, than X components, 544c than C components. 545c 546 local_c = 'local ' 547 nonlocal_c = 'non-local' 548 nineb_c = ' ' 549 do n = 1, numfunc 550 if (xccomb(n))write(LuOut,9223) xcname(n) 551 enddo 552c 553c Do exact exchange differently. 554c 555 if (lxfac(1).or.nlxfac(1)) 556 & write(LuOut,9224) xname(1), xfac(1), nineb_c 557c 558 do n = 2, numfunc 559 if (lxfac(n).and.nlxfac(n))then 560 write(LuOut,9224) xname(n), xfac(n), nineb_c 561 elseif (lxfac(n).and.(.not.nlxfac(n)))then 562 write(LuOut,9224) xname(n), xfac(n), local_c 563 elseif ((.not.lxfac(n)).and.nlxfac(n))then 564 write(LuOut,9224) xname(n), xfac(n), nonlocal_c 565 endif 566 enddo 567 do n = 1, numfunc 568 if (lcfac(n).and.nlcfac(n))then 569 write(LuOut,9224) cname(n), cfac(n), nineb_c 570 elseif (lcfac(n).and.(.not.nlcfac(n)))then 571 write(LuOut,9224) cname(n), cfac(n), local_c 572 elseif ((.not.lcfac(n)).and.nlcfac(n))then 573 write(LuOut,9224) cname(n), cfac(n), nonlocal_c 574 endif 575 enddo 576 577 if (libxcon) call nwchem_libxc_print() 578c 579c Check XC coefficients to make sure appropriate components 580c sum to 1.0 581c 582c ckfac = 0.0d0 583c do n = 1, numfunc 584c if (lcfac(n))ckfac = ckfac + cfac(n) 585c enddo 586c if (abs(ckfac-1.0d0).gt.1.d-8)then 587c write(LuOut,*) 588c & ' WARNING: Sum of local correlation is ',ckfac 589c write(LuOut,*)' Sum of components do not equal unity. ' 590c endif 591c ckfac = 0.0d0 592c do n = 1, numfunc 593c if (nlcfac(n))ckfac = ckfac + cfac(n) 594c enddo 595c if (abs(ckfac-1.0d0).gt.1.d-8.and.abs(ckfac).gt.1.d-8)then 596c write(LuOut,*) 597c & ' WARNING: Sum of nonlocal correlation is ',ckfac 598c write(LuOut,*) 599c & ' Sum of components do not equal unity or 0. ' 600c endif 601c ckfac = 0.0d0 602c do n = 1, numfunc 603c if (lxfac(n))ckfac = ckfac + xfac(n) 604c enddo 605c if (abs(ckfac-1.0d0).gt.1.d-8)then 606c write(LuOut,*) 607c & ' WARNING: Sum of local exchange is ',ckfac 608c write(LuOut,*)' Sum of components do not equal unity. ' 609c endif 610c ckfac = 0.0d0 611c do n = 1, numfunc 612c if (nlxfac(n))ckfac = ckfac + xfac(n) 613c enddo 614c if (abs(ckfac-1.0d0).gt.1.d-8.and.abs(ckfac).gt.1.d-8)then 615c write(LuOut,*) 616c & ' WARNING: Sum of nonlocal exchange is ',ckfac 617c write(LuOut,*) 618c & ' Sum of components do not equal unity or 0. ' 619c endif 620c 621c Check if asymptotic correction will be added to potential 622c If both LB94 and CS00 are .true., it is assumed that the 623c user meant to use CS00 (since CS00 uses LB94) 624c 625 if (cs00) then 626 if (delta_ac.gt.1.0d90) then 627 write(LuOut,*) 628 write(LuOut,9226) 629 & ' CS with a Zhan-Nichols-Dixon shift ' 630 else 631 write(LuOut,*) 632 write(LuOut,9227) 633 & ' Casida-Salahub correction with a shift ', 634 & delta_ac,'au' 635 endif 636 else if (lb94) then 637 write(LuOut,*) 638 write(LuOut,9226) 639 & ' van Leeuwen-Baerends correction ' 640 else if (ncap) then 641 write(LuOut,*) 642 write(LuOut,9227) 643 & ' NCAP with derivative discontinuity shift ', 644 & delta_ac,'au' 645 endif 646c 647c Print range-separation parameters 648 if (cam_exch) then 649 write(LuOut,*) 650 write(LuOut,8202) 651 & 'Range-Separation Parameters ' 652 write(LuOut,8203) 653 write(LuOut,8201)cam_alpha,cam_beta,cam_omega,cam_srhf 654 end if ! cam_exch 655c 656cc AJL/Begin/FDE 657 if (frozemb_fde) then 658 write(LuOut,*) 659 call util_print_centered 660 & (LuOut,'XC Information (FDE)',23,.true.) 661c 662c Write out XC info. Combo info first, than X components, 663c than C components. 664c 665 do n = 1, numfunc 666 if (xccomb_fde(n))write(LuOut,9223) xcname(n) 667 enddo 668c 669c Do exact exchange differently. 670c AJL: Hybrids will give you linear dependencies in the fde/qm combined 671c wavefunction. Therefore, we should never have non-local contributions 672c from exact exchange. I have left the functionality here, as it is 673c lifted from L510 above, but it should never be used. 674c 675c An error message is implemented in dft_rdinput.F 676c If someone thinks it is sensible to re-enable this they can do so themselves. 677c 678 if (lxfac_fde(1).or.nlxfac_fde(1)) 679 & write(LuOut,9224) xname(1), xfac_fde(1), nineb_c 680c 681 do n = 2, numfunc 682 if (lxfac_fde(n).and.nlxfac_fde(n))then 683 write(LuOut,9224) xname(n), xfac_fde(n), nineb_c 684 elseif (lxfac_fde(n).and.(.not.nlxfac_fde(n)))then 685 write(LuOut,9224) xname(n), xfac_fde(n), local_c 686 elseif ((.not.lxfac_fde(n)).and.nlxfac_fde(n))then 687 write(LuOut,9224) xname(n), xfac_fde(n), nonlocal_c 688 endif 689 enddo 690 do n = 1, numfunc 691 if (lcfac_fde(n).and.nlcfac_fde(n))then 692 write(LuOut,9224) cname(n), cfac_fde(n), nineb_c 693 elseif (lcfac_fde(n).and.(.not.nlcfac_fde(n)))then 694 write(LuOut,9224) cname(n), cfac_fde(n), local_c 695 elseif ((.not.lcfac_fde(n)).and.nlcfac_fde(n))then 696 write(LuOut,9224) cname(n), cfac_fde(n), nonlocal_c 697 endif 698 enddo 699c 700c The most sensible thing is again to do the error checks for 701c range-separated and self-interaction corrected setups. Due to the 702c complexity of implementation I'm just going to disable this for FDE, but 703c they could in theory be checked and implemented. 704c The error check is in frozemb_fde_init, which is called from dft_scf. 705c 706ccc DISABLED ccc 707c 708cc Lifted from above: 709cc Check if asymptotic correction will be added to potential 710cc If both LB94 and CS00 are .true., it is assumed that the 711cc user meant to use CS00 (since CS00 uses LB94) 712cc 713c if (cs00) then 714c if (delta_ac.gt.1.0d90) then 715c write(LuOut,*) 716c write(LuOut,9226) 717c & ' CS with a Zhan-Nichols-Dixon shift ' 718c else 719c write(LuOut,*) 720c write(LuOut,9227) 721c & ' Casida-Salahub correction with a shift ', 722c & delta_ac,'au' 723c endif 724c else if (lb94) then 725c write(LuOut,*) 726c write(LuOut,9226) 727c & ' van Leeuwen-Baerends correction ' 728c endif 729c 730c Print range-separation parameters 731c if (cam_exch) then 732c write(LuOut,*) 733c write(LuOut,8202) 734c & 'Range-Separation Parameters ' 735c write(LuOut,8203) 736c write(LuOut,8201)cam_alpha,cam_beta,cam_omega,cam_srhf 737c end if ! cam_exch 738ccc END DISABLED ccc 739c 740c Kinetic Energy Functional 741c 742 write(LuOut,*) 743 call util_print_centered 744 & (LuOut,'KE Information (FDE)',23,.true.) 745 do n = 1, numfunc_fde 746 if (tsfac(n).gt.0) 747 & write(LuOut,9224) tsname(n), tsfac(n), local_c 748 enddo 749 endif ! frozemb_fde 750cc AJL/End 751c 752 call util_flush(LuOut) 753c 754 endif ! oprint_xc 755c 756 if (oprint_sic) then 757 if (test_sic.eq.1) then 758 write(LuOut,'(/,14x,"SIC perturbative approximation")') 759 else 760 if (test_sic.eq.2) then 761 write(LuOut, 762 . '(/14x,"SIC/OEP without localized orbitals")') 763 else 764 if (test_sic.eq.4) then 765 write(LuOut, 766 . '(/14x,"SIC/OEP with localized orbitals")') 767 end if 768 end if 769 end if 770 end if 771 if (oprint_grid.and.xc_gotxc())then 772 write(LuOut,*) 773 call util_print_centered 774 & (LuOut,'Grid Information',20,.true.) 775 write(LuOut,9135)gridtype 776 write(LuOut,9136) rgridnames(wradgrid) 777 if (.not.leb)then 778 write(LuOut,9142) 779 else 780 write(LuOut,9143) 781 endif 782 write(LuOut,9144) 783 do n = 1, ntypes 784c 785c Find an atom of this kind in the complete list. 786c 787 do ictr = 1, ncenters 788 if (iatype(ictr).eq.n) then 789 if (.not. geom_cent_get(geom, ictr, tag, 790 & ictr_coord, ictr_chg))call errquit 791 & ('dft_inpana: geom_cent_get failed', 0, 792 & GEOM_ERR) 793 goto 40 794 endif 795 enddo 796 40 continue 797 bsr = bsrad_atom_type(n)*CAU2ANG 798 radang = dble(nint(cau2ang*dble(rad_cutoff(1,n)))) 799 if (leb)then 800 write(LuOut,9138)tag,bsr,nrad(n),radang, 801 & nang(n) 802 else 803 write(LuOut,9137)tag,bsr,nrad(n),radang, 804 & nang(n),2*nang(n) 805 endif 806 enddo 807 if (no_prune)then 808 on_off_1 = 'off' 809 else 810 on_off_1 = 'on' 811 endif 812 write(LuOut,4005)on_off_1 813 write(LuOut,4004)nqshells 814c 815 if (ldelley)then 816 write(LuOut,9140) 'Delley' 817 elseif(lssw) then 818 if(whichssw.eq.'ssf ') then 819 write(LuOut,9140) ' Stratmann-Scuseria-Frisch' 820 elseif(whichssw.eq.'erf1') then 821 write(LuOut,9140) ' Erf1' 822 elseif(whichssw.eq.'erf2') then 823 write(LuOut,9140) ' Erf2' 824 else 825 write(LuOut,9140) whichssw 826 endif 827 else 828 write(LuOut,9140) 'Becke' 829 endif 830 if (nquad_task.ne.1)then 831 write(LuOut,9141)nquad_task 832 endif 833 if (nq_chunk.ne.0)then 834 write(LuOut,9145)nq_chunk 835 endif 836 call util_flush(LuOut) 837 endif 838 if (oprint_convergence)then 839 write(LuOut,*) 840 call util_print_centered 841 & (LuOut,'Convergence Information',20,.true.) 842 write(LuOut,3231)hl_tol, nfock 843 write(LuOut,3232)ndamp, rlshift 844 asap_10c = ' ASAP ' 845 start_10c = ' start ' 846 NA_10c = ' N/A ' 847 if(ncydp.ne.0)then 848 strng1 = start_10c 849 write(strng4,'(i3,7h iters )')ncydp 850 elseif(nodamping)then 851 strng1 = NA_10c 852 strng4 = NA_10c 853 else 854 write(strng1,'(d10.2)')dampon 855 write(strng4,'(d10.2)')dampoff 856 endif 857c 858 if(ncysh.ne.0)then 859 strng2 = asap_10c 860 write(strng5,'(i3,7h iters )')ncysh 861 elseif(nolevelshifting)then 862 strng2 = NA_10c 863 strng5 = NA_10c 864 else 865 write(strng2,'(d10.2)')levlon 866 write(strng5,'(d10.2)')levloff 867 endif 868c 869 if(ncyds.ne.0)then 870 strng3 = start_10c 871 write(strng6,'(i3,7h iters )')ncyds 872 elseif(nodiis)then 873 strng3 = NA_10c 874 strng6 = NA_10c 875 else 876 write(strng3,'(d10.2)')diison 877 write(strng6,'(d10.2)')diisoff 878 endif 879 write(LuOut,3233)strng1,strng2,strng3,strng4,strng5,strng6 880 call util_flush(LuOut) 881 endif 882 if(lsigma .and. oprint_general) then 883 write(luout, 884 . "(10x,'Smearing applied: ',d9.2,' (hartree)')" 885 . ) smear_sigma 886 endif 887 888 889 if (oprint_tolerances)then 890 write(LuOut,*) 891 call util_print_centered 892 & (LuOut,'Screening Tolerance Information',20,.true.) 893 t1 = 10.d0**(-itol2e) 894c write(LuOut,9372)tol_rho, iaoacc, icdacc, ixcacc, t1, 895c & r1 896 write(LuOut,9372)tol_rho, iaoacc, icdacc, ixcacc, t1 897 call util_flush(LuOut) 898 endif 899 endif 900c 901 return 902 3231 format(10x,'Convergence aids based upon iterative change in ',/, 903 & 10x,'total energy or number of iterations. ',/, 904 & 10x,'Levelshifting, if invoked, occurs when the ',/, 905 & 10x,'HOMO/LUMO gap drops below (HL_TOL): ',1Pd9.2,/, 906 & 10x,'DIIS, if invoked, will attempt to extrapolate ',/, 907 & 10x,'using up to (NFOCK): ',i2,' stored Fock matrices.',/) 908 3232 format(10x, 909 & 10x,'Damping(',i2,'%) Levelshifting(',f3.1, 910 & ') DIIS',/, 911 & 10x,8x,15('-'),1x,19('-'),1x,15('-')) 912 3233 format(10x,'dE on:',2x,a10,7x,a10,10x,a10,/, 913 & 10x,'dE off:',2x,a10,7x,a10,10x,a10,/) 914 4001 format(10x,'AO basis - number of functions: ',i5,/, 915 & 10x,' number of shells: ',i5) 916 4002 format(10x,'CD basis - number of functions: ',i5,/, 917 & 10x,' number of shells: ',i5) 918 4003 format(10x,'XC basis - number of functions: ',i5,/, 919 & 10x,' number of shells: ',i5) 920 4004 format(10x,'Number of quadrature shells: ',i5) 921 4005 format(10x,'Grid pruning is: ',a3) 922 9020 format(10x,'SCF calculation type: DFT') 923 9030 format(10x,'Maximum number of iterations: ',I3) 924 9035 format(10x,'Convergence on energy requested: ',1Pd9.2) 925 9040 format(10x,'Convergence on density requested: ',1Pd9.2) 926 9045 format(10x,'Convergence on gradient requested: ',1Pd9.2) 927 9050 format(10x,'Wavefunction type: closed shell.') 928 9052 format(10x,'Wavefunction type: restricted open shell.') 929 9055 format(10x,'Wavefunction type: spin polarized.') 930 9056 format(10x,'Use of symmetry is: ',a3, 931 & '; symmetry adaption is: ',a3) 932 9110 format(10x,'This is a Direct SCF calculation.') 933 9120 format(10x,'An Exch-Corr fitting basis will be used.') 934 9130 format(10x,'A Charge density fitting basis will be used.') 935 9135 format(10x,'Grid used for XC integration: ',a) 936cedo 9136 format(10x,'Radial quadrature: Euler-MacLaurin. ') 937 9136 format(10x,'Radial quadrature: ',A) 938 9137 format(10x,a16,2x,f6.2,6x,i3,6x,f6.1,3x,i2,1x,'*',1x,i2) 939 9138 format(10x,a16,2x,f6.2,6x,i3,8x,f6.1,3x,2x,i5) 940c 9139 format(10x,'Spatial weights used: Delley. ') 941 9140 format(10x,'Spatial weights used: ',A) 942 9141 format(10x,'Parallel task size associated with evaluation of ',/, 943 & 10x,'grid based components has been modified to: ',i2) 944 9145 format(10x,'Chunking of the angular grid is being used; ', 945 & 'nq_chunk = ',i4) 946 9142 format(10x,'Angular quadrature: Gauss-Legendre. ') 947 9143 format(10x,'Angular quadrature: Lebedev. ') 948 9144 format(10x,'Tag',14x,'B.-S. Rad.',1x,'Rad. Pts.',1x,'Rad. Cut.', 949 & 1x,'Ang. Pts.',/, 950 & 10x,'---',14x,'----------',1x,'---------',1x,'---------', 951 & 1x,'---------') 952 9223 format(10x,a40) 953 9224 format(10x,a40,1x,f6.3,1x,a9) 954 9225 format(10x,a40,1x,f8.4) 955 9226 format(10x,a44) 956 9227 format(10x,a44,1x,f10.6,1x,a2) 957 8150 format(10x,'No. of atoms :',2x,i4) 958 8200 format(10x,'No. of electrons :',2x,i4,/, 959 & 10x,' Alpha electrons :',2x,i4,/, 960 & 10x,' Beta electrons :',2x,i4,/, 961 & 10x,'Charge :',2x,i4,/, 962 & 10x,'Spin multiplicity:',2x,i4) 963 8201 format(10x,' Alpha :',f6.2,/, 964 & 10x,' Beta :',f6.2,/, 965 & 10x,' Gamma :',f6.2,/, 966 & 10x,' Short-Range HF :',l6) 967 8202 format(2x,a44) 968 8203 format(11x,'---------------------------') 969 9372 format(10x,'Density screening/tol_rho: ',1Pd9.2,/, 970 & 10x,'AO Gaussian exp screening on grid/accAOfunc: ',i3,/, 971 & 10x,'CD Gaussian exp screening on grid/accCDfunc: ',i3,/, 972 & 10x,'XC Gaussian exp screening on grid/accXCfunc: ',i3,/, 973 & 10x,'Schwarz screening/accCoul: ',1Pd9.2,/) 974 end 975 subroutine dft_inpanae(rtdb) 976 implicit none 977#include "errquit.fh" 978#include "cdft.fh" 979#include "mafdecls.fh" 980#include "global.fh" 981#include "stdio.fh" 982#include "rtdb.fh" 983#include "util.fh" 984 integer rtdb 985c 986 logical oprint_general,oprint_tolgr 987 double precision tol2e 988c 989 oprint_general = util_print('general information',print_default) 990 oprint_tolgr = util_print('grid_tol_info', print_high) 991c 992c get e_conv from rtdb because hess_init 993c 994 if (.not.rtdb_get(rtdb,'dft:e_conv',mt_dbl,1,e_conv)) 995 . call errquit('dftinpanae: rtdbget econv failed',0, RTDB_ERR) 996c 997c check integral tolerances to make sure they match 998c requested convergence tolerances. 999c 1000c make sure itol2e is less than 0.1*e_conv or (0.01*g_conv**2) 1001c 1002 if (.not. rtdb_get(rtdb, 'dft:itol2e', 1003 & mt_int, 1, itol2e)) itol2e=8 1004 if (10.d0**(-itol2e).gt.(1.0d-1*e_conv))then 1005 itol2e = -nint(log10(1.0d-1*e_conv)) 1006 if (.not. rtdb_put(rtdb, 'dft:itol2e', 1007 & mt_int, 1, itol2e)) 1008 & call errquit('dft_inpanae: rtdb_put failed', 127, RTDB_ERR) 1009 if (ga_nodeid().eq.0.and.oprint_general)then 1010 write(LuOut,*)' itol2e modified to match energy' 1011 write(LuOut,*)' convergence criterion.' 1012 endif 1013 endif 1014c 1015c check density tolerance to make sure it matches 1016c requested convergence tolerances. 1017c 1018c make sure tol_rho is less than 0.01*e_conv or (0.01*g_conv**2) 1019c 1020 if (tol_rho.gt.(1.0d-3*e_conv))then 1021 tol_rho = 1.0d-3*e_conv 1022 if (.not. rtdb_put(rtdb, 'dft:tol_rho', 1023 & mt_dbl, 1, tol_rho)) 1024 & call errquit('dft_inpanae: rtdb_put failed', 127, RTDB_ERR) 1025 if (ga_nodeid().eq.0.and.oprint_general)then 1026 write(LuOut,*)' tol_rho modified to match energy' 1027 write(LuOut,*)' convergence criterion.' 1028 endif 1029 endif 1030 if (.not. rtdb_get(rtdb, 'dft:iAOacc', mt_int, 1, 1031 & iAOacc))then 1032 iAOacc=-nint(log(e_conv)) 1033 else 1034 iAOacc=max(iAOacc,-nint(log(e_conv))) 1035 endif 1036 if (.not. rtdb_put(rtdb, 'dft:iAOacc', 1037 & mt_int, 1, iAOacc)) 1038 & call errquit('dft_inpanae: rtdb_put failed', 124, RTDB_ERR) 1039 tol2e = 10.d0**(-itol2e) 1040 if (.not. rtdb_put(rtdb, 'dft:tol2e', MT_DBL, 1, tol2e)) 1041 . call errquit('dftinpanae:rtdbput failed',0, RTDB_ERR) 1042 if (.not. rtdb_put(rtdb, 'scf:tol2e', MT_DBL, 1, tol2e)) 1043 . call errquit('dftinpanae:rtdbput failed',0, RTDB_ERR) 1044 if(oprint_tolgr.and.ga_nodeid().eq.0) then 1045 write(luout,*) ' dftinpanae: itol2e ',itol2e, 1046 , ' iaoacc ',iaoacc, 1047 , ' tol_rho ',tol_rho 1048 call util_flush(luout) 1049 endif 1050 return 1051 end 1052