1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17#if defined(BUILD_PELIB) 18module pelib_interface 19 20 implicit none 21 22 private 23 24 public :: use_pelib, pelib_ifc_gspol 25 public :: pelib_ifc_do_mep, pelib_ifc_do_mep_noqm, pelib_ifc_do_cube 26 public :: pelib_ifc_do_lf 27 public :: pelib_ifc_activate, pelib_ifc_deactivate 28 public :: pelib_ifc_init, pelib_ifc_finalize, pelib_ifc_input_reader 29 public :: pelib_ifc_fock, pelib_ifc_energy, pelib_ifc_response, pelib_ifc_london 30 public :: pelib_ifc_molgrad, pelib_ifc_lf, pelib_ifc_localfield 31 public :: pelib_ifc_mep, pelib_ifc_mep_noqm, pelib_ifc_cube 32 public :: pelib_ifc_set_mixed, pelib_ifc_mixed 33 public :: pelib_ifc_do_savden, pelib_ifc_do_twoints 34 public :: pelib_ifc_save_density, pelib_ifc_twoints 35 public :: pelib_ifc_get_num_core_nuclei 36#if defined(VAR_MPI) 37 public :: pelib_ifc_slave 38#endif 39 ! TODO: update the following interface routines 40 public :: pelib_ifc_grad, pelib_ifc_lin, pelib_ifc_lr, pelib_ifc_qro 41 public :: pelib_ifc_cro 42 public :: pelib_ifc_pecc 43 public :: pelib_ifc_transformer, pelib_ifc_qrtransformer 44 45contains 46 47logical function use_pelib() 48 use pelib, only: pelib_enabled 49 if (pelib_enabled) then 50 use_pelib = .true. 51 else 52 use_pelib = .false. 53 end if 54end function use_pelib 55 56logical function pelib_ifc_gspol() 57 use pelib_options, only: pelib_gspol 58 if (pelib_gspol) then 59 pelib_ifc_gspol = .true. 60 else 61 pelib_ifc_gspol = .false. 62 end if 63end function pelib_ifc_gspol 64 65subroutine pelib_ifc_set_mixed(do_mixed) 66 use pelib_options, only: mixed 67 logical :: do_mixed 68 if (do_mixed) then 69 mixed = .true. 70 else 71 mixed = .false. 72 end if 73end subroutine pelib_ifc_set_mixed 74 75logical function pelib_ifc_do_mep() 76 use pelib_options, only: pelib_mep 77 if (pelib_mep) then 78 pelib_ifc_do_mep = .true. 79 else 80 pelib_ifc_do_mep = .false. 81 end if 82end function pelib_ifc_do_mep 83 84logical function pelib_ifc_do_mep_noqm() 85 use pelib_options, only: pelib_mep, mep_qmcube 86 if (pelib_mep .and. .not. mep_qmcube) then 87 pelib_ifc_do_mep_noqm = .true. 88 else 89 pelib_ifc_do_mep_noqm = .false. 90 end if 91end function pelib_ifc_do_mep_noqm 92 93logical function pelib_ifc_do_cube() 94 use pelib_options, only: pelib_cube 95 if (pelib_cube) then 96 pelib_ifc_do_cube = .true. 97 else 98 pelib_ifc_do_cube = .false. 99 end if 100end function pelib_ifc_do_cube 101 102logical function pelib_ifc_do_lf() 103 use pelib_options, only: pelib_lf 104 if (pelib_lf) then 105 pelib_ifc_do_lf = .true. 106 else 107 pelib_ifc_do_lf = .false. 108 end if 109end function pelib_ifc_do_lf 110 111logical function pelib_ifc_do_savden() 112 use pelib_options, only: pelib_savden 113 if (pelib_savden) then 114 pelib_ifc_do_savden = .true. 115 else 116 pelib_ifc_do_savden = .false. 117 end if 118end function pelib_ifc_do_savden 119 120logical function pelib_ifc_do_twoints() 121 use pelib_options, only: pelib_twoints 122 if (pelib_twoints) then 123 pelib_ifc_do_twoints = .true. 124 else 125 pelib_ifc_do_twoints = .false. 126 end if 127end function pelib_ifc_do_twoints 128 129subroutine pelib_ifc_activate() 130 use pelib, only: pelib_enabled 131 call qenter('pelib_ifc_activate') 132 if (use_pelib()) call quit('PElib already active') 133 pelib_enabled = .true. 134 call qexit('pelib_ifc_activate') 135end subroutine pelib_ifc_activate 136 137subroutine pelib_ifc_deactivate() 138 use pelib, only: pelib_enabled 139 call qenter('pelib_ifc_deactivate') 140 if (.not. use_pelib()) call quit('PElib already deactivated') 141 pelib_enabled = .false. 142 call qexit('pelib_ifc_deactivate') 143end subroutine pelib_ifc_deactivate 144 145subroutine pelib_ifc_input_reader(word) 146 use pelib_precision 147 use pelib_options 148 use pelib_constants 149 use pelib_utils, only: chcase 150 use pelib_cavity_generators, only: ntsatm 151#include "priunit.h" 152 character(len=*), intent(inout) :: word 153 154 character(len=80) :: option 155 character(len=2) :: auoraa 156 integer :: i, j 157 real(rp), dimension(2) :: temp 158 159 call qenter('pelib_ifc_input_reader') 160 161 potfile = 'POTENTIAL.INP' 162 h5pdefile = 'standard.h5' 163 164 do 165 read(lucmd, '(a80)') option 166 call chcase(option) 167 168 ! Read potential (optionally from potfile) 169 if (trim(option(2:7)) == 'POTENT') then 170 read(lucmd, '(a80)') option 171 backspace(lucmd) 172 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 173 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 174 read(lucmd, '(a240)') potfile 175 end if 176 ! direct solver for induced moments 177 else if (trim(option(2:7)) == 'DIRECT') then 178 read(lucmd, '(a80)') option 179 backspace(lucmd) 180 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 181 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 182 read(lucmd, '(a80)') option 183 call chcase(option) 184 if (trim(option(1:7)) == 'NOCHOL') then 185 chol = .false. 186 end if 187 end if 188 pelib_iter = .false. 189 ! iterative solver for induced moments (default) 190 else if (trim(option(2:7)) == 'ITERAT') then 191 read(lucmd, '(a80)') option 192 backspace(lucmd) 193 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 194 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 195 read(lucmd, *) thriter 196 end if 197 pelib_iter = .true. 198 else if (trim(option(2:7)) == 'DIIS T') then 199 read(lucmd, '(a80)') option 200 backspace(lucmd) 201 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 202 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 203 read(lucmd, *) thrdiis 204 end if 205 thriter = thrdiis 206 ! use reduced threshold in iterative induced moments solver 207 else if (trim(option(2:7)) == 'REDTHR') then 208 read(lucmd, '(a80)') option 209 backspace(lucmd) 210 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 211 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 212 read(lucmd, *) redlvl 213 end if 214 pelib_redthr = .true. 215 ! handling sites near quantum-classical border 216 else if (trim(option(2:7)) == 'BORDER') then 217 read(lucmd, '(a80)') option 218 backspace(lucmd) 219 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 220 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 221 read(lucmd, '(a)', advance='no') border_type 222 backspace(lucmd) 223 call chcase(border_type) 224 if ((trim(border_type) /= 'REMOVE') .and.& 225 & (trim(border_type) /= 'REDIST') .and.& 226 & (trim(border_type) /= 'CHGRED')) then 227 error stop 'unknown handling of border sites' 228 else if (trim(border_type) == 'REMOVE') then 229 read(lucmd, *) border_type, Rmin, auoraa 230 else if (trim(border_type) == 'REDIST') then 231 read(lucmd, *) border_type, redist_order, Rmin, auoraa, nredist 232 if ((nredist > 3) .or. (nredist < 1)) then 233 error stop 'parameters can only be distributed to& 234 & minimum one site and maximum three sites' 235 end if 236 else if (trim(border_type) == 'CHGRED') then 237 read(lucmd, *) border_type, Rmin, auoraa 238 else 239 error stop 'unrecognized input in .BORDER option' 240 end if 241 call chcase(auoraa) 242 if (trim(auoraa) == 'AA') Rmin = Rmin * aa2bohr 243 end if 244 pelib_border = .true. 245 ! damp electric field from induced multipoles 246 else if (trim(option(2:7)) == 'DAMP I') then 247 read(lucmd, '(a80)') option 248 backspace(lucmd) 249 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 250 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 251 read(lucmd, *) ind_damp 252 end if 253 pelib_ind_damp = .true. 254 ! damp electric field from permanent multipoles 255 else if (trim(option(2:7)) == 'DAMP M') then 256 read(lucmd, '(a80)') option 257 backspace(lucmd) 258 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 259 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 260 read(lucmd, *) mul_damp 261 end if 262 pelib_mul_damp = .true. 263 ! damp electric field from core region 264 else if (trim(option(2:7)) == 'DAMP C') then 265 read(lucmd, '(a80)') option 266 backspace(lucmd) 267 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 268 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 269 read(lucmd, *) core_damp 270 ! attempt to optionally read a custom specification of the polarizabilities 271 read(lucmd, '(a80)') option 272 backspace(lucmd) 273 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 274 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 275 read(lucmd, *) j 276 allocate(core_alphas(6,j)) 277 core_alphas = 0.0_rp 278 do i = 1, j 279 read(lucmd, *) core_alphas(1,i) 280 core_alphas(4,i) = core_alphas(1,i) 281 core_alphas(6,i) = core_alphas(1,i) 282 enddo 283 end if 284 end if 285 pelib_core_damp = .true. 286 ! damp electric fields using AMOEABA-style Thole damping 287 else if (trim(option(2:7)) == 'DAMP A') then 288 read(lucmd, '(a80)') option 289 backspace(lucmd) 290 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 291 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 292 read(lucmd, *) amoeba_damp 293 end if 294 pelib_amoeba_damp = .true. 295 ! the old deprecated DAMP option 296 else if (trim(option(2:7)) == 'DAMP') then 297 read(lucmd, '(a80)') option 298 backspace(lucmd) 299 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 300 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 301 read(lucmd, *) ind_damp 302 end if 303 pelib_ind_damp = .true. 304 write(luout, *) 'WARNING: the .DAMP option is deprecated, please use .DAMP INDUCED' 305 write(luerr, *) 'WARNING: the .DAMP option is deprecated, please use .DAMP INDUCED' 306 ! neglect dynamic response from environment 307 else if (trim(option(2:7)) == 'GSPOL') then 308 pelib_gspol = .true. 309 ! neglect many-body interactions 310 else if (trim(option(2:7)) == 'NOMB') then 311 pelib_nomb = .true. 312 ! Use existing files for restart 313 else if (trim(option(2:7)) == 'RESTAR') then 314 pelib_restart = .true. 315 ! calculate intermolecular two-electron integrals 316 else if (trim(option(2:7)) == 'TWOINT') then 317 read(lucmd, '(a80)') option 318 backspace(lucmd) 319 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 320 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 321 read(lucmd, '(a240)') h5pdefile 322 end if 323 pelib_twoints = .true. 324 ! save density matrix 325 else if (trim(option(2:7)) == 'SAVE D') then 326 read(lucmd, '(a80)') option 327 backspace(lucmd) 328 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 329 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 330 read(lucmd, '(a240)') h5pdefile 331 end if 332 pelib_savden = .true. 333 ! disable exchange repulsion 334 else if (trim(option(2:7)) == 'NO REP') then 335 pelib_repuls = .false. 336 ! electrostatics and exchange repulsion from fragment densities 337 else if (trim(option(2:7)) == 'PDE') then 338 read(lucmd, '(a80)') option 339 backspace(lucmd) 340 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 341 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 342 read(lucmd, '(a240)') h5pdefile 343 end if 344 pelib_fd = .true. 345 ! skip electrostatics from fragment densities 346 else if (trim(option(2:7)) == 'NO FD') then 347 pelib_fdes = .false. 348 ! request calculation of effective dipole integrals 349 else if (trim(option(2:7)) == 'EEF') then 350 read(lucmd, '(a80)') option 351 backspace(lucmd) 352 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 353 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 354 read(lucmd, *) ncrds 355 allocate(crds(3,ncrds)) 356 do i = 1, ncrds 357 read(lucmd, *) (crds(j,i), j = 1, 3) 358 end do 359 end if 360 pelib_lf = .true. 361 ! provide LJ parameters for the QM region 362 else if (trim(option(2:7)) == 'LJ') then 363 read(lucmd, '(a80)') option 364 backspace(lucmd) 365 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 366 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 367 read(lucmd, *) qmLJsites 368 allocate(qmLJs(2,qmLJsites)) 369 do i = 1, qmLJsites 370 read(lucmd, *) (temp(j), j = 1, 2) 371 qmLJs(1,i) = temp(1) * 2.0_rp 372 qmLJs(2,i) = temp(2) 373 end do 374 lvdw = .true. 375 end if 376 ! skip QM calculations, i.e. go directly into PE library 377 else if (trim(option(2:7)) == 'SKIPQM') then 378 pelib_skipqm = .true. 379 ! Write cube files 380 else if (trim(option(2:7)) == 'CUBE') then 381 read(lucmd, '(a80)') option 382 backspace(lucmd) 383 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 384 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 385 do 386 read(lucmd, '(a80)') option 387 call chcase(option) 388 if (trim(option(1:7)) == 'COARSE') then 389 xgrid = 3 390 ygrid = 3 391 zgrid = 3 392 else if (trim(option(1:7)) == 'MEDIUM') then 393 xgrid = 6 394 ygrid = 6 395 zgrid = 6 396 else if (trim(option(1:7)) == 'FINE') then 397 xgrid = 12 398 ygrid = 12 399 zgrid = 12 400 else if (trim(option(1:7)) == 'GRID') then 401 read(lucmd, *) xsize, xgrid, ysize, ygrid, zsize, zgrid 402 else if (trim(option(1:7)) == 'FIELD') then 403 cube_field = .true. 404 else if (option(1:1) == '.' .or. option(1:1) == '*') then 405 backspace(lucmd) 406 exit 407 else if (option(1:1) == '!' .or. option(1:1) == '#') then 408 cycle 409 else 410 error stop 'unknown input under .CUBE option' 411 end if 412 end do 413 end if 414 pelib_cube = .true. 415 ! evaluate molecular electrostatic potential 416 else if (trim(option(2:7)) == 'MEP') then 417 read(lucmd, '(a80)') option 418 backspace(lucmd) 419 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 420 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 421 do 422 read(lucmd, '(a80)') option 423 call chcase(option) 424 if (trim(option(1:7)) == 'COARSE') then 425 xgrid = 3 426 ygrid = 3 427 zgrid = 3 428 else if (trim(option(1:7)) == 'MEDIUM') then 429 xgrid = 6 430 ygrid = 6 431 zgrid = 6 432 else if (trim(option(1:7)) == 'FINE') then 433 xgrid = 12 434 ygrid = 12 435 zgrid = 12 436 else if (trim(option(1:7)) == 'GRID') then 437 read(lucmd, *) xsize, xgrid, ysize, ygrid, zsize, zgrid 438 else if (trim(option(1:7)) == 'FIELD') then 439 mep_field = .true. 440 else if (trim(option(1:7)) == 'EXTFLD') then 441 read(lucmd, *) (extfld(i), i = 1, 3) 442 mep_extfld = .true. 443 else if (trim(option(1:7)) == 'LOCFLD') then 444 read(lucmd, *) lf_component 445 mep_lf = .true. 446 else if (trim(option(1:7)) == 'SKIPQM') then 447 mep_qmcube = .false. 448 else if (trim(option(1:7)) == 'SKIPMUL') then 449 mep_mulcube = .false. 450 else if (trim(option(1:7)) == 'INPUT') then 451 mep_input = .true. 452 read(lucmd, '(a240)') h5gridfile 453 else if (option(1:1) == '.' .or. option(1:1) == '*') then 454 backspace(lucmd) 455 exit 456 else if (option(1:1) == '!' .or. option(1:1) == '#') then 457 cycle 458 else 459 error stop 'unknown option present in .MEP section.' 460 end if 461 end do 462 end if 463 pelib_mep = .true. 464 ! continuum solvation calculation 465 else if (trim(option(2:7)) == 'SOLVAT') then 466 read(lucmd, '(a80)') option 467 backspace(lucmd) 468 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 469 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 470 read(lucmd, *) solvent 471 else 472 solvent = 'H2O' 473 end if 474 pelib_sol = .true. 475 pelib_diis = .true. 476 pelib_fixsol = .true. 477 pelib_polar = .true. 478 chol = .false. 479 ! equilibrium solvation (non-equilibrium is default) 480 else if (trim(option(2:7)) == 'EQSOL') then 481 pelib_noneq = .false. 482 ! Turn off interaction between induced dipoles and surface charges 483 else if (trim(option(2:7)) == 'NOMUQ') then 484 pelib_nomuq = .true. 485 else if (trim(option(2:7)) == 'NOVMU') then 486 pelib_novmu = .true. 487 ! specify surface file 488 else if (trim(option(2:7)) == 'SURFAC') then 489 read(lucmd, '(a80)') option 490 backspace(lucmd) 491 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 492 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 493 read(lucmd, '(a240)') surfile 494 end if 495 pelib_read_surf = .true. 496 ! FixSol solvation using FIXPVA2 tesselation (optional: number of tessera per atom) 497 else if (trim(option(2:7)) == 'NTESS ') then 498 read(lucmd, '(a80)') option 499 backspace(lucmd) 500 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 501 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 502 read(lucmd, *) ntsatm 503 if ((ntsatm .ne. 60) .and. (ntsatm .ne. 240) .and. (ntsatm .ne. 960)) then 504 error stop 'number of tessera per atom must be equal to 60, 240 or 960' 505 end if 506 end if 507 ! apply external electric field 508 else if (trim(option(2:7)) == 'FIELD') then 509 pelib_field = .true. 510 read(lucmd, '(a80)') option 511 backspace(lucmd) 512 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 513 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 514 read(lucmd, *) (efields(i), i = 1, 3) 515 end if 516 ! verbose output 517 else if (trim(option(2:7)) == 'VERBOS') then 518 pelib_verbose = .true. 519 ! debug output 520 else if (trim(option(2:7)) == 'DEBUG') then 521 pelib_debug = .true. 522 pelib_verbose = .true. 523 ! isotropic polarizabilities 524 else if (trim(option(2:7)) == 'ISOPOL') then 525 pelib_isopol = .true. 526 ! zero out the polarizabilities 527 else if (trim(option(2:7)) == 'ZEROPO') then 528 pelib_zeropol = .true. 529 ! zero out higher-order multipoles 530 else if (trim(option(2:7)) == 'ZEROMU') then 531 read(lucmd, '(a80)') option 532 backspace(lucmd) 533 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 534 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 535 read(lucmd, *) zeromul_order 536 if (zeromul_order < 0) then 537 error stop 'ZEROMUL order cannot be negative' 538 end if 539 end if 540 pelib_zeromul = .true. 541 ! skip calculation of multipole-multipole interaction energy 542 else if (trim(option(2:7)) == 'SKIPMU') then 543 pelib_skipmul = .true. 544 else if (trim(option(2:7)) == 'GAUGE') then 545 read(lucmd, '(a80)') option 546 backspace(lucmd) 547 if ((option(1:1) /= '.') .and. (option(1:1) /= '*') .and.& 548 & (option(1:1) /= '!') .and. (option(1:1) /= '#')) then 549 allocate(gauge_input(3)) 550 gauge_input = 0.0_rp 551 do i = 1,3 552 read(lucmd, *) gauge_input(i) 553 enddo 554 end if 555 pelib_gauge = .true. 556 else if (option(1:1) == '*') then 557 word = option 558 exit 559 else if (option(1:1) == '!' .or. option(1:1) == '#') then 560 cycle 561 else 562 write(luout, *) 'unknown option:', option 563 error stop 'unknown option in *PELIB section' 564 end if 565 end do 566 567! check options 568 if (pelib_diis .and. .not. pelib_iter) then ! Assume that you want to use a direct solver for FixSol 569 pelib_diis = .false. 570 end if 571 572 if (pelib_sol .and. pelib_iter .and. .not. pelib_fixsol) then 573 error stop '.SOLV without .FIXSOL requires .DIRECT' 574 end if 575 576 if (pelib_mep .and. pelib_cube) then 577 error stop '.MEP and .CUBE are not compatible.' 578 end if 579 580 call qexit('pelib_ifc_input_reader') 581end subroutine pelib_ifc_input_reader 582 583subroutine pelib_ifc_init() 584 use pelib, only: pelib_init 585#include "priunit.h" 586#include "mxcent.h" 587#include "nuclei.h" 588 call qenter('pelib_ifc_init') 589 call pelib_init(lupri, cord(1:3,1:natoms), charge(1:natoms)) 590 call qexit('pelib_ifc_init') 591end subroutine pelib_ifc_init 592 593subroutine pelib_ifc_finalize() 594 use pelib, only: pelib_finalize 595 call qenter('pelib_ifc_finalize') 596 call pelib_finalize() 597 call qexit('pelib_ifc_finalize') 598end subroutine pelib_ifc_finalize 599 600subroutine pelib_ifc_fock(denmats, fckmats, energy) 601 use pelib, only: pelib_master 602#include "inforb.h" 603 real*8, dimension(nnbasx), intent(in) :: denmats 604 real*8, dimension(nnbasx), intent(out) :: fckmats 605 real*8, intent(out) :: energy 606 real*8, dimension(1) :: energies 607 call qenter('pelib_ifc_fock') 608 if (.not. use_pelib()) call quit('PElib not active') 609#if defined(VAR_MPI) 610 call pelib_ifc_start_slaves(1) 611#endif 612 call pelib_master(runtype='full_fock', & 613 triang=.true., & 614 ndim=nbast, & 615 nmats=1, & 616 denmats=denmats, & 617 fckmats=fckmats, & 618 expvals=energies) 619 energy = energies(1) 620 call qexit('pelib_ifc_fock') 621end subroutine pelib_ifc_fock 622 623subroutine pelib_ifc_mixed(denmats, fckmats, energy) 624 use pelib, only: pelib_master 625#include "inforb.h" 626 real*8, dimension(2*nnbasx), intent(in) :: denmats 627 real*8, dimension(nnbasx), intent(out) :: fckmats 628 real*8, intent(out) :: energy 629 real*8, dimension(1) :: energies 630 call qenter('pelib_ifc_mixed') 631 if (.not. use_pelib()) call quit('PElib not active') 632#if defined(VAR_MPI) 633 call pelib_ifc_start_slaves(1) 634#endif 635 call pelib_master(runtype='full_fock', & 636 triang=.true., & 637 ndim=nbast, & 638 nmats=1, & 639 denmats=denmats, & 640 fckmats=fckmats, & 641 expvals=energies) 642 energy = energies(1) 643 call qexit('pelib_ifc_mixed') 644end subroutine pelib_ifc_mixed 645 646subroutine pelib_ifc_energy(denmats, energy) 647 use pelib, only: pelib_master 648#include "inforb.h" 649 real*8, dimension(nnbasx), intent(in) :: denmats 650 real*8, intent(out), optional :: energy 651 real*8, dimension(1) :: energies 652 call qenter('pelib_ifc_energy') 653 if (.not. use_pelib()) call quit('PElib not active') 654#if defined(VAR_MPI) 655 call pelib_ifc_start_slaves(2) 656#endif 657 call pelib_master(runtype='print_energy', & 658 triang=.true., & 659 ndim=nbast, & 660 nmats=1, & 661 denmats=denmats, & 662 expvals=energies) 663 if (present(energy)) then 664 energy = energies(1) 665 end if 666 call qexit('pelib_ifc_energy') 667end subroutine pelib_ifc_energy 668 669subroutine pelib_ifc_molgrad(denmats, molgrad) 670 use pelib, only: pelib_master 671#include "inforb.h" 672#include "mxcent.h" 673#include "nuclei.h" 674 real*8, dimension(nnbasx), intent(in) :: denmats 675 real*8, dimension(3*natoms), intent(out) :: molgrad 676 call qenter('pelib_ifc_molgrad') 677 if (.not. use_pelib()) call quit('PElib not active') 678#if defined(VAR_MPI) 679 call pelib_ifc_start_slaves(5) 680#endif 681 call pelib_master(runtype='molecular_gradient', & 682 triang=.true., & 683 ndim=nbast, & 684 nmats=1, & 685 denmats=denmats, & 686 expvals=molgrad) 687 call qexit('pelib_ifc_molgrad') 688end subroutine pelib_ifc_molgrad 689 690subroutine pelib_ifc_response(denmats, fckmats, nmats) 691 use pelib, only: pelib_master 692#include "inforb.h" 693 integer, intent(in) :: nmats 694 real*8, dimension(nmats*nnbasx), intent(in) :: denmats 695 real*8, dimension(nmats*nnbasx), intent(out) :: fckmats 696 call qenter('pelib_ifc_response') 697 if (.not. use_pelib()) call quit('PElib not active') 698#if defined(VAR_MPI) 699 call pelib_ifc_start_slaves(3) 700#endif 701 call pelib_master(runtype='dynamic_response', & 702 triang=.true., & 703 ndim=nbast, & 704 nmats=nmats, & 705 denmats=denmats, & 706 fckmats=fckmats) 707 call qexit('pelib_ifc_response') 708end subroutine pelib_ifc_response 709 710subroutine pelib_ifc_london(fckmats) 711 use pelib, only: pelib_master 712#include "inforb.h" 713 real*8, dimension(3*n2basx), intent(out) :: fckmats 714 integer :: i, j, k, l, m 715 real*8, dimension(:), allocatable :: fckmats_packed 716 call qenter('pelib_ifc_london') 717 if (.not. use_pelib()) call quit('PElib not active') 718 allocate(fckmats_packed(3*nnbasx)) 719#if defined(VAR_MPI) 720 call pelib_ifc_start_slaves(4) 721#endif 722 call pelib_master('magnetic_gradient', & 723 triang=.true., & 724 ndim=nbast, & 725 fckmats=fckmats_packed) 726 do i = 1, 3 727 j = (i - 1) * nnbasx + 1 728 k = i * nnbasx 729 l = (i - 1) * n2basx + 1 730 m = i * n2basx 731 call daptge(nbas, fckmats_packed(j:k), fckmats(l:m)) 732 end do 733 deallocate(fckmats_packed) 734 call qexit('pelib_ifc_london') 735end subroutine pelib_ifc_london 736 737subroutine pelib_ifc_localfield(eefmats) 738 use pelib, only: pelib_master 739#include "inforb.h" 740 real*8, dimension(:), intent(out) :: eefmats 741 integer :: i, ndim 742 call qenter('pelib_ifc_localfield') 743 if (.not. use_pelib()) call quit('PElib not active') 744#if defined(VAR_MPI) 745 call pelib_ifc_start_slaves(8) 746#endif 747 call pelib_master(runtype='effdipole', & 748 triang=.true., & 749 ndim=nbast, & 750 fckmats=eefmats) 751 call qexit('pelib_ifc_localfield') 752end subroutine pelib_ifc_localfield 753 754subroutine pelib_ifc_lf() 755 use pelib, only: pelib_master 756#include "priunit.h" 757#include "inforb.h" 758#include "inftap.h" 759#include "orgcom.h" 760 real*8, dimension(nnbasx) :: dip 761 real*8, dimension(:),allocatable :: fckmats 762 integer :: i, j, k 763 logical :: lopen 764 character*8 :: lblinf(2) 765 allocate(fckmats(3*nnbasx)) 766 call qenter('pelib_ifc_lf') 767 if (.not. use_pelib()) call quit('PElib not active') 768 769 call flshfo(lupri) 770 lopen = .false. 771 dip = 0.0d0 772 fckmats = 0.0d0 773 774#if defined(VAR_MPI) 775 call pelib_ifc_start_slaves(8) 776#endif 777 call flshfo(lupri) 778 call pelib_master(runtype='effdipole', & 779 triang=.true., & 780 ndim=nbast, & 781 fckmats=fckmats) 782 783 if (luprop <= 0) then 784 call gpopen(luprop, 'AOPROPER', 'OLD', ' ', 'UNFORMATTED', 0, .false.) 785 lopen = .true. 786 end if 787 788! dipole integrals are stored in triangular form 789 rewind(luprop) 790 call mollb2('XDIPLEN ',lblinf,luprop,LUPRI) 791 call readt(luprop, nnbasx, dip) 792 dip(:) = dip(:) + fckmats(1:nnbasx) 793 call WRTPRO(dip,nnbasx,'XLFDIPLN',lblinf,0) 794 795 rewind(luprop) 796 call mollb2('YDIPLEN ',lblinf,luprop,LUPRI) 797 call readt(luprop, nnbasx, dip) 798 dip(:) = dip(:) + fckmats(nnbasx+1:2*nnbasx) 799 call WRTPRO(dip,nnbasx,'YLFDIPLN',lblinf,0) 800 801 rewind(luprop) 802 call mollb2('ZDIPLEN ',lblinf,luprop,LUPRI) 803 call readt(luprop, nnbasx, dip) 804 dip(:) = dip(:) + fckmats(2*nnbasx+1:3*nnbasx) 805 call WRTPRO(dip,nnbasx,'ZLFDIPLN',lblinf,0) 806 807 deallocate(fckmats) 808 if (lopen) call gpclose(luprop,'KEEP') 809 call qexit('pelib_ifc_lf') 810 811end subroutine pelib_ifc_lf 812 813subroutine pelib_ifc_mep(denmats) 814 use pelib, only: pelib_master 815 implicit none 816#include "inforb.h" 817 real*8, dimension(nnbasx), intent(in) :: denmats 818 call qenter('pelib_ifc_mep') 819#if defined(VAR_MPI) 820 call pelib_ifc_start_slaves(6) 821#endif 822 call pelib_master(runtype='mep', & 823 triang=.true., & 824 ndim=nbast, & 825 nmats=1, & 826 denmats=denmats) 827 call qexit('pelib_ifc_mep') 828end subroutine pelib_ifc_mep 829 830subroutine pelib_ifc_mep_noqm() 831 use pelib, only: pelib_master 832 implicit none 833 call qenter('pelib_ifc_mep_noqm') 834#if defined(VAR_MPI) 835 call pelib_ifc_start_slaves(6) 836#endif 837 call pelib_master(runtype='mep', & 838 triang=.true., & 839 ndim=0, & 840 nmats=0) 841 call qexit('pelib_ifc_mep_noqm') 842end subroutine pelib_ifc_mep_noqm 843 844subroutine pelib_ifc_cube(denmats, idx) 845 use pelib, only: pelib_master 846 implicit none 847#include "inforb.h" 848 real*8, dimension(nnbasx), intent(in) :: denmats 849 integer, intent(in) :: idx 850 call qenter('pelib_ifc_cube') 851#if defined(VAR_MPI) 852 call pelib_ifc_start_slaves(7) 853#endif 854 call pelib_master(runtype='cube', & 855 triang=.true., & 856 ndim=nbast, & 857 nmats=1, & 858 denmats=denmats, & 859 idx=idx) 860 call qexit('pelib_ifc_cube') 861end subroutine pelib_ifc_cube 862 863subroutine pelib_ifc_save_density(ao_denmat, mo_fckmat, mo_coefficients) 864 use pde_utils, only: pde_save_density 865#include "iprtyp.h" 866#include "maxorb.h" 867#include "infpar.h" 868#include "inforb.h" 869 real*8, dimension(:), intent(in) :: ao_denmat 870 real*8, dimension(:), intent(in) :: mo_fckmat 871 real*8, dimension(nbast,norbt), intent(in) :: mo_coefficients 872 real*8, dimension(:), allocatable :: mo_energies 873 real*8, dimension(:,:), allocatable :: ew_denmat 874 real*8, dimension(:,:), allocatable :: temp 875 integer, parameter :: iprtyp = POLARIZABLE_EMBEDDING 876 integer, parameter :: runtyp = 9 877 integer :: i 878 call qenter('pelib_ifc_save_density') 879 allocate(temp(nbast,nisht)) 880 allocate(mo_energies(nisht)) 881 do i = 1, nisht 882 mo_energies(i) = mo_fckmat(i*(i+1)/2) 883 temp(:,i) = mo_energies(i) * mo_coefficients(:,i) 884 end do 885 allocate(ew_denmat(nbast,nbast)) 886 ew_denmat = matmul(mo_coefficients(:,1:nisht), transpose(temp)) 887 deallocate(temp) 888#if defined(VAR_MPI) 889 if (nodtot >= 1) then 890 call mpixbcast(iprtyp, 1, 'INTEGER', master) 891 call mpixbcast(runtyp, 1, 'INTEGER', master) 892 end if 893#endif 894 call pde_save_density(ao_denmat, ew_denmat, nbast) 895 deallocate(ew_denmat) 896 deallocate(mo_energies) 897 call qexit('pelib_ifc_save_density') 898end subroutine pelib_ifc_save_density 899 900subroutine pelib_ifc_twoints(work, lwork) 901 use pde_utils, only: pde_twoints, pde_get_fragment_density 902#include "inforb.h" 903 real*8, dimension(:), intent(inout) :: work 904 integer, intent(in) :: lwork 905 real*8, dimension(:), allocatable :: overlap 906 real*8, dimension(:), allocatable :: core_fckmat 907 real*8, dimension(:), allocatable :: packed_frag_denmat 908 real*8, dimension(:,:), allocatable :: full_overlap 909 real*8, dimension(:,:), allocatable :: full_fckmat 910 real*8, dimension(:,:), allocatable :: full_denmat 911 real*8, dimension(:,:), allocatable :: frag_denmat 912 integer :: i, j, k 913 integer :: core_nbast 914 integer :: frag_nbast 915 integer, dimension(1) :: isymdm, ifctyp 916 call qenter('pelib_ifc_twoints') 917 call pde_get_fragment_density(packed_frag_denmat, frag_nbast) 918 allocate(frag_denmat(frag_nbast,frag_nbast)) 919 frag_denmat = 0.0d0 920 call dunfld(frag_nbast, packed_frag_denmat, frag_denmat) 921 core_nbast = nbast - frag_nbast 922 allocate(full_denmat(nbast,nbast)) 923 full_denmat = 0.0d0 924 full_denmat(core_nbast+1:nbast,core_nbast+1:nbast) = frag_denmat 925 deallocate(frag_denmat) 926 ! IFCTYP = +/-XY 927 ! X indicates symmetry about diagonal 928 ! X = 0 No symmetry 929 ! X = 1 Symmetric 930 ! X = 2 Anti-symmetric 931 ! Y indicates contributions 932 ! Y = 0 No contribution 933 ! Y = 1 Coulomb 934 ! Y = 2 Exchange 935 ! Y = 3 Coulomb + Exchange 936 ! + sign: alpha + beta matrix (singlet) 937 ! - sign: alpha - beta matrix (triplet) 938 ! SIRFCK(fckmat, denmat, ?, isymdm, ifctyp, direct, work, nwrk) 939 allocate(full_fckmat(nbast,nbast)) 940 full_fckmat = 0.0d0 941 isymdm = 1 942 ifctyp = 11 943 call sirfck(full_fckmat, full_denmat, 1, isymdm, ifctyp, .true., work(1), lwork) 944 deallocate(full_denmat) 945 allocate(core_fckmat(core_nbast*(core_nbast+1)/2)) 946 core_fckmat = 0.0d0 947 call dgetsp(core_nbast, full_fckmat(1:core_nbast,1:core_nbast), core_fckmat) 948 deallocate(full_fckmat) 949 allocate(overlap(nnbast)) 950 overlap = 0.0d0 951 call rdonel('OVERLAP', .true., overlap, nnbast) 952 allocate(full_overlap(nbast,nbast)) 953 full_overlap = 0.0d0 954 call dsptge(nbast, overlap, full_overlap) 955 deallocate(overlap) 956 call pde_twoints(core_fckmat, full_overlap(1:core_nbast,core_nbast+1:nbast), nbast) 957 deallocate(full_overlap, core_fckmat) 958 call qexit('pelib_ifc_twoints') 959end subroutine pelib_ifc_twoints 960 961integer function pelib_ifc_get_num_core_nuclei() 962 use pde_utils, only: pde_get_num_core_nuclei 963 integer :: num_nuclei 964 num_nuclei = pde_get_num_core_nuclei() 965 if (num_nuclei <= 0) then 966 call quit('Number of core nuclei must be more than zero') 967 end if 968 pelib_ifc_get_num_core_nuclei = num_nuclei 969end function pelib_ifc_get_num_core_nuclei 970 971#if defined(VAR_MPI) 972subroutine pelib_ifc_slave(runtype) 973 use pelib, only: pelib_slave 974 use pde_utils, only: pde_save_density 975 implicit none 976#include "inforb.h" 977 integer, intent(in) :: runtype 978 real*8, dimension(:), allocatable :: ao_denmat 979 real*8, dimension(:,:), allocatable :: dummy 980 call qenter('pelib_ifc_slave') 981 if (runtype == 1) then 982 call pelib_slave('full_fock') 983 else if (runtype == 2) then 984 call pelib_slave('print_energy') 985 else if (runtype == 3) then 986 call pelib_slave('dynamic_response') 987 else if (runtype == 4) then 988 call pelib_slave('magnetic_gradient') 989 else if (runtype == 5) then 990 call pelib_slave('molecular_gradient') 991 else if (runtype == 6) then 992 call pelib_slave('mep') 993 else if (runtype == 7) then 994 call pelib_slave('cube') 995 else if (runtype == 8) then 996 call pelib_slave('effdipole') 997 else if (runtype == 9) then 998 allocate(ao_denmat(nnbasx)) 999 allocate(dummy(1,1)) 1000 call pde_save_density(ao_denmat, dummy, nbast) 1001 deallocate(ao_denmat, dummy) 1002 end if 1003 call qexit('pelib_ifc_slave') 1004end subroutine pelib_ifc_slave 1005#endif 1006 1007subroutine pelib_ifc_grad(cref, cmo, cindx, dv, grd, energy, wrk, nwrk) 1008! Written by Erik Donovan Hedegård (edh) and Jogvan Magnus H. Olsen 1009! based on PCMGRAD 1010! 1011! Purpose: calculate (MCSCF) energy and gradient contribution 1012! from an embedding potential using the PE library 1013! 1014! Output: 1015! grd MCSCF gradient with PE contribution added 1016! energy total PE energy 1017! 1018! Used from common blocks: 1019! INFVAR: NCONF, NWOPT, NVAR, NVARH 1020! INFORB: NNASHX, NNBASX, NNORBX, etc. 1021! INFTAP: LUIT2 1022 implicit none 1023#include "priunit.h" 1024#include "infvar.h" 1025#include "inforb.h" 1026#include "inftap.h" 1027 1028 integer :: nwrk 1029 real*8 :: energy 1030 real*8, dimension(*) :: cref, cmo, cindx, dv, grd 1031 real*8, dimension(nwrk) :: wrk 1032 character*8 :: star8 = '********' 1033 character*8 :: solvdi = 'SOLVDIAG' 1034 character*8 :: eodata = 'EODATA ' 1035 1036 logical :: fndlab 1037 integer :: nc4, nw4, i 1038 real*8 :: solelm, ddot 1039 real*8 :: tmo, tac, test 1040 real*8, dimension(1) :: etmp 1041 real*8, dimension(:), allocatable :: fckmo, fckac 1042 real*8, dimension(:), allocatable :: pegrd, diape 1043 real*8, dimension(:), allocatable :: dcao, dvao, fdtao, fckao 1044 1045 call qenter('pelib_ifc_grad') 1046 if (.not. use_pelib()) call quit('PElib not active') 1047 1048 allocate(dcao(n2basx), dvao(n2basx)) 1049 call fckden((nisht > 0), (nasht > 0), dcao, dvao, cmo, dv, wrk, nwrk) 1050 if (nisht == 0) dcao = 0.0d0 1051 if (nasht > 0) dcao = dcao + dvao 1052 deallocate(dvao) 1053 allocate(fdtao(nnbasx)) 1054 call dgefsp(nbast, dcao, fdtao) 1055 deallocate(dcao) 1056 allocate(fckao(nnbasx)) 1057 call pelib_ifc_fock(fdtao, fckao, energy) 1058 deallocate(fdtao) 1059 allocate(fckmo(nnorbx)) 1060 call uthu(fckao, fckmo, cmo, wrk, nbast, norbt) 1061 deallocate(fckao) 1062 1063 allocate(fckac(nnashx)) 1064 fckac = 0.0d0 1065 if (nasht > 0) call getac2(fckmo, fckac) 1066 tmo = solelm(dv, fckac, fckmo, tac) 1067 1068 allocate(pegrd(nvarh)) 1069 pegrd = 0.0d0 1070 if (nconf > 1) then 1071 ! edh: SOLGC calc. < u | Fg | 0 > + < 0 | Fg | 0 > c_u 1072 call solgc(cref, fckac, tac, pegrd, cindx, wrk, nwrk) 1073 end if 1074 if (nwopt > 0) then 1075 ! edh: SOLGO calc. 2 < 0 | [Ers, Fg] | 0 > 1076 call solgo(2.0d0, dv, fckmo, pegrd(1+nconf:nvarh)) 1077 end if 1078 1079 allocate(diape(nvar)) 1080 diape = 0.0d0 1081 call soldia(tac, fckac, cindx, fckmo, dv, diape, wrk, nwrk) 1082 diape = - diape 1083 deallocate(fckmo, fckac) 1084 1085 !--------------- Orthogonality test ---------------- 1086 test = ddot(nconf, cref, 1, pegrd, 1) 1087 if (abs(test) > 1.0d-8) then 1088 nwarn = nwarn + 1 1089 write(lupri,*) ' --- PE GRADIENT WARNING --- ' 1090 write(lupri,*) ' < CREF | GRAD > =', test 1091 end if 1092 1093 ! Add PE gradient contribution to MCSCF gradient 1094 call daxpy(nvarh, 1.0d0, pegrd, 1, grd, 1) 1095 deallocate(pegrd) 1096 if (luit2 > 0) then 1097 nc4 = max(nconf, 4) 1098 nw4 = max(nwopt, 4) 1099 rewind luit2 1100 if (fndlab(eodata,luit2)) backspace luit2 1101 write(luit2) star8, star8, star8, solvdi 1102 if (nconf > 1) call writt(luit2, nc4, diape) 1103 write(luit2) star8, star8, star8, eodata 1104 end if 1105 1106 call qexit('pelib_ifc_grad') 1107 1108end subroutine pelib_ifc_grad 1109 1110subroutine pelib_ifc_lin(ncsim, nosim, bcvecs, bovecs, cref, cmo, cindx, dv, dtv,& 1111 scvecs, sovecs, orblin, wrk, nwrk) 1112! 1113! Written by Erik Donovan Hedegård and Jogvan Magnus H. Olsen 1114! after original code by Hans Joergen Aa. Jensen 1115! 1116! Common driver for pelib_lnc and pelib_lno 1117! 1118! Used from common blocks: 1119! INFLIN : NWOPPT,NVARPT 1120 implicit none 1121#include "priunit.h" 1122#include "inflin.h" 1123#include "infvar.h" 1124#include "inforb.h" 1125 1126 logical :: orblin 1127 integer :: ncsim, nosim, nwrk 1128 integer, dimension(*) :: cindx 1129 real*8, dimension(*) :: bcvecs, bovecs, scvecs, sovecs 1130 real*8, dimension(*) :: cmo, cref, dv, dtv 1131 real*8, dimension(nwrk) :: wrk 1132 1133 integer :: nso 1134 1135 call qenter('pelib_ifc_lin') 1136 if (.not. use_pelib()) call quit('PElib not active') 1137 1138 if (ncsim > 0) then 1139 call pelib_lnc(ncsim, bcvecs, cref, cmo, cindx, dv, dtv, scvecs, wrk, nwrk) 1140 end if 1141 1142 if (nosim > 0) then 1143 if (.not. orblin) then 1144 nso = nvarpt 1145 else 1146 nso = nwoppt 1147 end if 1148 call pelib_lno(nosim, bovecs, cref, cmo, cindx, dv, sovecs, nso, wrk, nwrk) 1149 end if 1150 1151 call qexit('pelib_ifc_lin') 1152 1153end subroutine pelib_ifc_lin 1154 1155subroutine pelib_lnc(ncsim, bcvecs, cref, cmo, cindx, dv, dtv, scvecs, wrk, nwrk) 1156! 1157! Written by Erik Donovan Hedegaard and Jogvan Magnus H. Olsen 1158! after original routine by Hans Jørgen Aa. Jensen 1159! 1160! Purpose: Calculate Hessian contribution from a polarizable 1161! embedding potantial to a csf trial vector. 1162! 1163! Used from common blocks: 1164! INFORB : NNASHX, NNORBX, NNBASX, etc. 1165! INFVAR : NWOPH 1166! INFLIN : NCONST, NVARPT, NWOPPT 1167! 1168 use pelib_options, only: pelib_polar 1169 implicit none 1170#include "priunit.h" 1171#include "dummy.h" 1172#include "inforb.h" 1173#include "infvar.h" 1174#include "inflin.h" 1175#include "infdim.h" 1176 1177 integer :: ncsim, nwrk 1178 integer, dimension(*) :: cindx 1179 real*8, dimension(*) :: bcvecs, cref, cmo, dv 1180 real*8, dimension(nnashx,*) :: dtv 1181 real*8, dimension(nvarpt,*) :: scvecs 1182 real*8, dimension(nwrk) :: wrk 1183 1184 logical :: fndlab 1185 integer :: i, j, jscvec, mwoph 1186 real*8 :: tfxc, tfyc, tfycac, solelm, energy 1187 real*8, dimension(:), allocatable :: dcao, dvao, fdtao, fycao 1188 real*8, dimension(:), allocatable :: dtvao, fdtvaos, fxcaos 1189 real*8, dimension(:), allocatable :: tfxcacs, fyc, fycac 1190 real*8, dimension(:,:), allocatable :: fxcs, fxcacs 1191 1192 call qenter('pelib_lnc') 1193 if (.not. use_pelib()) call quit('PElib not active') 1194 1195 allocate(fxcs(nnorbx,ncsim)) 1196 allocate(fxcacs(nnashx,ncsim)) 1197 allocate(tfxcacs(ncsim)) 1198 fxcs = 0.0d0 1199 fxcacs = 0.0d0 1200 tfxcacs = 0.0d0 1201 1202 if (pelib_polar) then 1203 ! Fxc = -R<0|Fe|B>Fe in fxcaos 1204 allocate(dtvao(n2basx)) 1205 allocate(fdtvaos(ncsim*nnbasx)) 1206 do i = 1, ncsim 1207 j = (i - 1) * nnbasx + 1 1208 call fckden(.false., .true., dummy, dtvao, cmo, dtv(:,i), wrk, nwrk) 1209 call dgefsp(nbast, dtvao, fdtvaos(j)) 1210 end do 1211 deallocate(dtvao) 1212 allocate(fxcaos(ncsim*nnbasx)) 1213 call pelib_ifc_response(fdtvaos, fxcaos, ncsim) 1214 deallocate(fdtvaos) 1215 do i = 1, ncsim 1216 j = (i - 1) * nnbasx + 1 1217 call uthu(fxcaos(j), fxcs(:,i), cmo, wrk, nbast, norbt) 1218 if (nasht > 0) call getac2(fxcs(:,i), fxcacs(:,i)) 1219 tfxc = solelm(dv, fxcacs(:,i), fxcs(:,i), tfxcacs(i)) 1220 end do 1221 deallocate(fxcaos) 1222 end if 1223 1224 ! Fg = Vmul -R<0|Fe|0>Fe in fyc 1225 allocate(dcao(n2basx), dvao(n2basx)) 1226 call fckden((nisht > 0), (nasht > 0), dcao, dvao, cmo, dv, wrk, nwrk) 1227 if (nisht == 0) dcao = 0.0d0 1228 if (nasht > 0) dcao = dcao + dvao 1229 deallocate(dvao) 1230 allocate(fdtao(nnbasx)) 1231 call dgefsp(nbast, dcao, fdtao) 1232 deallocate(dcao) 1233 allocate(fycao(nnbasx)) 1234 call pelib_ifc_fock(fdtao, fycao, energy) 1235 deallocate(fdtao) 1236 allocate(fyc(nnorbx)) 1237 call uthu(fycao, fyc, cmo, wrk, nbast, norbt) 1238 deallocate(fycao) 1239 allocate(fycac(nnashx)) 1240 if (nasht > 0) call getac2(fyc, fycac) 1241 tfyc = solelm(dv, fycac, fyc, tfycac) 1242 1243 ! ...CSF part of sigma vectors 1244 call solsc(ncsim, 0, bcvecs, cref, scvecs, fxcacs, fycac, tfxcacs, tfycac,& 1245 cindx, wrk, nwrk) 1246 1247 if (nwoppt > 0) then 1248 mwoph = nwoph 1249 nwoph = nwoppt 1250 jscvec = 1 + nconst 1251 do i = 1, ncsim 1252 if (pelib_polar) then 1253 call solgo(2.0d0, dv, fxcs(:,i), scvecs(jscvec,i)) 1254 end if 1255 call solgo(0.0d0, dtv(:,i), fyc, scvecs(jscvec,i)) 1256 end do 1257 nwoph = mwoph 1258 end if 1259 1260 deallocate(fxcacs, fycac, tfxcacs) 1261 deallocate(fyc, fxcs) 1262 1263 call qexit('pelib_lnc') 1264 1265end subroutine pelib_lnc 1266 1267subroutine pelib_lno(nosim, bovecs, cref, cmo, cindx, dv, sovecs, nso,& 1268 wrk, nwrk) 1269! 1270! Written by Erik Donovan Hedegaard and Jogvan Magnus H. Olsen 1271! after original code by Hans Jorgen Aa. Jensen 1272! 1273! Purpose: Calculate Hessian contribution from a 1274! PE potential to an orbital trial vector. 1275! 1276! NSVEC may be NVAR or NWOPT, dependent on LINTRN 1277! 1278! Used from common blocks: 1279! INFORB : NNASHX, NNORBX, NNBASX, etc. 1280! INFVAR : JWOP 1281! INFLIN : NWOPPT, NVARPT, NCONST, NCONRF 1282! 1283 use pelib_options, only: pelib_polar 1284 implicit none 1285#include "priunit.h" 1286#include "dummy.h" 1287#include "inforb.h" 1288#include "infvar.h" 1289#include "inflin.h" 1290 1291 integer :: nosim, nso, nwrk 1292 integer, dimension(*) :: cindx 1293 real*8, dimension(*) :: cref, cmo, dv 1294 real*8, dimension(nwrk) :: wrk 1295 real*8, dimension(nwoppt,*) :: bovecs 1296 real*8, dimension(nso,*) :: sovecs 1297 1298 integer :: i, j, jsovec, mwoph, ncolim 1299 logical :: fulhes, fndlab 1300 real*8 :: solelm 1301 real*8 :: txyo 1302 real*8 :: energy 1303 real*8, dimension(:), allocatable :: txyoacs 1304 real*8, dimension(:), allocatable :: ubodcao, ubodvao 1305 real*8, dimension(:), allocatable :: bodtaos, fxoaos 1306 real*8, dimension(:), allocatable :: fckmo, fyo, ufyo 1307 real*8, dimension(:), allocatable :: dcao, dvao, fdtao, fckao 1308 real*8, dimension(:,:), allocatable :: ubovecs, fxos 1309 real*8, dimension(:,:), allocatable :: fxyos, fxyoacs 1310 1311 1312 call qenter('pelib_lno') 1313 if (.not. use_pelib()) call quit('PElib not active') 1314 1315 allocate(ubovecs(n2orbx,nosim)) 1316 if (nosim > 0) then 1317 do i = 1, nosim 1318 call upkwop(nwoppt, jwop, bovecs(:,i), ubovecs(:,i)) 1319 end do 1320 end if 1321 1322 ! 1. Calculation of Fxo = R*<0|Fe(k)|O>Fe 1323 ! Store in fxos 1324 if (pelib_polar) then 1325 allocate(ubodcao(n2basx), ubodvao(n2basx)) 1326 allocate(bodtaos(nosim*nnbasx)) 1327 do i = 1, nosim 1328 j = (i - 1) * nnbasx + 1 1329 call tr1den(cmo, ubovecs(:,i), dv, ubodcao, ubodvao, wrk, nwrk) 1330 if (nasht > 0) ubodcao = ubodcao + ubodvao 1331 call dgefsp(nbast, ubodcao, bodtaos(j)) 1332 end do 1333 deallocate(ubodcao, ubodvao) 1334 allocate(fxoaos(nosim*nnbasx)) 1335 call pelib_ifc_response(bodtaos, fxoaos, nosim) 1336 deallocate(bodtaos) 1337 allocate(fxos(nnorbx,nosim)) 1338 do i = 1, nosim 1339 j = (i - 1) * nnbasx + 1 1340 call uthu(fxoaos(j), fxos(:,i), cmo, wrk, nbast, norbt) 1341 end do 1342 deallocate(fxoaos) 1343 end if 1344 1345 ! 2. Calculation of Fyo = V(k) + R<0|F|0>Fe(k) 1346 ! Store in fyos 1347 allocate(dcao(n2basx), dvao(n2basx)) 1348 call fckden((nisht > 0), (nasht > 0), dcao, dvao, cmo, dv, wrk, nwrk) 1349 if (nisht == 0) dcao = 0.0d0 1350 if (nasht > 0) dcao = dcao + dvao 1351 deallocate(dvao) 1352 allocate(fdtao(nnbasx)) 1353 call dgefsp(nbast, dcao, fdtao) 1354 deallocate(dcao) 1355 allocate(fckao(nnbasx)) 1356 call pelib_ifc_fock(fdtao, fckao, energy) 1357 deallocate(fdtao) 1358 allocate(fckmo(nnorbx)) 1359 call uthu(fckao, fckmo, cmo, wrk, nbast, norbt) 1360 deallocate(fckao) 1361 allocate(fyo(n2orbx)) 1362 call dsptsi(norbt, fckmo, fyo) 1363 deallocate(fckmo) 1364 1365 allocate(ufyo(n2orbx), txyoacs(nosim)) 1366 allocate(fxyos(nnorbx,nosim), fxyoacs(nnashx,nosim)) 1367 do i = 1, nosim 1368 ufyo = 0.0d0 1369 call tr1uh1(ubovecs(:,i), fyo, ufyo, 1) 1370 call dgetsp(norbt, ufyo, fxyos(:,i)) 1371 if (pelib_polar) then 1372 call daxpy(nnorbx, 1.0d0, fxos(:,i), 1, fxyos(:,i), 1) 1373 end if 1374 if (nasht > 0) then 1375 call getac2(fxyos(:,i), fxyoacs(:,i)) 1376 end if 1377 txyo = solelm(dv, fxyoacs(:,i), fxyos(:,i), txyoacs(i)) 1378 end do 1379 ! 3. / <0[Epq,Fxo + Fyo]|0> \ orbital part 1380 ! \ 2<0|Fyo + Fxo|mu> - <0|Fyo|0>*c0 / CSF part 1381 ! ... CSF part of sigma vectors 1382 if (lsymrf == lsymst) then 1383 ncolim = 1 1384 else 1385 ncolim = 0 1386 end if 1387 1388 ! Determine if full Hessian or only orbital Hessian 1389 fulhes = (nso == nvarpt) 1390 if (fulhes) then 1391 jsovec = 1 + nconst 1392 else 1393 jsovec = 1 1394 end if 1395 1396 if (fulhes .and. (nconst > ncolim)) then 1397 call solsc(0, nosim, dummy, cref, sovecs, fxyoacs, dummy, txyoacs,& 1398 dummy, cindx, wrk, nwrk) 1399 end if 1400 1401 ! ... orbital part of sigma vectors 1402 mwoph = nwoph 1403 nwoph = nwoppt 1404 ! ... tell SOLGO only to use the NWOPPT first JWOP entries 1405 do i = 1, nosim 1406 call solgo(2.0d0, dv, fxyos(:,i), sovecs(jsovec,i)) 1407 end do 1408 nwoph = mwoph 1409 1410 call qexit('pelib_lno') 1411 1412end subroutine pelib_lno 1413 1414subroutine pelib_ifc_lr(ncsim, nosim, bcvecs, bovecs, cref, cmo, cindx, udv,& 1415 dv, udvtr, dvtr, dtv, dtvtr, scvecs, sovecs, wrk, nwrk) 1416 implicit none 1417#include "priunit.h" 1418#include "dummy.h" 1419#include "wrkrsp.h" 1420#include "infrsp.h" 1421#include "inftap.h" 1422 1423 integer :: ncsim, nosim, nwrk 1424 integer, dimension(*) :: cindx 1425 real*8, dimension(*) :: bcvecs, bovecs 1426 real*8, dimension(*) :: cref, cmo, udv, dv 1427 real*8, dimension(*) :: udvtr, dvtr, dtv, dtvtr 1428 real*8, dimension(*) :: scvecs, sovecs 1429 real*8, dimension(nwrk) :: wrk 1430 1431 call qenter('pelib_ifc_lr') 1432 if (.not. use_pelib()) call quit('PElib not active') 1433 1434 if (ncsim > 0 .and. .not. soppa) then 1435 call pelib_rsplnc(ncsim, bcvecs, cref, cmo, cindx, udv, dv,& 1436 udvtr, dvtr, dtv, dtvtr, scvecs, wrk, nwrk) 1437 end if 1438 1439 if (nosim > 0) then 1440 call pelib_rsplno(nosim, bovecs, cref, cmo, cindx,& 1441 udv, dv, udvtr, dvtr, sovecs, wrk, nwrk) 1442 end if 1443 1444 call qexit('pelib_ifc_lr') 1445 1446end subroutine pelib_ifc_lr 1447 1448subroutine pelib_rsplnc(ncsim, bcvecs, cref, cmo, cindx, udv, dv,& 1449 udvtr, dvtr, dtv, dtvtr, scvecs, wrk, nwrk) 1450 use pelib_options, only: pelib_polar 1451 implicit none 1452#include "priunit.h" 1453#include "dummy.h" 1454#include "infrsp.h" 1455#include "inftap.h" 1456#include "wrkrsp.h" 1457#include "inforb.h" 1458#include "qrinf.h" 1459#include "infvar.h" 1460 1461 integer :: i, j 1462 integer :: ncsim, nwrk 1463 integer, dimension(*) :: cindx 1464 1465 real*8, dimension(*) :: bcvecs, cref, cmo, udv, dv 1466 real*8, dimension(*) :: udvtr, dvtr 1467 real*8, dimension(n2ashx,*) :: dtv 1468 real*8, dimension(n2ashx,*) :: dtvtr 1469 real*8, dimension(kzyvar,*) :: scvecs 1470 real*8, dimension(nwrk) :: wrk 1471 1472 real*8 :: ovlap, solelm, tfpeac, tfpe 1473 real*8, dimension(:,:), allocatable :: udtv 1474 real*8, dimension(:,:), allocatable :: fxcs, fxcacs 1475 real*8, dimension(:), allocatable :: dtvao, fuxcs 1476 real*8, dimension(:), allocatable :: fdtvaos, fxcaos 1477 real*8, dimension(:), allocatable :: tfxc, tfxcacs 1478 real*8, dimension(:), allocatable :: fpe, fupe, fpeac 1479 1480 logical :: lexist, lopen, locdeb 1481 logical :: fndlab 1482 logical :: tdm, norho2 1483 1484 call qenter('pelib_rsplnc') 1485 if (.not. use_pelib()) call quit('PElib not active') 1486 1487 locdeb = .false. 1488 1489 lopen = .false. 1490 tdm = .true. 1491 norho2 = .true. 1492 1493 allocate(fxcs(nnorbx,ncsim)) 1494 allocate(fuxcs(n2orbx)) 1495 allocate(fxcacs(nnashx,ncsim)) 1496 allocate(tfxc(ncsim)) 1497 allocate(tfxcacs(ncsim)) 1498 fxcs = 0.0d0 1499 fuxcs = 0.0d0 1500 fxcacs = 0.0d0 1501 tfxc = 0.0d0 1502 tfxcacs = 0.0d0 1503 1504 ! Fxc = R*(<0(L)|Fe|0> + <0|Fe|0(R)>)Fe 1505 if (pelib_polar .or. .not. trplet) then 1506 call getref(cref, ncref) 1507 ! ...Construct <0(L)|...|0> + <0|...|0(R)> 1508 allocate(udtv(n2ashx,ncsim)) 1509 udtv = 0.0d0 1510 call rsptdm(ncsim, irefsy, ksymst, ncref, kzconf, cref, bcvecs,& 1511 udtv, dummy, 0, 0, tdm, norho2, cindx, wrk, 1, nwrk) 1512 udtv = - 1.0d0 * udtv 1513 1514 if ( ncsim > 0 ) then 1515 allocate(fdtvaos(nnbasx*ncsim)) 1516 fdtvaos = 0.0d0 1517 allocate(dtvao(n2basx)) 1518 dtvao = 0.0d0 1519 do i = 1, ncsim 1520 j = (i - 1) * nnbasx + 1 1521 call fckden2(.false., .true., dummy, dtvao, cmo,& 1522 udtv(:,i), wrk, nwrk) 1523 call dgefsp(nbast, dtvao, fdtvaos(j)) 1524 end do 1525 deallocate(udtv, dtvao) 1526 end if 1527 1528 allocate(fxcaos(ncsim*nnbasx)) 1529 fxcaos = 0.0d0 1530 call pelib_ifc_response(fdtvaos, fxcaos, ncsim) 1531 deallocate(fdtvaos) 1532 1533 do i = 1, ncsim 1534 j = (i - 1) * nnbasx + 1 1535 call uthu(fxcaos(j), fxcs(:,i), cmo, wrk, nbast, norbt) 1536 if (nasht > 0) call getac2(fxcs(:,i), fxcacs(:,i)) 1537 tfxc = solelm(dv, fxcacs(:,i), fxcs(:,i), tfxcacs(i)) 1538 end do 1539 deallocate(fxcaos) 1540 end if 1541 1542 ! Fg = V - <0|F|0>Fe -unpack into fupe 1543 if (.not. tdhf) then 1544 allocate(fpe(nnorbx)) 1545 fpe = 0.0d0 1546 if (lusifc <= 0) then 1547 call gpopen(lusifc, 'SIRIFC', 'OLD', ' ', 'UNFORMATTED',& 1548 idummy, .false.) 1549 lopen = .true. 1550 end if 1551 rewind(lusifc) 1552 call mollab('PEFMAT ', lusifc, lupri) 1553 call readt(lusifc, nnorbx, fpe) 1554 if (lopen) call gpclose(lusifc, 'KEEP') 1555 allocate(fupe(n2orbx), fpeac(nnashx)) 1556 fupe = 0.0d0 1557 fpeac = 0.0d0 1558 call dsptsi(norbt, fpe, fupe) 1559 if (nasht > 0) call getac2(fpe, fpeac) 1560 tfpe = solelm(dv, fpeac, fpe, tfpeac) 1561 deallocate(fpe) 1562 end if 1563 1564 ! Calculate Fxc(Rxc) and Fg(Ryc) contributions to SCVECS(NVAR,NCSIM) 1565 ! ... CSF part of sigma vectors 1566 if (locdeb) then 1567 write(lupri,*)' Linear transformed configuration vector' 1568 write(lupri,*)' **** Before slvsc in pelib_rsplnc **** ' 1569 call output(scvecs,1,kzyvar,1,ncsim,kzyvar,ncsim,1,lupri) 1570 endif 1571 1572 call slvsc(ncsim, 0, nnashx, bcvecs, cref, scvecs, fxcacs,& 1573 fpeac, tfxcacs, tfpeac, cindx, wrk, nwrk) 1574 deallocate(fxcacs, tfxcacs, fpeac) 1575 1576 if (locdeb) then 1577 write(lupri,*)' Linear transformed configuration vector' 1578 write(lupri,*)' **** After slvsc in pelib_rsplnc **** ' 1579 call output(scvecs,1,kzyvar,1,ncsim,kzyvar,ncsim,1,lupri) 1580 end if 1581 1582 ! edh: The triplet will only work for response from a single reference 1583 ! state, where the term from Fxc is 0. Can be generalized (at least to 1584 ! dublet) by including Fxc... 1585 ! ... orbital part of sigma vector(s) 1586 if (kzwopt .gt. 0) then 1587 do i = 1,ncsim 1588 fuxcs = 0.0d0 1589 call dsptsi(norbt,fxcs(:,i), fuxcs) 1590 if (trplet) then 1591 ! zero for singlet reference state 1592 else 1593 call slvsor(.true.,.true., 1, udv, scvecs(1,i), fuxcs) 1594 end if 1595 if (locdeb) then 1596 write(lupri,*) ' *** After slvsor in pelib_rsplnc *** ' 1597 write(lupri,*) 'Orbital part of lin transf conf vec no ', i 1598 write(lupri,*) ' Txc contribution' 1599 call output(scvecs(1,i), 1, kzyvar, 1, 1, kzyvar, 1, 1,& 1600 lupri) 1601 end if 1602 if (trplet) then 1603 call slvsor(.false.,.false.,1, dtvtr(1,i),scvecs(1,i),fupe) 1604 else 1605 call slvsor(.false., .false., 1, dtv(1,i), scvecs(1,i), fupe) 1606 end if 1607 if (locdeb) then 1608 write(lupri,*) 'Orbital part of lin transf conf vec no ', i 1609 write(lupri,*)' Tg contribution' 1610 call output(scvecs(1,i), 1, kzyvar, 1, 1, kzyvar, 1, 1,& 1611 lupri) 1612 end if 1613 end do 1614 deallocate(fupe, fuxcs) 1615 1616 if (locdeb) then 1617 write(lupri,*)' linear transformed conf. vector' 1618 write(lupri,*)' *** after slvsor in pelib_rsplnc *** ' 1619 call output(scvecs, 1, kzyvar, 1, ncsim, kzyvar, ncsim, 1,& 1620 lupri) 1621 end if 1622 end if 1623 1624 if (ncref /= kzconf) then 1625 call quit('ERROR in pelib_rsplnc: ncref /= kzconf') 1626 end if 1627 1628 call qexit('pelib_rsplnc') 1629 1630end subroutine pelib_rsplnc 1631 1632subroutine pelib_rsplno(nosim, bovecs, cref, cmo, cindx, udv, dv,& 1633 udvtr, dvtr, sovecs, wrk, nwrk) 1634 use pelib_options, only: pelib_polar, pelib_gspol 1635 implicit none 1636#include "priunit.h" 1637#include "dummy.h" 1638#include "wrkrsp.h" 1639#include "inforb.h" 1640#include "infrsp.h" 1641#include "inftap.h" 1642 1643 integer :: nosim, nwrk 1644 integer, dimension(*) :: cindx 1645 real*8, dimension(*) :: bovecs 1646 real*8, dimension(kzyvar,*) :: sovecs 1647 real*8, dimension(*) :: cref, cmo, udv, dv, udvtr, dvtr 1648 real*8, dimension(nwrk) :: wrk 1649 1650 integer :: i, j 1651 real*8 :: txyo 1652 real*8 :: ddot, slvqlm 1653 real*8, dimension(:), allocatable :: dcao, dvao 1654 real*8, dimension(:), allocatable :: daos, fckaos 1655 real*8, dimension(:), allocatable :: daotrs 1656 real*8, dimension(:), allocatable :: evec 1657 real*8, dimension(:,:), allocatable :: ubovecs, evecs, eacs 1658 real*8, dimension(:), allocatable :: fpemo,fupemo 1659 real*8, dimension(:), allocatable :: txyoacs 1660 real*8, dimension(:), allocatable :: ovlp 1661 logical :: lexist, lopen 1662 1663 call qenter('pelib_rsplno') 1664 if (.not. use_pelib()) call quit('PElib not active') 1665 1666 ! return if no polarization and not MCSCF 1667 if (tdhf .and. .not. pelib_polar) then 1668 call qexit('pelib_rsplno') 1669 return 1670 ! no polarization for triplet excitations in closed shell SCF 1671 else if ((nasht == 0) .and. trplet) then 1672 !write(lupri,*)'WARNING: Triplet PE-response experimental' 1673 call qexit('pelib_rsplno') 1674 return 1675 ! ground state polarization approximation 1676 else if (pelib_gspol) then 1677 call qexit('pelib_rsplno') 1678 return 1679 ! triplet response for open shell (and MCSCF) not ready yet 1680 else if (tdhf .and. (nasht > 0) .and. trplet) then 1681 !write(lupri,*)'WARNING: Triplet code experimental' 1682 call quit('ERROR: triplet operators for open shell'//& 1683 ' systems not implemented') 1684 end if 1685 1686 lopen = .false. 1687 1688 if (.not. tdhf) then 1689 ! Read Fg = V - <0|F|0>Fe from file 1690 allocate(fpemo(nnorbx)) 1691 if (lusifc <= 0) then 1692 call gpopen(lusifc, 'SIRIFC', 'OLD', ' ', 'UNFORMATTED',& 1693 idummy, .false.) 1694 lopen = .true. 1695 end if 1696 rewind(lusifc) 1697 call mollab('PEFMAT ', lusifc, lupri) 1698 call readt(lusifc, nnorbx, fpemo) 1699 if (lopen) call gpclose(lusifc, 'KEEP') 1700 allocate(fupemo(n2orbx)) 1701 call dsptsi(norbt, fpemo, fupemo) 1702 deallocate(fpemo) 1703 end if 1704 1705 allocate(ubovecs(n2orbx,nosim)) 1706 call rspzym(nosim, bovecs, ubovecs) 1707 1708 ubovecs = - ubovecs 1709 1710 if (.not. trplet) then 1711 allocate(dcao(n2basx), dvao(n2basx), daos(nosim*nnbasx)) 1712 ! Calculate Fxo = <0|Fe(k)|0>Fe 1713 do i = 1, nosim 1714 j = (i - 1) * nnbasx + 1 1715 call deq27(cmo, ubovecs(:,i), udv, dcao, dvao, wrk, nwrk) 1716 if (nasht > 0) then 1717 dcao = dcao + dvao 1718 end if 1719 call dgefsp(nbast, dcao, daos(j)) 1720 end do 1721 deallocate(dcao, dvao) 1722 1723 allocate(fckaos(nosim*nnbasx)) 1724 call pelib_ifc_response(daos, fckaos, nosim) 1725 deallocate(daos) 1726 end if 1727 1728 allocate(evec(nnorbx)) 1729 allocate(evecs(n2orbx,nosim)) 1730 evecs = 0.0d0 1731 if (.not. tdhf) then 1732 allocate(eacs(n2ashx,nosim)) 1733 allocate(txyoacs(nosim)) 1734 eacs = 0.0d0 1735 txyoacs = 0.0d0 1736 end if 1737 do i = 1, nosim 1738 j = (i - 1) * nnbasx + 1 1739 if (.not. trplet) then 1740 call uthu(fckaos(j), evec, cmo, wrk, nbast, norbt) 1741 call dsptsi(norbt, evec, evecs(:,i)) 1742 end if 1743 ! Fyo = V(k) - <0|F|0>Fe(k) 1744 if (.not. tdhf) then 1745 call onexh1(ubovecs(:,i), fupemo, evecs(:,i)) 1746 call getacq(evecs(:,i), eacs(:,i)) 1747 if (trplet) then 1748 ! do nothing (txyoacs should be zero for triplet with singlet cref) 1749 else 1750 txyo = slvqlm(udv, eacs(:,i), evecs(:,i), txyoacs(i)) 1751 end if 1752 end if 1753 end do 1754 1755 deallocate(evec) 1756 if (.not. tdhf) then 1757 deallocate(fupemo) 1758 ! Special triplet handling is taken care of inside slvsc! 1759 call slvsc(0, nosim, n2ashx, dummy, cref, sovecs, eacs,& 1760 dummy, txyoacs, dummy, cindx, wrk, nwrk) 1761 deallocate(eacs) 1762 deallocate(txyoacs) 1763 end if 1764 ! Note: triplet and singlet calls to slvsor get identical. 1765 ! Singlet case: singlet evecs and singlet orbital operator. 1766 ! Triplet case: triplet evecs and triplet orbital operator. 1767 ! The two cases give the same expressions for gradient elements. 1768 call slvsor(.true., .true., nosim, udv, sovecs, evecs) 1769 deallocate(evecs) 1770 1771 call qexit('pelib_rsplno') 1772 1773end subroutine pelib_rsplno 1774 1775subroutine pelib_ifc_qro(vecb, vecc, etrs, xindx, zymb, zymc, udv, wrk, nwrk,& 1776 kzyva, kzyvb, kzyvc, isyma, isymb, isymc, cmo, mjwop) 1777 implicit none 1778#include "inforb.h" 1779#include "infvar.h" 1780#include "infrsp.h" 1781#include "infdim.h" 1782#include "qrinf.h" 1783 1784 integer :: kzyva, kzyvb, kzyvc 1785 integer :: isyma, isymb, isymc 1786 integer :: nwrk 1787 real*8, dimension(nwrk) :: wrk 1788 real*8, dimension(kzyva) :: etrs 1789 real*8, dimension(kzyvb) :: vecb 1790 real*8, dimension(kzyvc) :: vecc 1791 real*8, dimension(ncmot) :: cmo 1792 real*8, dimension(norbt,norbt) :: zymb, zymc 1793 real*8, dimension(nashdi,nashdi) :: udv 1794 real*8, dimension(lcindx) :: xindx 1795 integer, dimension(2,maxwop,8) :: mjwop 1796 1797 integer :: i, j, k 1798 integer :: idum = 1 1799 real*8, dimension(:), allocatable :: udcao, ufcmo 1800 real*8, dimension(:), allocatable :: dcaos, fcaos 1801 real*8, dimension(:), allocatable :: fcmo 1802 1803 call qenter('pelib_ifc_qro') 1804 if (.not. use_pelib()) call quit('PElib not active') 1805 1806 call gtzymt(1, vecb, kzyvb, isymb, zymb, mjwop) 1807 call gtzymt(1, vecc, kzyvc, isymc, zymc, mjwop) 1808 1809 allocate(udcao(n2basx)) 1810 allocate(ufcmo(n2orbx)) 1811 allocate(dcaos(4*nnbasx)) 1812 dcaos = 0.0d0 1813 1814 ! D(1k) 1815 udcao = 0.0d0 1816 call cdens1(isymb, cmo, zymb, udcao, wrk, nwrk) 1817 call dgefsp(nbast, udcao, dcaos(1:nnbasx)) 1818 1819 ! D(1k,2k) 1820 udcao = 0.0d0 1821 call cdens2(isymb, isymc, cmo, zymb, zymc, udcao,& 1822 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1823 call dgefsp(nbast, udcao, dcaos(nnbasx+1:2*nnbasx)) 1824 1825 ! D(2k) 1826 udcao = 0.0d0 1827 call cdens1(isymc, cmo, zymc, udcao, wrk, nwrk) 1828 call dgefsp(nbast, udcao, dcaos(2*nnbasx+1:3*nnbasx)) 1829 1830 ! D(2k,1k) 1831 udcao = 0.0d0 1832 call cdens2(isymc, isymb, cmo, zymc, zymb, udcao,& 1833 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1834 call dgefsp(nbast, udcao, dcaos(3*nnbasx+1:4*nnbasx)) 1835 1836 deallocate(udcao) 1837 1838 allocate(fcaos(4*nnbasx)) 1839 call pelib_ifc_response(dcaos, fcaos, 4) 1840 deallocate(dcaos) 1841 1842 allocate(fcmo(nnorbx)) 1843 ufcmo = 0.0d0 1844 1845 i = 1 1846 j = nnbasx 1847 call uthu(2.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 1848 wrk(1:n2orbx) = 0.0d0 1849 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 1850 call oith1(isymc, zymc, wrk(1:n2orbx), ufcmo, isyma) 1851 1852 i = i + nnbasx 1853 j = j + nnbasx 1854 call uthu(fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 1855 wrk(1:n2orbx) = 0.0d0 1856 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 1857 ufcmo = ufcmo + wrk(1:n2orbx) 1858 1859 i = i + nnbasx 1860 j = j + nnbasx 1861 call uthu(2.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 1862 wrk(1:n2orbx) = 0.0d0 1863 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 1864 call oith1(isymb, zymb, wrk(1:n2orbx), ufcmo, isyma) 1865 1866 i = i + nnbasx 1867 j = j + nnbasx 1868 call uthu(fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 1869 wrk(1:n2orbx) = 0.0d0 1870 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 1871 ufcmo = ufcmo + wrk(1:n2orbx) 1872 1873 call rsp1gr(1, kzyva, idum, 0, isyma, 0, 1, etrs,& 1874 wrk, idum, idum, 1.0d0, 1, udv, ufcmo, xindx,& 1875 mjwop, wrk, nwrk, .true., .false., .false.) 1876 1877 deallocate(fcaos, fcmo, ufcmo) 1878 1879 call qexit('pelib_ifc_qro') 1880 1881end subroutine pelib_ifc_qro 1882 1883subroutine pelib_ifc_cro(vecb, vecc, vecd, etrs, xindx, zymb, zymc, zymd, udv,& 1884 wrk, nwrk, kzyva, kzyvb, kzyvc, kzyvd, isyma, isymb,& 1885 isymc, isymd, cmo,mjwop) 1886 implicit none 1887#include "inforb.h" 1888#include "infvar.h" 1889#include "infrsp.h" 1890#include "infdim.h" 1891#include "qrinf.h" 1892 1893 integer :: kzyva, kzyvb, kzyvc, kzyvd 1894 integer :: isyma, isymb, isymc, isymd 1895 integer :: nwrk 1896 real*8, dimension(nwrk) :: wrk 1897 real*8, dimension(kzyva) :: etrs 1898 real*8, dimension(kzyvb) :: vecb 1899 real*8, dimension(kzyvc) :: vecc 1900 real*8, dimension(kzyvd) :: vecd 1901 real*8, dimension(ncmot) :: cmo 1902 real*8, dimension(norbt,norbt) :: zymb, zymc, zymd 1903 real*8, dimension(nashdi,nashdi) :: udv 1904 real*8, dimension(lcindx) :: xindx 1905 integer, dimension(2,maxwop,8) :: mjwop 1906 1907 integer :: i, j, k 1908 integer :: idum = 1 1909 real*8, dimension(:), allocatable :: udcao, ufcmo 1910 real*8, dimension(:), allocatable :: dcaos, fcaos 1911 real*8, dimension(:), allocatable :: fcmo 1912 1913 call qenter('pelib_ifc_cro') 1914 if (.not. use_pelib()) call quit('PElib not active') 1915 1916 call gtzymt(1, vecb, kzyvb, isymb, zymb, mjwop) 1917 call gtzymt(1, vecc, kzyvc, isymc, zymc, mjwop) 1918 call gtzymt(1, vecd, kzyvd, isymd, zymd, mjwop) 1919 1920 allocate(udcao(n2basx)) 1921 allocate(ufcmo(n2orbx)) 1922 allocate(dcaos(15*nnbasx)) 1923 dcaos = 0.0d0 1924 1925 udcao = 0.0d0 1926 call cdens1(isymb, cmo, zymb, udcao, wrk, nwrk) 1927 call dgefsp(nbast, udcao, dcaos(1:nnbasx)) 1928 udcao = 0.0d0 1929 call cdens1(isymc, cmo, zymc, udcao, wrk, nwrk) 1930 call dgefsp(nbast, udcao, dcaos(nnbasx+1:2*nnbasx)) 1931 udcao = 0.0d0 1932 call cdens1(isymd, cmo, zymd, udcao, wrk, nwrk) 1933 call dgefsp(nbast, udcao, dcaos(2*nnbasx+1:3*nnbasx)) 1934 1935 udcao = 0.0d0 1936 call cdens2(isymb, isymc, cmo, zymb, zymc, udcao,& 1937 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1938 call dgefsp(nbast, udcao, dcaos(3*nnbasx+1:4*nnbasx)) 1939 udcao = 0.0d0 1940 call cdens2(isymc, isymb, cmo, zymc, zymb, udcao,& 1941 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1942 call dgefsp(nbast, udcao, dcaos(4*nnbasx+1:5*nnbasx)) 1943 udcao = 0.0d0 1944 call cdens2(isymb, isymd, cmo, zymb, zymd, udcao,& 1945 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1946 call dgefsp(nbast, udcao, dcaos(5*nnbasx+1:6*nnbasx)) 1947 udcao = 0.0d0 1948 call cdens2(isymd, isymb, cmo, zymd, zymb, udcao,& 1949 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1950 call dgefsp(nbast, udcao, dcaos(6*nnbasx+1:7*nnbasx)) 1951 udcao = 0.0d0 1952 call cdens2(isymc, isymd, cmo, zymc, zymd, udcao,& 1953 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1954 call dgefsp(nbast, udcao, dcaos(7*nnbasx+1:8*nnbasx)) 1955 udcao = 0.0d0 1956 call cdens2(isymd, isymc, cmo, zymd, zymc, udcao,& 1957 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1958 call dgefsp(nbast, udcao, dcaos(8*nnbasx+1:9*nnbasx)) 1959 1960 udcao = 0.0d0 1961 call cdens3(isymb, isymc, isymd, cmo, zymb, zymc, zymd, udcao,& 1962 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1963 call dgefsp(nbast, udcao, dcaos(9*nnbasx+1:10*nnbasx)) 1964 udcao = 0.0d0 1965 call cdens3(isymd, isymb, isymc, cmo, zymd, zymb, zymc, udcao,& 1966 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1967 call dgefsp(nbast, udcao, dcaos(10*nnbasx+1:11*nnbasx)) 1968 udcao = 0.0d0 1969 call cdens3(isymc, isymd, isymb, cmo, zymc, zymd, zymb, udcao,& 1970 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1971 call dgefsp(nbast, udcao, dcaos(11*nnbasx+1:12*nnbasx)) 1972 udcao = 0.0d0 1973 call cdens3(isymb, isymd, isymc, cmo, zymb, zymd, zymc, udcao,& 1974 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1975 call dgefsp(nbast, udcao, dcaos(12*nnbasx+1:13*nnbasx)) 1976 udcao = 0.0d0 1977 call cdens3(isymc, isymb, isymd, cmo, zymc, zymb, zymd, udcao,& 1978 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1979 call dgefsp(nbast, udcao, dcaos(13*nnbasx+1:14*nnbasx)) 1980 udcao = 0.0d0 1981 call cdens3(isymd, isymc, isymb, cmo, zymd, zymc, zymb, udcao,& 1982 wrk(1:n2basx), wrk(n2basx+1:2*n2basx), ufcmo) 1983 call dgefsp(nbast, udcao, dcaos(14*nnbasx+1:15*nnbasx)) 1984 1985 deallocate(udcao) 1986 1987 allocate(fcaos(15*nnbasx)) 1988 call pelib_ifc_response(dcaos, fcaos, 15) 1989 deallocate(dcaos) 1990 1991 allocate(fcmo(nnorbx)) 1992 ufcmo = 0.0d0 1993 1994 i = 1 1995 j = nnbasx 1996 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 1997 wrk(1:2*n2orbx) = 0.0d0 1998 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 1999 call oith1(isymc, zymc, wrk(1:n2orbx), wrk(n2orbx+1:2*n2orbx), isyma) 2000 call oith1(isymd, zymd, wrk(n2orbx+1:2*n2orbx), ufcmo, isyma) 2001 i = i + nnbasx 2002 j = j + nnbasx 2003 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2004 wrk(1:2*n2orbx) = 0.0d0 2005 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2006 call oith1(isymb, zymb, wrk(1:n2orbx), wrk(n2orbx+1:2*n2orbx), isyma) 2007 call oith1(isymd, zymd, wrk(n2orbx+1:2*n2orbx), ufcmo, isyma) 2008 i = i + nnbasx 2009 j = j + nnbasx 2010 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2011 wrk(1:2*n2orbx) = 0.0d0 2012 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2013 call oith1(isymb, zymb, wrk(1:n2orbx), wrk(n2orbx+1:2*n2orbx), isyma) 2014 call oith1(isymc, zymc, wrk(n2orbx+1:2*n2orbx), ufcmo, isyma) 2015 i = 1 2016 j = nnbasx 2017 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2018 wrk(1:2*n2orbx) = 0.0d0 2019 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2020 call oith1(isymd, zymd, wrk(1:n2orbx), wrk(n2orbx+1:2*n2orbx), isyma) 2021 call oith1(isymc, zymc, wrk(n2orbx+1:2*n2orbx), ufcmo, isyma) 2022 i = i + nnbasx 2023 j = j + nnbasx 2024 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2025 wrk(1:2*n2orbx) = 0.0d0 2026 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2027 call oith1(isymd, zymd, wrk(1:n2orbx), wrk(n2orbx+1:2*n2orbx), isyma) 2028 call oith1(isymb, zymb, wrk(n2orbx+1:2*n2orbx), ufcmo, isyma) 2029 i = i + nnbasx 2030 j = j + nnbasx 2031 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2032 wrk(1:2*n2orbx) = 0.0d0 2033 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2034 call oith1(isymc, zymc, wrk(1:n2orbx), wrk(n2orbx+1:2*n2orbx), isyma) 2035 call oith1(isymb, zymb, wrk(n2orbx+1:2*n2orbx), ufcmo, isyma) 2036 2037 i = i + nnbasx 2038 j = j + nnbasx 2039 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2040 wrk(1:n2orbx) = 0.0d0 2041 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2042 call oith1(isymd, zymd, wrk(1:n2orbx), ufcmo, isyma) 2043 i = i + nnbasx 2044 j = j + nnbasx 2045 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2046 wrk(1:n2orbx) = 0.0d0 2047 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2048 call oith1(isymd, zymd, wrk(1:n2orbx), ufcmo, isyma) 2049 i = i + nnbasx 2050 j = j + nnbasx 2051 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2052 wrk(1:n2orbx) = 0.0d0 2053 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2054 call oith1(isymc, zymc, wrk(1:n2orbx), ufcmo, isyma) 2055 i = i + nnbasx 2056 j = j + nnbasx 2057 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2058 wrk(1:n2orbx) = 0.0d0 2059 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2060 call oith1(isymc, zymc, wrk(1:n2orbx), ufcmo, isyma) 2061 i = i + nnbasx 2062 j = j + nnbasx 2063 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2064 wrk(1:n2orbx) = 0.0d0 2065 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2066 call oith1(isymb, zymb, wrk(1:n2orbx), ufcmo, isyma) 2067 i = i + nnbasx 2068 j = j + nnbasx 2069 call uthu(1.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2070 wrk(1:n2orbx) = 0.0d0 2071 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2072 call oith1(isymb, zymb, wrk(1:n2orbx), ufcmo, isyma) 2073 2074 i = i + nnbasx 2075 j = j + nnbasx 2076 call uthu(1.0d0/3.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2077 wrk(1:n2orbx) = 0.0d0 2078 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2079 ufcmo = ufcmo + wrk(1:n2orbx) 2080 i = i + nnbasx 2081 j = j + nnbasx 2082 call uthu(1.0d0/3.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2083 wrk(1:n2orbx) = 0.0d0 2084 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2085 ufcmo = ufcmo + wrk(1:n2orbx) 2086 i = i + nnbasx 2087 j = j + nnbasx 2088 call uthu(1.0d0/3.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2089 wrk(1:n2orbx) = 0.0d0 2090 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2091 ufcmo = ufcmo + wrk(1:n2orbx) 2092 i = i + nnbasx 2093 j = j + nnbasx 2094 call uthu(1.0d0/3.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2095 wrk(1:n2orbx) = 0.0d0 2096 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2097 ufcmo = ufcmo + wrk(1:n2orbx) 2098 i = i + nnbasx 2099 j = j + nnbasx 2100 call uthu(1.0d0/3.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2101 wrk(1:n2orbx) = 0.0d0 2102 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2103 ufcmo = ufcmo + wrk(1:n2orbx) 2104 i = i + nnbasx 2105 j = j + nnbasx 2106 call uthu(1.0d0/3.0d0*fcaos(i:j), fcmo, cmo, wrk, nbast, norbt) 2107 wrk(1:n2orbx) = 0.0d0 2108 call dsptsi(norbt, fcmo, wrk(1:n2orbx)) 2109 ufcmo = ufcmo + wrk(1:n2orbx) 2110 2111 call rsp1gr(1, kzyva, idum, 0, isyma, 0, 1, etrs,& 2112 wrk, idum, idum, 1.0d0, 1, udv, ufcmo, xindx,& 2113 mjwop, wrk, nwrk, .true., .false., .false.) 2114 2115 deallocate(fcaos, fcmo, ufcmo) 2116 2117 call qexit('pelib_ifc_cro') 2118 2119end subroutine pelib_ifc_cro 2120 2121subroutine pelib_ifc_pecc(aoden, aodencc, converged, t_or_tbar) 2122 implicit none 2123#include "priunit.h" 2124#include "ccsdinp.h" 2125#include "ccsdsym.h" 2126#include "ccorb.h" 2127#include "qm3.h" 2128#include "ccfdgeo.h" 2129#include "ccslvinf.h" 2130#include "ccinftap.h" 2131 2132 real*8, parameter :: zero = 0.0d0 2133 integer, intent(in) :: t_or_tbar 2134 logical :: converged 2135 real*8, intent(in) :: aoden(n2bst(isymop)), aodencc(n2bst(isymop)) 2136 real*8 :: energy, ecmcu, eccgrs1 2137 real*8, allocatable :: denmat(:), fockmat(:) 2138 character*10 :: model, label 2139 character*6 :: modelpri 2140 2141 call qenter('PELIB_IFC_PECC') 2142 2143 write (lupri,'(1X,A,/)') ' ' 2144 write (lupri,'(1X,A,/)') & 2145 '******************************************************************' 2146 write (lupri,'(1X,A,/)') & 2147 '**** Output from Polarizable Embedding Coupled Cluster module ****' 2148 write (lupri,'(1X,A,/)') & 2149 '******************************************************************' 2150 2151 model = 'CCSD' 2152 if (ccsd) then 2153 call around('The Coupled Cluster model is CCSD') 2154 model = 'CCSD' 2155 modelpri = 'CCSD' 2156 end if 2157 if ((cc2) .and. (.not. mcc2)) then 2158 call around('The Coupled Cluster model is CC2') 2159 model = 'CC2' 2160 modelpri = ' CC2' 2161 end if 2162 if (.not. (ccsd .or. cc2)) then 2163 call quit('PECC only implemented for CCSD and CC2 parametrizations') 2164 end if 2165 allocate(denmat(nnbasx), fockmat(nnbasx)) 2166 fockmat = 0.0d0 2167 call dgefsp(nbas, aoden, denmat) 2168 if (hffld) then 2169 call pelib_ifc_energy(denmat, energy) 2170 else 2171 call pelib_ifc_fock(denmat, fockmat, energy) 2172 call pelib_ifc_energy(denmat) 2173 call around('Writing the Fock matrix to FOCKMAT file') 2174 call put_to_file('FOCKMAT', nnbasx, fockmat) 2175 end if 2176 deallocate(denmat, fockmat) 2177 ecmcu = eccgrs + energy 2178 if (abs(ecmcu-eccpr).lt.cvgesol) lslecvg = .true. 2179 write(lupri,*) 'E(PECC) contribution in iteration', iccslit, ': ', energy 2180 eccpr = ecmcu 2181 converged = .false. 2182 if (lslecvg .and. lsltcvg .and. lsllcvg) then 2183 converged = .true. 2184 if (loiter) then 2185 write(lupri,'(1X,A,I3,A)') 'Maximum inner iterations for '// & 2186 't set to', mxtinit, 'in each outer iteration.' 2187 write(lupri,'(1X,A,I3,A)') 'Maximum inner iterations for '//& 2188 't-bar set to', mxlinit, 'in each outer iteration.' 2189 end if 2190 write(lupri,*) 'PECC equations are converged in ', iccslit, & 2191 ' outer iterations' 2192 write(lupri,*) 'PECC equations are converged in ', nslvinit, & 2193 ' inner iterations' 2194 write(lures,'(12X,A4,A,F20.10)') modelpri, ' Total energy: ',& 2195 ecmcu 2196 write(lures,'(12X,A4,A,F20.10)') modelpri, ' E(PE-CC) : ',& 2197 energy 2198 eccgrs1 = ecmcu 2199 label = 'PE-'//modelpri//' ' 2200 call wripro(eccgrs1,model,0,label,label,label,label,zero,zero,zero,1,0,0,0) 2201 label = 'E(PE-CC) ' 2202 call wripro(energy,model,0,label,label,label,label,zero,zero,zero,1,0,0,0) 2203 else 2204 iccslit = iccslit + 1 2205 if (iccslit .gt. mxccslit) then 2206 write(lupri,*) 'Maximum number of PE-CC iterations:', mxccslit, & 2207 'is reached!' 2208 call quit('Maximum number of PE-CC iterations reached') 2209 end if 2210 end if 2211 2212 call qexit('PELIB_IFC_PECC') 2213 2214end subroutine pelib_ifc_pecc 2215 2216subroutine pelib_ifc_transformer(rho1,rho2,ctr1,ctr2,model,isymtr,lr,work,lwork) 2217 implicit none 2218 2219#include "mxcent.h" 2220#include "qmmm.h" 2221#include "qm3.h" 2222#include "ccsdsym.h" 2223#include "priunit.h" 2224#include "ccsdinp.h" 2225#include "ccslvinf.h" 2226#include "ccorb.h" 2227 2228 real*8, parameter :: zero = 0.0d0, one = 1.0d0, half = 0.5d0, two = 2.0d0 2229 integer, parameter :: izero = 0 2230 integer :: lwork, ndim, idldum, isydum, idlino, isymtr 2231 character*1, intent(in) :: lr 2232 character*2 :: list 2233 real*8, dimension(lwork) :: work 2234 real*8 :: rho1(nt1am(isymtr)), rho2(nt2am(isymtr)), ctr1(nt1am(isymtr)),& 2235 ctr2(nt2am(isymtr)), ddot, rho1n, rho2n, tal1, tal2, fact 2236 real*8, allocatable :: denmats(:), gmat(:), eta(:), denmattemp(:), gmattemp(:) 2237 character*2 :: lisdum 2238 character*10 :: model 2239 logical, parameter :: locdeb = .false. 2240 character*8 :: label 2241 2242 call qenter('PELIB_IFC_TRANSFORMER') 2243 2244 if (ipqmmm .gt. 10) then 2245 write(lupri,*) 'PECC TRANSFORMER : ISYMTR : ', isymtr, ' and LR : ', lr 2246 end if 2247 2248 if (ccs) call quit('PECC TRANSFORMER not implemented for CCS') 2249 2250 if (discex) call quit('PECC TRANSFORMER not implemented for DISCEX') 2251 if (hffld) call quit('HFFLD should not be here') 2252 if (ccfixf) then 2253 call qexit('PELIB_IFC_TRANSFORMER') 2254 return 2255 end if 2256 2257 if (lwork .lt. 0) then 2258 write(lupri,*) 'Available LWORK: ', lwork 2259 call quit('Too little work in PECC TRANSFORMER') 2260 end if 2261 2262 allocate(denmats(nnbasx),gmat(n2bst(isymtr)),eta(nt1am(isymtr)+nt2am(isymtr)), & 2263 denmattemp(n2bst(isymtr)),gmattemp(nnbasx)) 2264 eta = zero 2265 denmattemp = zero 2266 gmat = zero 2267 if ((ipqmmm .gt. 10) .or. (locdeb)) then 2268 rho1n = ddot(nt1am(isymtr),rho1,1,rho1,1) 2269 rho2n = ddot(nt2am(isymtr),rho2,1,rho2,1) 2270 write(lupri,*) 'Norm of RHO1 in PECC TRANSFORMER on input: ', rho1n 2271 write(lupri,*) 'Norm of RHO2 in PECC TRANSFORMER on input: ', rho2n 2272 rho1n = ddot(nt1am(isymtr),ctr1,1,ctr1,1) 2273 rho2n = ddot(nt2am(isymtr),ctr2,1,ctr2,1) 2274 write(lupri,*) 'Norm of C1AM in PECC TRANSFORMER on input: ', rho1n 2275 write(lupri,*) 'Norm of C2AM in PECC TRANSFORMER on input: ', rho2n 2276 end if 2277 call ccmm_d1ao(denmattemp,ctr1,ctr2,trim(model),lr,lisdum,idldum, & 2278 isydum,work,lwork) 2279 call dgefsp(nbas,denmattemp,denmats) 2280 call pelib_ifc_response(denmats,gmattemp,1) 2281 call dsptsi(nbas,gmattemp,gmat) 2282 if ((lr .eq. 'L') .or. (lr .eq. 'F')) then 2283 label = 'GIVE INT' 2284 list = 'L0' 2285 idlino = 1 2286 call cc_etac(isymtr,label,eta,list,idlino,0,gmat,work,lwork) 2287 else if ((lr .eq. 'R') .or. (lr .eq. 'P')) then 2288 label = 'GIVE INT' 2289 call cc_xksi(eta,label,isymtr,0,gmat,work,lwork) 2290 if (lr .eq. 'R') then 2291 call cclr_diascl(eta(nt1am(isymtr)+1),two,isymtr) 2292 end if 2293 end if 2294 2295 if ((locdeb) .or. (ipqmmm .gt. 14)) then 2296 tal1 = ddot(nt1am(isymtr),eta,1,eta,1) 2297 tal2 = ddot(nt2am(isymtr),eta(nt1am(isymtr)+1),1,eta(nt1am(isymtr)+1),1) 2298 write(lupri,*) 'Printing TRANSFORMATION contribution. & 2299 Norm2 of singles: ', tal1, 'Norm2 of doubles: ', tal2 2300 end if 2301 2302 fact=one 2303 2304 call daxpy(nt1am(isymtr),fact,eta,1,rho1,1) 2305 call daxpy(nt2am(isymtr),fact,eta(nt1am(isymtr)+1),1,rho2,1) 2306 2307 if ((locdeb) .or. (ipqmmm .gt. 14)) then 2308 tal1 = ddot(nt1am(isymtr),rho1,1,rho1,1) 2309 tal2 = ddot(nt2am(isymtr),rho2,1,rho2,1) 2310 write(lupri,*) 'Printing RHO: & 2311 Norm2 of singles: ', tal1, 'Norm2 of doubles: ', tal2 2312 end if 2313 deallocate(denmats,gmat,eta,denmattemp,gmattemp) 2314 2315 call qexit('PELIB_IFC_TRANSFORMER') 2316 2317end subroutine pelib_ifc_transformer 2318 2319subroutine pelib_ifc_qrtransformer(rho1,rho2,isyres,listb,idlstb,isymtb, & 2320 listc,idlstc,isymtc,model,rsptyp, & 2321 work,lwork) 2322 implicit none 2323 2324#include "mxcent.h" 2325#include "qmmm.h" 2326#include "qm3.h" 2327#include "ccorb.h" 2328#include "ccsdinp.h" 2329#include "ccslvinf.h" 2330#include "ccsdsym.h" 2331#include "priunit.h" 2332 2333 real*8, parameter :: zero = 0.0d0, half = 0.5d0, one = 1.0d0, two = 2.0d0 2334 integer, parameter :: izero = 0 2335 integer :: lwork, ndim, isycb 2336 real*8 :: work(lwork), rho1(*), rho2(*), ddot,rho1n, rho2n, tal1, tal2, & 2337 factks, fact 2338 real*8, allocatable :: b1am(:), b2am(:), c1am(:), c2am(:), denmattemp(:), & 2339 denmats(:), eta1(:), eta2(:), eta3(:), gbmat(:), & 2340 gbmattemp(:) 2341 character*(*) :: listb, listc 2342 integer :: idlstb, isymtb, idlstc, isymtc, isyres, iopt 2343 logical, parameter :: locdeb = .false. 2344 logical :: lsame 2345 character*5 :: model, moddum 2346 character*8 :: label 2347 character*1 :: dentyp, rsptyp 2348 character*2 :: listl 2349 2350 call qenter('PELIB_IFC_QRTRANSFORMER') 2351 2352 listl = 'L0' 2353 2354 if (ccs) call quit('PECC QR TRANSFORMER: Not implemented for CCS') 2355 if (discex) call quit('PECC QR TRANSFORMER: Not implemented for DISCEX') 2356 if ((rsptyp .ne. 'F') .and. (rsptyp .ne. 'G') .and. (rsptyp .ne. 'B') .and. & 2357 (rsptyp .ne. 'K')) then 2358 write(lupri,*) 'Response flag in QRTRANSFORMER is : ', rsptyp 2359 write(lupri,*) 'Response flag in QRTRANSFORMER must be F, G, B or K!' 2360 call quit('Wrong flag in PECC QR TRANSFORMER') 2361 end if 2362 2363! if B=C, just multiply by 2 for second contribution 2364 if (rsptyp .eq. 'F') then 2365 factks = one 2366 lsame = .false. 2367 else 2368 lsame = ((listc .eq. listb) .and. (idlstc .eq. idlstb)) 2369 if (lsame) then 2370 factks = two 2371 else 2372 factks = one 2373 end if 2374 end if 2375 2376 if ((locdeb) .or. (ipqmmm .gt. 10)) then 2377 write(lupri,*) 'PECC QR TRANSFORMER: RSPTYP = ', rsptyp 2378 write(lupri,*) 'PECC QR TRANSFORMER: ISYRES = ', isyres 2379 write(lupri,*) 'PECC QR TRANSFORMER: ISYMTB = ', isymtb 2380 write(lupri,*) 'PECC QR TRANSFORMER: ISYMTC = ', isymtc 2381 write(lupri,*) 'PECC QR TRANSFORMER: LISTB = ', listb 2382 write(lupri,*) 'PECC QR TRANSFORMER: LISTC = ', listc 2383 write(lupri,*) 'PECC QR TRANSFORMER: IDLISTB = ', idlstb 2384 write(lupri,*) 'PECC QR TRANSFORMER: IDLISTC = ', idlstc 2385 write(lupri,*) 'PECC QR TRANSFORMER: LSAME = ', lsame 2386 call flshfo(lupri) 2387 end if 2388 2389 isycb = muld2h(isymtb,isymtc) 2390 if (isycb .ne. isyres) then 2391 call quit('Symmetry problem in PECC QR TRANSFORMER') 2392 end if 2393 if (lwork .lt. 0) then 2394 write(lupri,*) 'Available memory: ', lwork 2395 call quit('Too little work in PECC QR TRANSFORMER (1).') 2396 end if 2397 if ((.not. lsame) .and. (rsptyp .ne. 'K')) then 2398 allocate(b1am(nt1am(isymtb)),b2am(nt2am(isymtb)), & 2399 c1am(nt1am(isymtc)),c2am(nt2am(isymtc)), & 2400 denmattemp(n2bst(isymtb)+n2bst(isymtc)+n2bst(isycb)), & 2401 gbmat(n2bst(isymtb)+n2bst(isymtc)+n2bst(isycb)), & 2402 denmats(3*nnbasx),gbmattemp(3*nnbasx), & 2403 eta1(nt1am(isycb)+nt2am(isycb)), & 2404 eta2(nt1am(isycb)+nt2am(isycb)), & 2405 eta3(nt1am(isycb)+nt2am(isycb))) 2406 ndim = 3 2407 denmattemp = zero 2408 gbmattemp = zero 2409 eta1 = zero 2410 eta2 = zero 2411 eta3 = zero 2412 fact = one 2413 else if ((lsame) .and. (rsptyp .ne. 'K')) then 2414 allocate(b1am(nt1am(isymtb)),b2am(nt2am(isymtb)), & 2415 denmattemp(n2bst(isymtb)+n2bst(isymtc)), & 2416 denmats(2*nnbasx),gbmattemp(2*nnbasx),& 2417 gbmat(n2bst(isymtb)+n2bst(isymtc)), & 2418 eta1(nt1am(isycb)+nt2am(isycb)), & 2419 eta2(nt1am(isycb)+nt2am(isycb))) 2420 ndim = 2 2421 denmattemp = zero 2422 gbmattemp = zero 2423 eta1 = zero 2424 eta2 = zero 2425 fact = one 2426 else if ((.not.lsame) .and. ((rsptyp .eq. 'K'))) then 2427 allocate(b1am(nt1am(isymtb)),b2am(nt2am(isymtb)), & 2428 c1am(nt1am(isymtc)),c2am(nt2am(isymtc)), & 2429 denmattemp(n2bst(isymtb)+n2bst(isymtc)), & 2430 denmats(2*nnbasx),gbmattemp(2*nnbasx), & 2431 gbmat(n2bst(isymtb)+n2bst(isymtc)), & 2432 eta1(nt1am(isycb)+nt2am(isycb)), & 2433 eta2(nt1am(isycb)+nt2am(isycb))) 2434 ndim = 2 2435 denmattemp = zero 2436 gbmattemp = zero 2437 eta1 = zero 2438 eta2 = zero 2439 fact = one 2440 else if ((lsame) .and. ((rsptyp .eq. 'K'))) then 2441 allocate(b1am(nt1am(isymtb)),b2am(nt2am(isymtb)), & 2442 denmattemp(n2bst(isymtb)),denmats(nnbasx), & 2443 gbmat(n2bst(isycb)),gbmattemp(nnbasx), & 2444 eta1(nt1am(isycb)+nt2am(isycb))) 2445 ndim = 1 2446 denmattemp = zero 2447 gbmattemp = zero 2448 eta1 = zero 2449 fact = one 2450 end if 2451 2452 iopt = 3 2453 call cc_rdrsp(listb,idlstb,isymtb,iopt,model,b1am,b2am) 2454 if (.not. lsame) then 2455 call cc_rdrsp(listc,idlstc,isymtc,iopt,model,c1am,c2am) 2456 end if 2457 2458 if ((locdeb) .or. (ipqmmm .gt. 10)) then 2459 rho1n = ddot(nt1am(isycb),rho1,1,rho1,1) 2460 rho2n = ddot(nt2am(isycb),rho2,1,rho2,1) 2461 write(lupri,*) 'Norm of RHO1 in PECC QM TRANSFORMER on input (1): ', & 2462 rho1n 2463 write(lupri,*) 'Norm of RHO2 in PECC QM TRANSFORMER on input (1): ', & 2464 rho2n 2465 rho1n = ddot(nt1am(isycb),b1am,1,b1am,1) 2466 rho2n = ddot(nt2am(isycb),b2am,1,b2am,1) 2467 write(lupri,*) 'Norm of B1AM in PECC QM TRANSFORMER on input (1): ', & 2468 rho1n 2469 write(lupri,*) 'Norm of B2AM in PECC QM TRANSFORMER on input (1): ', & 2470 rho2n 2471 if (.not. lsame) then 2472 rho1n = ddot(nt1am(isycb),c1am,1,c1am,1) 2473 rho2n = ddot(nt2am(isycb),c2am,1,c2am,1) 2474 write(lupri,*) 'Norm of C1AM in PECC QM TRANSFORMER on input (1): ', & 2475 rho1n 2476 write(lupri,*) 'Norm of C2AM in PECC QM TRANSFORMER on input (1): ', & 2477 rho2n 2478 end if 2479 end if 2480 2481! check kind of response (for F start with left density) 2482 if ((rsptyp .eq. 'G') .or. (rsptyp .eq. 'B')) then 2483 dentyp = 'R' 2484 else if ((rsptyp .eq. 'K') .or. (rsptyp .eq. 'F')) then 2485 dentyp = 'L' 2486 end if 2487 call ccmm_d1ao(denmattemp,b1am,b2am,model,dentyp,listc,idlstc,isymtc,& 2488 work,lwork) 2489 if (.not. lsame) then 2490 if (rsptyp .eq. 'F') dentyp = 'R' 2491 call ccmm_d1ao(denmattemp(n2bst(isymtb)+1),c1am,c2am,model,dentyp, & 2492 listb,idlstb,isymtb,work,lwork) 2493 if ((rsptyp .eq. 'G').or.(rsptyp .eq. 'B')) then 2494 call ccmm_d2ao(denmattemp(n2bst(isymtb)+n2bst(isymtc)+1),isyres, & 2495 listb,idlstb,isymtb,listc,idlstc,isymtc,model, & 2496 work,lwork) 2497 else if (rsptyp .eq. 'F') then 2498 dentyp = 'Q' 2499 call ccmm_d1ao(denmattemp(n2bst(isymtb)+n2bst(isymtc)+1),c1am,c2am, & 2500 model,dentyp,listb,idlstb,isymtb,work,lwork) 2501 end if 2502 else if (lsame) then 2503 if ((rsptyp .eq. 'G') .or. (rsptyp .eq. 'B')) then 2504 call ccmm_d2ao(denmattemp(n2bst(isymtb)+1),isyres,listb,idlstb, & 2505 isymtb,listc,idlstc,isymtc,model,work,lwork) 2506 end if 2507 end if 2508! construct effective G operator 2509 call dgefsp(nbas,denmattemp,denmats) 2510 if (.not. lsame) then 2511 call dgefsp(nbas,denmattemp(n2bst(isymtb)+1),denmats(nnbasx+1)) 2512 if (rsptyp .eq. 'B') then 2513 call dgefsp(nbas,denmattemp(n2bst(isymtb)+n2bst(isymtc)+1), & 2514 denmats(2*nnbasx+1)) 2515 else if ((rsptyp .eq. 'F').or.(rsptyp.eq.'G')) then 2516 call dgefsp(nbas,denmattemp(n2bst(isymtb)+n2bst(isymtc)+1), & 2517 denmats(2*nnbasx+1)) 2518 end if 2519 else if (lsame) then 2520 if ((rsptyp .eq. 'G') .or. (rsptyp .eq. 'B')) then 2521 call dgefsp(nbas,denmattemp(n2bst(isymtb)+1),denmats(nnbasx+1)) 2522 end if 2523 end if 2524 call pelib_ifc_response(denmats,gbmattemp,ndim) 2525 call dsptsi(nbas,gbmattemp,gbmat) 2526 if (.not. lsame) then 2527 call dsptsi(nbas,gbmattemp(nnbasx+1),gbmat(n2bst(isymtb)+1)) 2528 if ((rsptyp .eq. 'B')) then 2529 call dsptsi(nbas,gbmattemp(2*nnbasx+1),gbmat(n2bst(isymtb)+ & 2530 n2bst(isymtc)+1)) 2531 else if ((rsptyp .eq. 'F').or.(rsptyp.eq.'G')) then 2532 call dsptsi(nbas,gbmattemp(2*nnbasx+1),gbmat(n2bst(isymtb)+ & 2533 n2bst(isymtc)+1)) 2534 end if 2535 else if (lsame) then 2536 if ((rsptyp .eq. 'G') .or. (rsptyp .eq. 'B')) then 2537 call dsptsi(nbas,gbmattemp(nnbasx+1),gbmat(n2bst(isymtb)+1)) 2538 end if 2539 end if 2540 if ((locdeb) .or. (ipqmmm .gt. 14)) then 2541 tal1 = ddot(n2bst(isymtb),gbmat,1,gbmat,1) 2542 write(lupri,*) 'Print Norm2 GBMAT (1): ', tal1 2543 if (.not. lsame) then 2544 tal1 = ddot(n2bst(isymtc),gbmat(n2bst(isymtb)+1),1, & 2545 gbmat(n2bst(isymtb)+1),1) 2546 write(lupri,*) 'Print Norm2 GBMAT (2): ', tal1 2547 if ((rsptyp.eq.'B').or.(rsptyp.eq.'G').or.(rsptyp.eq.'F')) then 2548 tal1 = ddot(n2bst(isymtb),gbmat(n2bst(isymtb)+n2bst(isymtc)+1), & 2549 1,gbmat(n2bst(isymtb)+n2bst(isymtc)+1),1) 2550 write(lupri,*) 'Print Norm2 GBCMAT (3): ', tal1 2551 end if 2552 else if (lsame) then 2553 if ((rsptyp .eq. 'G') .or. (rsptyp .eq. 'B')) then 2554 tal1 = ddot(n2bst(isymtc),gbmat(n2bst(isymtb)+1),1, & 2555 gbmat(n2bst(isymtb)+1),1) 2556 write(lupri,*) 'Print Norm2 GBMAT (2): ', tal1 2557 end if 2558 end if 2559 end if 2560 label = 'GIVE INT' 2561 if ((rsptyp .eq. 'G') .or. (rsptyp .eq. 'F')) then 2562 call cclr_fa(label,isymtb,listc,idlstc,listl,0,gbmat, & 2563 work,lwork) 2564 call daxpy(nt1am(isycb)+nt2am(isycb),1.0d0,work,1,eta1,1) 2565 else if (rsptyp .eq. 'B') then 2566 call cccr_aa(label,isymtb,listc,idlstc,gbmat,work,lwork) 2567 call daxpy(nt1am(isycb)+nt2am(isycb),1.0d0,work,1,eta1,1) 2568 call cclr_diascl(eta1(nt1am(isycb)+1),two,isycb) 2569 else if (rsptyp .eq. 'K') then 2570 call cc_etac(isymtb,label,eta1,listc,idlstc,0,gbmat,work,lwork) 2571 end if 2572 if (locdeb .or.(ipqmmm .gt. 14)) then 2573 tal1 = ddot(nt1am(isycb),eta1,1,eta1,1) 2574 tal2 = ddot(nt2am(isycb),eta1(nt1am(isycb)+1),1,eta1(nt1am(isycb)+1),1) 2575 write(lupri,*) 'Printing transformed G^B*C contribution.' 2576 write(lupri,*) 'Norm2 of singles : ', tal1 2577 write(lupri,*) 'Norm2 of doubles : ', tal2 2578 end if 2579 if (.not. lsame) then 2580 if (rsptyp .eq. 'G') then 2581 call cclr_fa(label,isymtc,listb,idlstb,listl,0,gbmat(n2bst(isymtb)+1), & 2582 work,lwork) 2583 call daxpy(nt1am(isycb)+nt2am(isycb),1.0d0,work,1,eta2,1) 2584 else if (rsptyp .eq. 'B') then 2585 call cccr_aa(label,isymtc,listb,idlstb,gbmat(n2bst(isymtb)+1),work,lwork) 2586 call daxpy(nt1am(isycb)+nt2am(isycb),1.0d0,work,1,eta2,1) 2587 call cclr_diascl(eta2(nt1am(isycb)+1),two,isycb) 2588 else if ((rsptyp .eq. 'K') .or. (rsptyp .eq. 'F')) then 2589 call cc_etac(isymtc,label,eta2,listb,idlstb,0,gbmat(n2bst(isymtb)+1), & 2590 work,lwork) 2591 end if 2592 if (locdeb .or.(ipqmmm .gt. 14)) then 2593 tal1 = ddot(nt1am(isycb),eta2,1,eta2,1) 2594 tal2 = ddot(nt2am(isycb),eta2(nt1am(isycb)+1),1,eta2(nt1am(isycb)+1),1) 2595 write(lupri,*) 'Printing transformed G^C*B contribution.' 2596 write(lupri,*) 'Norm2 of singles : ', tal1 2597 write(lupri,*) 'Norm2 of doubles : ', tal2 2598 end if 2599 if ((rsptyp .eq. 'G').or.(rsptyp .eq. 'F')) then 2600 call cc_etac(isycb,label,eta3,listl,0,0,gbmat(n2bst(isymtb)+ & 2601 n2bst(isymtc)+1),work,lwork) 2602 else if ((rsptyp .eq. 'B')) then 2603 call cc_xksi(eta3,label,isycb,0,gbmat(n2bst(isymtb)+n2bst(isymtc)+1), & 2604 work,lwork) 2605 call cclr_diascl(eta3(nt1am(isycb)+1),two,isycb) 2606 end if 2607 if (locdeb .or.(ipqmmm .gt. 14)) then 2608 if ((rsptyp .eq. 'G').or.(rsptyp .eq. 'B').or.(rsptyp.eq.'F')) then 2609 tal1 = ddot(nt1am(isycb),eta3,1,eta3,1) 2610 tal2 = ddot(nt2am(isycb),eta3(nt1am(isycb)+1),1,eta3(nt1am(isycb)+1),1) 2611 write(lupri,*) 'Printing transformed G^{BC} contribution.' 2612 write(lupri,*) 'Norm2 of singles : ', tal1 2613 write(lupri,*) 'Norm2 of doubles : ', tal2 2614 end if 2615 end if 2616 else if (lsame) then 2617 if (rsptyp .eq. 'G') then 2618 call cc_etac(isycb,label,eta2,listl,0,0,gbmat(n2bst(isymtb)+1), & 2619 work,lwork) 2620 else if ((rsptyp .eq. 'B')) then 2621 call cc_xksi(eta2, label, isycb, 0, gbmat(n2bst(isymtb)+1), work,lwork) 2622 call cclr_diascl(eta2(nt1am(isycb)+1),two,isycb) 2623 end if 2624 if (locdeb .or.(ipqmmm .gt. 14)) then 2625 if ((rsptyp .eq. 'G').or.(rsptyp .eq. 'B')) then 2626 tal1 = ddot(nt1am(isycb),eta2,1,eta2,1) 2627 tal2 = ddot(nt2am(isycb),eta2(nt1am(isycb)+1),1,eta2(nt1am(isycb)+1),1) 2628 write(lupri,*) 'Printing transformed G^{BC} contribution.' 2629 write(lupri,*) 'Norm2 of singles : ', tal1 2630 write(lupri,*) 'Norm2 of doubles : ', tal2 2631 end if 2632 end if 2633 end if 2634 2635 if ((locdeb) .or. (ipqmmm .gt. 14)) then 2636 tal1 = ddot(nt1am(isycb),eta1,1,eta1,1) 2637 tal2 = ddot(nt2am(isycb),eta1(nt1am(isycb)+1),1,eta1(nt1am(isycb)+1),1) 2638 write(lupri,*) 'Printing transformed G^B*C contribution.' 2639 write(lupri,*) 'Norm2 of singles: ', tal1 2640 write(lupri,*) 'Norm2 of doubles: ', tal2 2641 end if 2642 call daxpy(nt1am(isycb),factks*fact,eta1,1,rho1,1) 2643 call daxpy(nt2am(isycb),factks*fact,eta1(nt1am(isycb)+1),1,rho2,1) 2644 if (.not. lsame) then 2645 call daxpy(nt1am(isycb),one*fact,eta2,1,rho1,1) 2646 call daxpy(nt2am(isycb),one*fact,eta2(nt1am(isycb)+1),1,rho2,1) 2647 if ((rsptyp.eq.'G').or.(rsptyp.eq.'B').or.(rsptyp.eq.'F')) then 2648 call daxpy(nt1am(isycb),one*fact,eta3,1,rho1,1) 2649 call daxpy(nt2am(isycb),one*fact,eta3(nt1am(isycb)+1),1,rho2,1) 2650 end if 2651 else if (lsame) then 2652 if ((rsptyp .eq. 'G') .or. (rsptyp .eq. 'B')) then 2653 call daxpy(nt1am(isycb),one*fact,eta2,1,rho1,1) 2654 call daxpy(nt2am(isycb),one*fact,eta2(nt1am(isycb)+1),1,rho2,1) 2655 end if 2656 end if 2657 2658 if ((locdeb) .or. (ipqmmm .gt. 14)) then 2659 tal1 = ddot(nt1am(isycb),rho1,1,rho1,1) 2660 tal2 = ddot(nt2am(isycb),rho2,1,rho2,1) 2661 write(lupri,*) 'Printing RHO.' 2662 write(lupri,*) 'Norm2 of singles: ', tal1 2663 write(lupri,*) 'Norm2 of doubles: ', tal2 2664 end if 2665 2666 deallocate(denmattemp,denmats,gbmat,gbmattemp,b1am,b2am,eta1) 2667 if (allocated(c1am)) deallocate(c1am) 2668 if (allocated(c2am)) deallocate(c2am) 2669 if (allocated(eta2)) deallocate(eta2) 2670 if (allocated(eta3)) deallocate(eta3) 2671 2672 call qexit('PELIB_IFC_QRTRANSFORMER') 2673 2674end subroutine pelib_ifc_qrtransformer 2675 2676end module pelib_interface 2677 2678subroutine pelib_ifc_start_slaves(runtyp) 2679 use pelib_interface, only: use_pelib 2680 integer :: runtyp 2681#include "iprtyp.h" 2682#include "maxorb.h" 2683#include "infpar.h" 2684 integer, parameter :: iprtyp = POLARIZABLE_EMBEDDING 2685 call qenter('pelib_ifc_start_slaves') 2686 if (nodtot >= 1) then 2687 call mpixbcast(iprtyp, 1, 'INTEGER', master) 2688 call mpixbcast(runtyp, 1, 'INTEGER', master) 2689 end if 2690 call qexit('pelib_ifc_start_slaves') 2691end subroutine pelib_ifc_start_slaves 2692 2693#else 2694 2695module pelib_interface 2696 2697 implicit none 2698 2699 private 2700 2701 public :: use_pelib, pelib_ifc_gspol 2702 public :: pelib_ifc_do_mep, pelib_ifc_do_mep_noqm, pelib_ifc_do_cube 2703 public :: pelib_ifc_do_lf 2704 public :: pelib_ifc_activate, pelib_ifc_deactivate 2705 public :: pelib_ifc_init, pelib_ifc_finalize, pelib_ifc_input_reader 2706 public :: pelib_ifc_fock, pelib_ifc_energy, pelib_ifc_response, pelib_ifc_london 2707 public :: pelib_ifc_molgrad, pelib_ifc_lf, pelib_ifc_localfield 2708 public :: pelib_ifc_mep, pelib_ifc_mep_noqm, pelib_ifc_cube 2709 public :: pelib_ifc_set_mixed, pelib_ifc_mixed 2710 public :: pelib_ifc_do_savden, pelib_ifc_do_twoints 2711 public :: pelib_ifc_save_density, pelib_ifc_twoints 2712 public :: pelib_ifc_get_num_core_nuclei 2713#if defined(VAR_MPI) 2714 public :: pelib_ifc_slave 2715#endif 2716 public :: pelib_ifc_grad, pelib_ifc_lin, pelib_ifc_lr, pelib_ifc_qro, pelib_ifc_cro 2717 public :: pelib_ifc_pecc, pelib_ifc_transformer, pelib_ifc_qrtransformer 2718 2719contains 2720 2721logical function use_pelib() 2722 use_pelib = .false. 2723end function use_pelib 2724 2725logical function pelib_ifc_gspol() 2726 pelib_ifc_gspol = .false. 2727end function pelib_ifc_gspol 2728 2729subroutine pelib_ifc_set_mixed(do_mixed) 2730 logical :: do_mixed 2731 call quit('using dummy PElib interface routines') 2732end subroutine pelib_ifc_set_mixed 2733 2734logical function pelib_ifc_do_mep() 2735 pelib_ifc_do_mep = .false. 2736end function pelib_ifc_do_mep 2737 2738logical function pelib_ifc_do_mep_noqm() 2739 pelib_ifc_do_mep_noqm = .false. 2740end function pelib_ifc_do_mep_noqm 2741 2742logical function pelib_ifc_do_cube() 2743 pelib_ifc_do_cube = .false. 2744end function pelib_ifc_do_cube 2745 2746logical function pelib_ifc_do_lf() 2747 pelib_ifc_do_lf = .false. 2748end function pelib_ifc_do_lf 2749 2750logical function pelib_ifc_do_savden() 2751 pelib_ifc_do_savden = .false. 2752end function pelib_ifc_do_savden 2753 2754logical function pelib_ifc_do_twoints() 2755 pelib_ifc_do_twoints = .false. 2756end function pelib_ifc_do_twoints 2757 2758subroutine pelib_ifc_activate() 2759 call qenter('pelib_ifc_activate') 2760 call quit('PElib not compiled, please use -DENABLE_PELIB=ON to enable PElib') 2761 call qexit('pelib_ifc_activate') 2762end subroutine pelib_ifc_activate 2763 2764subroutine pelib_ifc_deactivate() 2765 call qenter('pelib_ifc_deactivate') 2766 call quit('using dummy PElib interface routines') 2767 call qexit('pelib_ifc_deactivate') 2768end subroutine pelib_ifc_deactivate 2769 2770subroutine pelib_ifc_input_reader(word) 2771 character(len=7), intent(in) :: word 2772 call qenter('pelib_ifc_input_reader') 2773 call quit('using dummy PElib interface routines') 2774 call qexit('pelib_ifc_input_reader') 2775end subroutine pelib_ifc_input_reader 2776 2777subroutine pelib_ifc_init() 2778 call qenter('pelib_ifc_init') 2779 call quit('using dummy PElib interface routines') 2780 call qexit('pelib_ifc_init') 2781end subroutine pelib_ifc_init 2782 2783subroutine pelib_ifc_finalize() 2784 call qenter('pelib_ifc_finalize') 2785 call quit('using dummy PElib interface routines') 2786 call qexit('pelib_ifc_finalize') 2787end subroutine pelib_ifc_finalize 2788 2789subroutine pelib_ifc_fock(denmats, fckmats, energy) 2790 real*8, dimension(*), intent(in) :: denmats 2791 real*8, dimension(*), intent(out) :: fckmats 2792 real*8, intent(out) :: energy 2793 call qenter('pelib_ifc_fock') 2794 call quit('using dummy PElib interface routines') 2795 call qexit('pelib_ifc_fock') 2796end subroutine pelib_ifc_fock 2797 2798subroutine pelib_ifc_mixed(denmats, fckmats, energy) 2799 real*8, dimension(*), intent(in) :: denmats 2800 real*8, dimension(*), intent(in) :: fckmats 2801 real*8, intent(in) :: energy 2802 call qenter('pelib_ifc_mixed') 2803 call quit('using dummy PElib interface routines') 2804 call qexit('pelib_ifc_mixed') 2805end subroutine pelib_ifc_mixed 2806 2807subroutine pelib_ifc_energy(denmats) 2808 real*8, dimension(*), intent(in) :: denmats 2809 call qenter('pelib_ifc_energy') 2810 call quit('using dummy PElib interface routines') 2811 call qexit('pelib_ifc_energy') 2812end subroutine pelib_ifc_energy 2813 2814subroutine pelib_ifc_molgrad(denmats, molgrad) 2815 real*8, dimension(*), intent(in) :: denmats 2816 real*8, dimension(*), intent(out) :: molgrad 2817 call qenter('pelib_ifc_molgrad') 2818 call quit('using dummy PElib interface routines') 2819 call qexit('pelib_ifc_molgrad') 2820end subroutine pelib_ifc_molgrad 2821 2822subroutine pelib_ifc_response(denmats, fckmats, nmats) 2823 integer, intent(in) :: nmats 2824 real*8, dimension(*), intent(in) :: denmats 2825 real*8, dimension(*), intent(out) :: fckmats 2826 call qenter('pelib_ifc_response') 2827 call quit('using dummy PElib interface routines') 2828 call qexit('pelib_ifc_response') 2829end subroutine pelib_ifc_response 2830 2831subroutine pelib_ifc_london(fckmats) 2832 real*8, dimension(*), intent(out) :: fckmats 2833 call qenter('pelib_ifc_london') 2834 call quit('using dummy PElib interface routines') 2835 call qexit('pelib_ifc_london') 2836end subroutine pelib_ifc_london 2837 2838subroutine pelib_ifc_localfield(eefmats) 2839 real*8, dimension(:), intent(in) :: eefmats 2840 call qenter('pelib_ifc_localfield') 2841 call quit('using dummy PElib interface routines') 2842 call qexit('pelib_ifc_localfield') 2843end subroutine pelib_ifc_localfield 2844 2845subroutine pelib_ifc_lf() 2846 call qenter('pelib_ifc_lf') 2847 call quit('using dummy PElib interface routines') 2848 call qexit('pelib_ifc_lf') 2849end subroutine pelib_ifc_lf 2850 2851subroutine pelib_ifc_mep(denmats) 2852 real*8, dimension(*), intent(in) :: denmats 2853 call qenter('pelib_ifc_mep') 2854 call quit('using dummy PElib interface routines') 2855 call qexit('pelib_ifc_mep') 2856end subroutine pelib_ifc_mep 2857 2858subroutine pelib_ifc_mep_noqm() 2859 call qenter('pelib_ifc_mep_noqm') 2860 call quit('using dummy PElib interface routines') 2861 call qexit('pelib_ifc_mep_noqm') 2862end subroutine pelib_ifc_mep_noqm 2863 2864subroutine pelib_ifc_cube(denmats, idx) 2865 real*8, dimension(*), intent(in) :: denmats 2866 integer, intent(in) :: idx 2867 call qenter('pelib_ifc_cube') 2868 call quit('using dummy PElib interface routines') 2869 call qexit('pelib_ifc_cube') 2870end subroutine pelib_ifc_cube 2871 2872subroutine pelib_ifc_save_density(ao_denmat, mo_fckmat, mo_coefficients) 2873 real*8, dimension(*), intent(in) :: ao_denmat 2874 real*8, dimension(*), intent(in) :: mo_fckmat 2875 real*8, dimension(*), intent(in) :: mo_coefficients 2876 call qenter('pelib_ifc_save_density') 2877 call quit('using dummy PElib interface routines') 2878 call qexit('pelib_ifc_save_density') 2879end subroutine pelib_ifc_save_density 2880 2881subroutine pelib_ifc_twoints(work, lwork) 2882 real*8, dimension(*), intent(in) :: work 2883 integer, intent(in) :: lwork 2884 call qenter('pelib_ifc_twoints') 2885 call quit('using dummy PElib interface routines') 2886 call qexit('pelib_ifc_twoints') 2887end subroutine pelib_ifc_twoints 2888 2889integer function pelib_ifc_get_num_core_nuclei() 2890 call quit('using dummy PElib interface routines') 2891end function pelib_ifc_get_num_core_nuclei 2892 2893#if defined(VAR_MPI) 2894subroutine pelib_ifc_slave(runtype) 2895 integer, intent(in) :: runtype 2896 call qenter('pelib_ifc_slave') 2897 call quit('using dummy PElib interface routines') 2898 call qexit('pelib_ifc_slave') 2899end subroutine pelib_ifc_slave 2900#endif 2901 2902subroutine pelib_ifc_grad(cref, cmo, cindx, dv, grd, energy, wrk, nwrk) 2903 integer :: nwrk 2904 real*8 :: energy 2905 real*8, dimension(*) :: cref, cmo, cindx, dv, grd 2906 real*8, dimension(*) :: wrk 2907 call qenter('pelib_ifc_grad') 2908 call quit('using dummy PElib interface routines') 2909 call qexit('pelib_ifc_grad') 2910end subroutine pelib_ifc_grad 2911 2912subroutine pelib_ifc_lin(ncsim, nosim, bcvecs, bovecs, cref, cmo, cindx, dv, dtv,& 2913 scvecs, sovecs, orblin, wrk, nwrk) 2914 logical :: orblin 2915 integer :: ncsim, nosim, nwrk 2916 integer, dimension(*) :: cindx 2917 real*8, dimension(*) :: bcvecs, bovecs, scvecs, sovecs 2918 real*8, dimension(*) :: cmo, cref, dv, dtv 2919 real*8, dimension(*) :: wrk 2920 call qenter('pelib_ifc_lin') 2921 call quit('using dummy PElib interface routines') 2922 call qexit('pelib_ifc_lin') 2923end subroutine pelib_ifc_lin 2924 2925subroutine pelib_ifc_lr(ncsim, nosim, bcvecs, bovecs, cref, cmo, cindx, udv,& 2926 dv, udvtr, dvtr, dtv, dtvtr, scvecs, sovecs, wrk, nwrk) 2927 integer :: ncsim, nosim, nwrk 2928 integer, dimension(*) :: cindx 2929 real*8, dimension(*) :: bcvecs, bovecs 2930 real*8, dimension(*) :: cref, cmo, udv, dv 2931 real*8, dimension(*) :: udvtr, dvtr, dtv, dtvtr 2932 real*8, dimension(*) :: scvecs, sovecs 2933 real*8, dimension(nwrk) :: wrk 2934 call qenter('pelib_ifc_lr') 2935 call quit('using dummy PElib interface routines') 2936 call qexit('pelib_ifc_lr') 2937end subroutine pelib_ifc_lr 2938 2939subroutine pelib_ifc_qro(vecb, vecc, etrs, xindx, zymb, zymc, udv, wrk, nwrk,& 2940 kzyva, kzyvb, kzyvc, isyma, isymb, isymc, cmo, mjwop) 2941 integer :: kzyva, kzyvb, kzyvc 2942 integer :: isyma, isymb, isymc 2943 integer :: nwrk 2944 real*8, dimension(*) :: wrk 2945 real*8, dimension(*) :: etrs 2946 real*8, dimension(*) :: vecb 2947 real*8, dimension(*) :: vecc 2948 real*8, dimension(*) :: cmo 2949 real*8, dimension(*) :: zymb, zymc 2950 real*8, dimension(*) :: udv 2951 real*8, dimension(*) :: xindx 2952 integer, dimension(*) :: mjwop 2953 call qenter('pelib_ifc_qro') 2954 call quit('using dummy PElib interface routines') 2955 call qexit('pelib_ifc_qro') 2956end subroutine pelib_ifc_qro 2957 2958subroutine pelib_ifc_cro(vecb, vecc, vecd, etrs, xindx, zymb, zymc, zymd, udv,& 2959 wrk, nwrk, kzyva, kzyvb, kzyvc, kzyvd, isyma, isymb,& 2960 isymc, isymd, cmo,mjwop) 2961 integer :: kzyva, kzyvb, kzyvc, kzyvd 2962 integer :: isyma, isymb, isymc, isymd 2963 integer :: nwrk 2964 real*8, dimension(*) :: wrk 2965 real*8, dimension(*) :: etrs 2966 real*8, dimension(*) :: vecb 2967 real*8, dimension(*) :: vecc 2968 real*8, dimension(*) :: vecd 2969 real*8, dimension(*) :: cmo 2970 real*8, dimension(*) :: zymb, zymc, zymd 2971 real*8, dimension(*) :: udv 2972 real*8, dimension(*) :: xindx 2973 integer, dimension(*) :: mjwop 2974 call qenter('pelib_ifc_cro') 2975 call quit('using dummy PElib interface routines') 2976 call qexit('pelib_ifc_cro') 2977end subroutine pelib_ifc_cro 2978 2979subroutine pelib_ifc_pecc(aoden, aodencc, converged, t_or_tbar) 2980 real*8, intent(in) :: aoden(*), aodencc(*) 2981 logical :: converged 2982 integer, intent(in) :: t_or_tbar 2983 call qenter('PELIB_IFC_PECC') 2984 call quit('using dummy PElib interface routines') 2985 call qexit('PELIB_IFC_PECC') 2986end subroutine pelib_ifc_pecc 2987 2988subroutine pelib_ifc_transformer(rho1,rho2,ctr1,ctr2,model,isymtr,lr,work,lwork) 2989 integer :: lwork, isymtr 2990 character*1, intent(in) :: lr 2991 real*8, dimension(*) :: work 2992 real*8 :: rho1(*), rho2(*), ctr1(*), ctr2(*) 2993 character*10 :: model 2994 call qenter('PELIB_IFC_TRANSFORMER') 2995 call quit('using dummy PElib interface routines') 2996 call qexit('PELIB_IFC_TRANSFORMER') 2997end subroutine pelib_ifc_transformer 2998 2999subroutine pelib_ifc_qrtransformer(rho1,rho2,isyres,listb,idlstb,isymtb, & 3000 & listc,idlstc,isymtc,model,rsptyp, & 3001 & work,lwork) 3002 integer :: lwork 3003 real*8 :: work(*), rho1(*), rho2(*) 3004 character*(*) :: listb, listc 3005 integer :: idlstb, isymtb, idlstc, isymtc, isyres 3006 character*5 :: model 3007 character*1 :: rsptyp 3008 call qenter('PELIB_IFC_QRTRANSFORMER') 3009 call quit('using dummy PElib interface routines') 3010 call qexit('PELIB_IFC_QRTRANSFORMER') 3011end subroutine pelib_ifc_qrtransformer 3012 3013end module pelib_interface 3014 3015subroutine pelib_ifc_start_slaves(runtyp) 3016 integer :: runtyp 3017 call qenter('pelib_ifc_start_slaves') 3018 call quit('using dummy PElib interface routines') 3019 call qexit('pelib_ifc_start_slaves') 3020end subroutine pelib_ifc_start_slaves 3021 3022subroutine pelib_ifc_pecc(aoden,converged,work,lwork) 3023 implicit none 3024 integer :: lwork 3025 real*8 :: work(lwork), aoden(*) 3026 logical :: converged 3027 call qenter('PELIB_IFC_PECC') 3028 call quit('using dummy PElib interface routines') 3029 call qexit('PELIB_IFC_PECC') 3030end subroutine pelib_ifc_pecc 3031 3032#endif 3033