1 block data cscfdata 2 implicit none 3#include "cscf.fh" 4c 5 data oinitialized /.false./ 6 data g_movecs /0,0/ 7 data k_occ /10000000/ 8 data k_eval /10000000/ 9c 10 end 11 subroutine scf_init(rtdb) 12 implicit none 13#include "errquit.fh" 14#include "mafdecls.fh" 15#include "global.fh" 16#include "bas.fh" 17#include "geom.fh" 18#include "rtdb.fh" 19#include "inp.fh" 20#include "sym.fh" 21#include "util.fh" 22#include "cscf.fh" 23#include "cosmo.fh" 24c 25 integer rtdb ! database handle 26c 27 double precision nuclear_charge 28 character*255 name 29 integer len_occ 30 external cscfdata ! For T3D linker 31c 32 logical osome 33c 34 logical hf_job 35 character*30 tag 36 character*255 theory 37 integer mult 38 logical xc_gotxc,xc_gothfx 39 external xc_gotxc,xc_gothfx 40c 41 hf_job = (.not. xc_gotxc()).and.(.not. xc_gothfx()) 42c 43 if (.not. rtdb_cget(rtdb, 'title', 1, title)) 44 $ title = ' ' 45c 46c load geometry and symmetry info 47c 48 if (.not. geom_create(geom, 'geometry')) 49 $ call errquit('scf_init: geom_create?', 0, GEOM_ERR) 50 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 51 $ call errquit('scf_init: no geometry ', 0, RTDB_ERR) 52 if (.not.rtdb_get(rtdb, 'scf:skeleton',MT_LOG, 1, oskel)) then 53 oskel = sym_number_ops(geom) .gt. 0 54 endif 55 if (.not.rtdb_get(rtdb, 'scf:adapt',MT_LOG, 1, oadapt)) then 56 oadapt = sym_number_ops(geom) .gt. 0 57 endif 58 if (.not.rtdb_get(rtdb, 'scf:lock',MT_LOG, 1, olock)) then 59 olock = .false. 60 endif 61c 62c load the basis set and get info about it 63c 64 if (.not. bas_create(basis, 'ao basis')) 65 $ call errquit('scf_init: bas_create?', 0, BASIS_ERR) 66 if (.not. bas_rtdb_load(rtdb, geom, basis, 'ao basis')) 67 $ call errquit('scf_init: no ao basis set', 0, RTDB_ERR) 68c 69c For debug ... call int_init and do the 2-e 70c 71 if (util_print('texas init debug',print_never)) then 72 call int_init(rtdb, 1, basis) 73 write(6,*) ' DONE INIT' 74* call schwarz_init(geom, basis) 75* call schwarz_tidy() 76 call int_terminate 77 endif 78c 79 if (.not. bas_name(basis, name, trans)) 80 $ call errquit('scf_init: bas_name?', 0, BASIS_ERR) 81c 82 if (.not. bas_numbf(basis, nbf)) call errquit 83 $ ('scf_init: basis info',0, BASIS_ERR) 84c 85c Is RI approximation to be used? If so get fitting basis set. 86c 87 if (rtdb_get(rtdb, 'scf:ri', MT_INT, 1, nriscf)) then 88 if (.not. bas_create(riscf_basis, 'riscf basis')) 89 $ call errquit('scf_init: bas_create?', 0, BASIS_ERR) 90 if (.not. bas_rtdb_load(rtdb, geom, riscf_basis, 'riscf basis')) 91 $ call errquit('scf_init: no riscf basis set', 0, RTDB_ERR) 92 else 93 nriscf = 0 94 riscf_basis = 0 95 endif 96c 97c Figure input/output MO vectors 98c 99 if (hf_job) then 100 tag = 'scf:input vectors' 101 else 102 tag = 'dft:input vectors' 103 endif 104 if (.not. rtdb_cget(rtdb, tag, 1, movecs_in)) 105 $ movecs_in = 'atomic' 106 if (hf_job) then 107 tag = 'scf:output vectors' 108 else 109 tag = 'dft:output vectors' 110 endif 111 if (.not. rtdb_cget(rtdb, tag, 1, movecs_out)) 112 $ movecs_out = ' ' 113 if (movecs_out.eq.' ') then 114 if (movecs_in.eq.'atomic' .or. movecs_in.eq.'hcore' .or. 115 $ movecs_in.eq.'project' .or. movecs_in.eq.'fragment'.or. 116 $ movecs_in.eq.'molden' 117 $ .or.movecs_in.eq.'rotate') then 118 call util_file_name('movecs', .false.,.false.,movecs_out) 119 else 120 movecs_out = movecs_in 121 endif 122 endif 123c 124c Resolve names of MO files to full paths defaulting to the 125c permanent directory 126c 127 if (movecs_in.eq.'atomic' .or. movecs_in.eq.'hcore' .or. 128 $ movecs_in.eq.'project' .or. movecs_in.eq.'fragment'.or. 129 $ movecs_in.eq.'molden' 130 $ .or.movecs_in.eq.'rotate') then 131 continue 132 else 133 call util_file_name_resolve(movecs_in, .false.) 134 endif 135 call util_file_name_resolve(movecs_out, .false.) 136c 137c Figure out the number of electrons from the required total 138c charge and the sum of nuclear charges 139c 140 if (.not. rtdb_get(rtdb, 'charge', MT_DBL, 1, charge)) 141 $ charge = 0.0d0 142 if (.not. geom_nuc_charge(geom, nuclear_charge)) 143 $ call errquit('scf: geom_nuc_charge failed', 0, GEOM_ERR) 144 nelec = nint(nuclear_charge - charge) 145 if (nelec .le. 0) call errquit 146 $ ('scf: negative no. of electrons ?', nelec, INPUT_ERR) 147 if (abs(nuclear_charge - charge - dble(nelec)) .gt. 1d-8) 148 $ call errquit('scf: non-integral no. of electrons ?', 0, 149 & INPUT_ERR) 150c 151c Determine no. of open and closed shells ... default is to run closed 152c shell unless told otherwise 153c 154 if(.not.rtdb_cget(rtdb,'task:theory',1,theory)) 155 + call errquit('task: no task input for theory?',0, INPUT_ERR) 156c if (.not. rtdb_get(rtdb, 'scf:nopen', MT_INT, 1, nopen)) 157c $ nopen = 0 158 if(theory .eq. 'dft')then 159 if (.not. rtdb_get(rtdb, 'dft:mult', MT_INT, 1,mult)) 160 * mult = 1 161 nopen = mult - 1 162 endif 163 if (.not. rtdb_get(rtdb, 'scf:nopen', MT_INT, 1, nopen)) 164 $ nopen = 0 165 if (nopen .gt. nelec) call errquit 166 $ ('scf: nopen > nelec ', nopen, INPUT_ERR) 167 if (mod(nelec-nopen,2) .ne. 0) call errquit 168 $ ('scf: no. of closed-shell electrons is not even!',nopen, 169 & INPUT_ERR) 170 nclosed = (nelec-nopen) / 2 171c 172 if (.not. rtdb_cget(rtdb, 'scf:scftype', 1, scftype)) then 173 if (nopen .eq. 0) then 174 scftype = 'RHF' 175 else 176 scftype = 'ROHF' 177 endif 178 endif 179c 180 call inp_ucase(scftype) 181c 182c Take care of holes in the input routines 183c 184 if (scftype.eq.'RHF' .and. nopen.gt.0) then 185 scftype = 'ROHF' 186 else if (scftype.eq.'ROHF' .and. nopen.eq.0) then 187 scftype = 'RHF' 188 endif 189c 190 if ( scftype.ne.'ROHF' .and. scftype.ne.'RHF' .and. 191 $ scftype.ne.'UHF' ) call errquit 192 $ ('scf: only ROHF, RHF, and UHF currently supported', 0, 193 & INPUT_ERR) 194c 195c Dump lagrangian? Yes by default now since if the SCF 196c has converged for an energy it will not be rerun for the gradient. 197c 198 if (.not.rtdb_get(rtdb, 'scf:lagrangian',MT_LOG, 1, olagr)) 199 $ olagr = scftype.eq.'ROHF' 200c 201 nalpha = nclosed + nopen 202 nbeta = nclosed 203c 204c DIIS toggle 205c 206 if (.not.rtdb_get(rtdb, 'scf:diis',MT_LOG, 1, odiis)) 207 $ odiis = .false. 208c 209 call ga_sync() 210c 211c For now set NMO = NBF, however this may change later when the 212c linear dependency analysis is done just before the starting guess 213c 214 nmo = nbf 215c 216c Store the derived info in the database for other wavefunction 217c modules and/or restart to access 218c 219 if (.not. rtdb_cput(rtdb, 'scf:scftype', 1, scftype)) 220 $ call errquit('scf_init: put of scftyp failed',0, RTDB_ERR) 221 if (.not. rtdb_put(rtdb, 'scf:nopen', MT_INT, 1, nopen)) 222 $ call errquit('scf_init: put of nopen failed',0, RTDB_ERR) 223 if (.not. rtdb_put(rtdb, 'scf:nclosed', MT_INT, 1, nclosed)) 224 $ call errquit('scf_init: put of nclosed failed',0, RTDB_ERR) 225 if (.not. rtdb_put(rtdb, 'scf:nelec', MT_INT, 1, nelec)) 226 $ call errquit('scf_init: put of nelec failed',0, RTDB_ERR) 227 if (.not. rtdb_put(rtdb, 'scf:nmo', MT_INT, 1, nmo)) 228 $ call errquit('scf_init: put of nmo failed',0, RTDB_ERR) 229 if (scftype .eq. 'UHF') then 230 if (.not. rtdb_put(rtdb, 'scf:nalpha', MT_INT, 1, nalpha)) 231 $ call errquit('scf_init: put of nalpha failed',0, RTDB_ERR) 232 if (.not. rtdb_put(rtdb, 'scf:nbeta', MT_INT, 1, nbeta)) 233 $ call errquit('scf_init: put of nbeta failed',0, RTDB_ERR) 234 endif 235c 236c Allocate persistent local and global arrays ... these may 237c be reallocated later when the dependency analysis is done 238c 239 if (scftype .eq. 'UHF') then 240 if (.not. ga_create(MT_DBL, nbf, nmo, 'scf_init: alpha MOs', 241 $ 32, 32, g_movecs)) call errquit('scf_init: alpha MOs', 0, 242 & GA_ERR) 243 if (.not. ga_create(MT_DBL, nbf, nmo, 'scf_init: beta MOs', 244 $ 32, 32, g_movecs(2))) call errquit('scf_init: beta MOs',0, 245 & GA_ERR) 246 else 247 if (.not. ga_create(MT_DBL, nbf, nmo, 'scf_init: MOs', 248 $ 32, 32, g_movecs)) call errquit('scf_init: MOs', 0, 249 & GA_ERR) 250 endif 251c 252 len_occ = nmo 253 if (scftype .eq. 'UHF') len_occ = nbf * 2 254 if (.not. ma_push_get(mt_dbl, len_occ, 'scf_init: mo evals', 255 $ l_eval, k_eval)) call errquit 256 $ ('scf_init: insufficient memory?', len_occ, MA_ERR) 257c 258 if (.not. ma_push_get(mt_dbl, len_occ, 'scf_init: mo occ', 259 $ l_occ, k_occ)) call errquit 260 $ ('scf_init: insufficient memory?', len_occ, MA_ERR) 261c 262 if (.not. ma_push_get(mt_int, len_occ, 'scf_init: mo irs', 263 $ l_irs, k_irs)) call errquit 264 $ ('scf_init: insufficient memory?', len_occ, MA_ERR) 265c 266 call ifill(len_occ, 1, int_mb(k_irs), 1) ! In case not adapting 267c 268c Fill in the SCF convergence info 269c 270 call scf_get_conv_info(rtdb) 271 call scf_get_fock_param(rtdb, tol2e) 272c 273c ----- cosmo initialization ---- 274c 275 cosmo_last = .false. 276 if ( rtdb_get(rtdb,'slv:cosmo',mt_log,1,cosmo_on)) then 277 if(cosmo_on) then 278 osome = util_print('information', print_low) 279c 280 call cosmo_initialize(rtdb,geom,basis,osome) 281c 282c Turn cosmo on, we want to run the calculation 283c Start with gas_phase run 284c 285 cosmo_last = .true. 286 cosmo_on = .true. 287 if(.not.rtdb_get(rtdb,'cosmo_phase',mt_int,1,cosmo_phase)) 288 > cosmo_phase = 1 289 endif 290 endif 291 oinitialized = .true. 292c 293 end 294 subroutine scf_get_conv_info(rtdb) 295C$Id$ 296 implicit none 297#include "mafdecls.fh" 298#include "global.fh" 299#include "rtdb.fh" 300#include "cscf.fh" 301c 302 integer rtdb 303c 304 double precision mone 305 parameter (mone = -1.0d0) 306c 307 ouser_changed_conv = .false. 308c 309 if (.not. rtdb_get(rtdb, 'scf:maxiter', MT_INT, 1, maxiter)) then 310 if (scftype .ne. 'UHF') then 311 maxiter = 30 312 else 313 maxiter = 30 314 endif 315 endif 316 if (.not.rtdb_get(rtdb, 'scf:thresh', MT_DBL, 1, gnorm_tol)) 317 $ gnorm_tol = 1.d-4 318c 319c Ensure that the default integral selection is sufficient 320c for the request accuracy of the SCF. However, allow user override. 321c 322 if (.not. rtdb_get(rtdb, 'scf:tol2e', MT_DBL, 1, tol2e)) 323 $ tol2e = min(1.0d-7,gnorm_tol*1d-2) 324c 325 if (rtdb_get(rtdb, 'scf:level shift info', MT_DBL, 6,shifts))then 326 ouser_changed_conv = .true. 327 else 328 call dfill(6, mone, shifts, 1) 329 endif 330 if (rtdb_get(rtdb, 'scf:full hessian switch', MT_DBL, 1, 331 $ nr_gswitch)) then 332 ouser_changed_conv = .true. 333 else 334 nr_gswitch = 0.1d0 335 endif 336c 337c Apply defaults 338c 339 if (shifts(1) .eq. -1.0d0) shifts(1) = 5.0d0 340 if (shifts(2) .eq. -1.0d0) shifts(2) = 0.5d0 341 if (shifts(3) .eq. -1.0d0) shifts(3) = 0.0d0 342 if (shifts(4) .eq. -1.0d0) shifts(4) = 0.0d0 343 if (shifts(5) .eq. -1.0d0) shifts(5) = 0.0d0 344 if (shifts(6) .eq. -1.0d0) shifts(6) = 0.0d0 345c 346 end 347 subroutine scf_get_fock_param(rtdb, tol2e) 348 implicit none 349#include "mafdecls.fh" 350#include "rtdb.fh" 351#include "cfock.fh" 352#include "global.fh" 353#include "util.fh" 354 integer rtdb 355 double precision tol2e 356c 357 logical orestart 358 logical int2e_file_open 359 external int2e_file_open 360c 361c Load balancing information 362c 363 if (.not. rtdb_get(rtdb, 'fock:task_bf', MT_INT, 1, task_bf)) 364 $ task_bf = -1 365c 366c Dentolmax constrains dentol=tol2e/denmax which is used 367c to screen integrals before screening integrals*density 368c against tol2e. Making dentolmax smaller will, up to a 369c point, make the hessian-vector products more stable, 370c though slower. In strongly convergent cases dentolmax 371c could be increased to perhaps 1d-3 with some speed gain. 372c Dentolmax is still used even if density 373c screening is turned off since otherwise the hessian 374c vector products would be disastrously expensive. 375c 376c C4H10 shows worse convergence with dentolmax > 1d-5 377c C60 shows worsse convergence with dentolmax > 1d-6 378c 379c Some poor convergence in the CPHF and line search can 380c be attributed to dentolmax being too high when preconditioning 381c with Hessian-vector products. I think there must be a problem 382c in the screening algorithm somewhere ... but where? 383c 384 if (.not. rtdb_get(rtdb, 'fock:dentolmax', mt_dbl, 1, dentolmax)) 385 $ dentolmax = 1d-6 386c 387c odensityscreen enables/disables density screening inside 388c the fock build. Turning it off can increase stability but will 389c slow things down. 390c 391 if (.not. rtdb_get(rtdb, 'fock:densityscreen', mt_log, 1, 392 $ odensityscreen)) odensityscreen = .true. 393c 394c Integral caching/file information 395c 396#if defined(CRAYXT) || defined(BGP) || defined(BGQ) || defined(NOIO) 397c direct by default to avoid io 398 if (.not. rtdb_get(rtdb, 'int2e:filesize', 399 $ MT_INT, 1, filesize)) filesize = -1 400 if (.not. rtdb_get(rtdb, 'int2e:memsize', 401 $ MT_INT, 1, memsize)) memsize = -1 402#else 403 if (.not. rtdb_get(rtdb, 'int2e:filesize', 404 $ MT_INT, 1, filesize)) filesize = 0 405 if (.not. rtdb_get(rtdb, 'int2e:memsize', 406 $ MT_INT, 1, memsize)) memsize = 0 407#endif 408 if (.not. rtdb_get(rtdb, 'int2e:restart', 409 $ MT_LOG, 1, orestart)) orestart = .false. 410c 411c The opening routine will put the .pid on the integral filename 412c (hence even tho' parallel file, open as sequential) 413c 414 if (.not. rtdb_cget(rtdb, 'int2e:filename', 1, int2efilename)) 415 $ call util_file_name('aoints',.true.,.false.,int2efilename) 416c 417 oreadfile = .false. 418 owritefile= .false. 419c 420 if (filesize.gt.0 .or. memsize.gt.0) then 421 if (.not. int2e_file_open(int2efilename, memsize, filesize, 422 $ tol2e*0.1d0, orestart)) then 423c Probably memsize/filesize too small to be useful ... just force direct 424 filesize = -1 425 memsize = -1 426 else 427c Managed to open a file or memory cache 428 if (orestart) then 429 oreadfile = .true. 430 else 431 owritefile = .true. 432 endif 433 endif 434 endif 435c 436c Blocking integral interface 437c 438 if (.not. rtdb_get(rtdb, 'fock:maxquartet', 439 $ MT_INT, 1, maxquartet)) maxquartet = 10000 440 if (.not. rtdb_get(rtdb, 'fock:maxpair', 441 $ MT_INT, 1, maxquartet)) maxpair = sqrt(dble(maxquartet)) 442 if (.not. rtdb_get(rtdb, 'fock:maxeri', 443 $ MT_INT, 1, maxeri)) maxeri = 1296*100 ! 100 d(6) quartets 444 if (.not. rtdb_get(rtdb, 'fock:maxscr', 445 $ MT_INT, 1, maxscr)) maxscr = 0 ! <=0 implies use texas estimate 446 if (.not. rtdb_get(rtdb, 'fock:intacc', 447 $ MT_DBL, 1, intacc)) intacc = 0.0d0 ! Default is variable 448c 449 if (.not. rtdb_get(rtdb, 'fock:replicated', 450 $ mt_log, 1, oreplicated)) oreplicated = .true. 451 if (.not. rtdb_get(rtdb, 'fock:noblock', 452 $ mt_log, 1, onoblock)) onoblock = .false. 453c 454c Use labels by default 455c 456 if (.not. rtdb_get( rtdb, 'fock:uselabels', MT_LOG, 1, 457 $ oerilabel)) oerilabel = .true. 458c 459 if (ga_nodeid().eq.0 .and. 460 $ util_print('fock param',print_high)) then 461 write(6,2) task_bf, maxquartet, maxeri, maxscr, intacc, 462 $ odensityscreen, dentolmax 463 2 format(/ 464 $ ' Setting fock-build task_bf :', i8/ 465 $ ' maxquartet:', i8/ 466 $ ' maxeri :', i8/ 467 $ ' maxscr :', i8/ 468 $ ' intacc :', 1p,d8.1/ 469 $ ' denscreen : ',l1/ 470 $ ' dentol :', d8.1/) 471 call util_flush(6) 472 endif 473c 474 end 475 subroutine fock_force_direct(rtdb) 476 implicit none 477#include "cfock.fh" 478 integer rtdb 479c 480c Call after scf_get_fock_param to force a fock build 481c to be direct 482c 483 call fock_2e_tidy(rtdb) 484c 485 filesize = -1 486 memsize = -1 487c 488 end 489 subroutine scf_tidy(rtdb) 490 implicit none 491#include "errquit.fh" 492#include "cscf.fh" 493#include "mafdecls.fh" 494#include "global.fh" 495#include "geom.fh" 496#include "bas.fh" 497#include "cfock.fh" 498#include "rtdb.fh" 499#include "cosmo.fh" 500 integer rtdb 501c 502 logical status 503c 504 if (oinitialized) then 505 if (.not. geom_destroy(geom)) call errquit 506 $ ('scf_tidy: geom destroy failed', 0, GEOM_ERR) 507 if (.not. bas_destroy(basis)) call errquit 508 $ ('scf_tidy: basis destroy failed',0, BASIS_ERR) 509 status = ma_pop_stack(l_irs) 510 status = ma_pop_stack(l_occ) .and. status 511 status = ma_pop_stack(l_eval) .and. status 512 if (.not. status) call errquit 513 $ ('scf_tidy: failed to free irs/occupation/evals',0, 514 & MA_ERR) 515 if (.not. ga_destroy(g_movecs)) call errquit 516 $ ('scf_tidy: failed to free movecs',0, GA_ERR) 517 if (scftype .eq. 'UHF') then 518 if (.not. ga_destroy(g_movecs(2))) call errquit 519 $ ('scf_tidy: failed to free beta movecs',0, GA_ERR) 520 endif 521 oinitialized = .false. 522 endif 523c 524 call fock_2e_tidy(rtdb) 525c 526c ----- cosmo cleanup and reset ----- 527c 528 if ( rtdb_get(rtdb,'slv:cosmo',mt_log,1,cosmo_on)) then 529 if(cosmo_on) then 530 call cosmo_tidy(rtdb) 531 cosmo_on = .false. 532 cosmo_phase = 1 533 endif 534 endif 535c 536 end 537 subroutine fock_2e_tidy(rtdb) 538 implicit none 539#include "errquit.fh" 540#include "cfock.fh" 541#include "rtdb.fh" 542#include "mafdecls.fh" 543 integer rtdb 544c 545 logical int2e_file_close 546 external int2e_file_close 547 logical okeep 548c 549 if (.not. rtdb_get(rtdb, 'int2e:keep', 550 $ MT_LOG, 1, okeep)) okeep = .false. 551 if (.not. int2e_file_close(okeep)) 552 $ call errquit('scf_tidy: closing aoints?', 0, INT_ERR) 553c 554 end 555