1! Dalton, a molecular electronic structure program 2! Copyright (C) by the authors of Dalton. 3! 4! This program is free software; you can redistribute it and/or 5! modify it under the terms of the GNU Lesser General Public 6! License version 2.1 as published by the Free Software Foundation. 7! 8! This program is distributed in the hope that it will be useful, 9! but WITHOUT ANY WARRANTY; without even the implied warranty of 10! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 11! Lesser General Public License for more details. 12! 13! If a copy of the GNU LGPL v2.1 was not distributed with this 14! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 15 16module qmcmm 17 18 implicit none 19 20 public read_pot_qmnpmm 21 public getdim_relmat 22 public getdim_mmmat 23 public comp_relmat 24 public read_relmat 25 public write_relmat 26 public comp_mmrelmat 27 public comp_dampvmat 28 29 private 30 31contains 32 33 34 SUBROUTINE READ_POT_QMNPMM() 35! 36! Purpose: 37! reads NP and MM subsystems geometry, total charges, and 38! force field data from POTENTIAL.INP file 39! 40! Last updated: 22/03/2013 by Z. Rinkevicius. 41! 42 43#include "codata.h" 44#include "priunit.h" 45#include "qmnpmm.h" 46! 47 real(8), parameter :: xfact = 1.0d0/xtang 48 integer :: i, istart, j, joff, luqmnp, idummy 49! 50 CHARACTER UNITS*2, NPWORD*2, FFWORD*4, NPLBL(MXNPATM)*2 51 CHARACTER MMWORD*2, MMLBL(MXMMATM)*2 52! 53! Open POTENTIAL.INP file 54 LUQMNP=-1 55 CALL GPOPEN(LUQMNP,'POTENTIAL.INP','OLD',' ', & 56 'FORMATTED',IDUMMY,.FALSE.) 57 REWIND(LUQMNP) 58! Read geometry input units, number of NP and MM 59! separate subsystems 60 READ(LUQMNP,*) UNITS, TNPBLK, TMMBLK 61 CALL UPCASE(UNITS) 62! Check input consistency 63 IF ((UNITS.NE.'AA').AND.(UNITS.NE.'AU')) THEN 64 CALL QUIT('Unknown units in POTENTIAL.INP') 65 ENDIF 66 IF ((DONPSUB).AND.(TNPBLK.LT.1)) THEN 67 CALL QUIT('NP system is missing in POTENTIAL.INP') 68 END IF 69 IF ((DOMMSUB).AND.(TMMBLK.LT.1)) THEN 70 CALL QUIT('MM system is missing in POTENTIAL.INP') 71 END IF 72 IF (TNPBLK.GT.MAXBLK) THEN 73 WRITE(LUPRI,'(/2X,A)') & 74 'Maximum number of NP subsystems exceeded!' 75 WRITE(LUPRI,'(2X,A,I3,2X,A,I3)') 'Input TNPBLK=', & 76 TNPBLK, 'Maximum allowed:', MAXBLK 77 CALL QUIT('Increase MAXBLK in qmnpmm.h') 78 END IF 79 IF (TMMBLK.GT.MAXBLK) THEN 80 WRITE(LUPRI,'(/2X,A)') & 81 'Maximum number of MM subsystems exceeded!' 82 WRITE(LUPRI,'(2X,A,I3,2X,A,I3)') 'Input TNPBLK=', & 83 TNPBLK, 'Maximum allowed:', MAXBLK 84 CALL QUIT('Increase MAXBLK in qmnpmm.h') 85 END IF 86! Read NP subsystems data 87 IF (DONPSUB) THEN 88 ISTART = 0 89 DO I=1,TNPBLK 90! Read NP subsystem header 91 READ(LUQMNP,*) NPWORD, NPCHRG(I), NPATOM(I) 92! Check input consistency 93 CALL UPCASE(NPWORD) 94 IF (NPWORD.NE.'NP') THEN 95 WRITE(LUPRI,'(/2X,A,I2,1X,A)') & 96 'Incorrect NP subsystem I=', I, & 97 'header in POTENTIAL.INP file!' 98 CALL QUIT('Corrupted POTENTIAL.INP') 99 END IF 100 IF (NPATOM(I).LE.0) THEN 101 WRITE(LUPRI,'(/2X,A,I5,A,I2,A)') & 102 'Incorrect number of atoms ', NPATOM(I), & 103 'in NP subsystem', I, '!' 104 CALL QUIT('Corrupted POTENTIAL.INP') 105 END IF 106! Read NP subsystem data 107 DO J=1,NPATOM(I) 108 JOFF = ISTART + J 109 READ(LUQMNP,*) NPLBL(JOFF), NPCORD(1,JOFF), & 110 NPCORD(2,JOFF), NPCORD(3,JOFF), & 111 NPFTYP(JOFF) 112 END DO 113 ISTART = ISTART + NPATOM(I) 114 TNPATM = TNPATM + NPATOM(I) 115 IF (TNPATM.GT.MXNPATM) THEN 116 WRITE(LUPRI,'(/2X,A)') & 117 'Maximum number of NP atoms exceeded!' 118 WRITE(LUPRI,'(2X,A,I3,2X,A,I3)') & 119 'Input TNPATM=', TNPATM, & 120 'Maximum allowed:', MXNPATM 121 CALL QUIT('Increase MXNPATM in qmnpmm.h') 122 END IF 123 END DO 124! Convert to atomic units if neeeded 125 IF (UNITS.EQ.'AA') THEN 126 CALL DSCAL(3*TNPATM,XFACT,NPCORD,1) 127 END IF 128! Read force field data 129 READ(LUQMNP,*) FFWORD, TNPFF 130 CALL UPCASE(FFWORD) 131! Check NP force field data consistency 132 IF (FFWORD.NE.'NPFF') THEN 133 WRITE(LUPRI,'(/2X,A,A)') 'Corrupted NP force field ', & 134 'header in POTENTIAL.INP file!' 135 CALL QUIT('Corrupted POTENTIAL.INP') 136 END IF 137 IF (TNPFF.GT.MXNPFF) THEN 138 WRITE(LUPRI,'(/2X,A)') & 139 'Maximum number of NP force field types exceeded!' 140 WRITE(LUPRI,'(2X,A,I3,2X,A,I3)') & 141 'Input TNPFF=', TNPFF, & 142 'Maximum allowed:', MXNPFF 143 CALL QUIT('Increase MXNPFF in qmnpmm.h') 144 END IF 145! Read force field data 146 DO I=1,TNPFF 147 READ(LUQMNP,*) NPFPOL(I), NPFCAP(I), NPFOMG1(I), & 148 NPFGAM1(I), NPFOMG2(I), NPFGAM2(I), & 149 NPFFAC(I) 150 END DO 151! Check force field definitions 152 DO I=1,TNPATM 153 IF (NPFTYP(I).GT.TNPFF) THEN 154 WRITE(LUPRI,'(/2X,A,I4,1X,A)') & 155 'Unknown force field requested for atom', I, & 156 'in NP region!' 157 CALL QUIT('Corrupted POTENTIAL.INP') 158 END IF 159 END DO 160 IF (IPRTLVL.GE.5) THEN 161 CALL PRINT_NPREGION(NPLBL) 162 END IF 163 END IF 164! MM subsystems input 165 IF (DOMMSUB) THEN 166 ISTART = 0 167 DO I=1,TMMBLK 168! Read MM subsystem header 169 READ(LUQMNP,*) MMWORD, MMATOM(I) 170! Check input consistency 171 CALL UPCASE(MMWORD) 172 IF (MMWORD.NE.'MM') THEN 173 WRITE(LUPRI,'(/2X,A,I2,1X,A)') & 174 'Incorrect MM subsystem I=', I, & 175 'header in POTENTIAL.INP file!' 176 CALL QUIT('Corrupted POTENTIAL.INP') 177 END IF 178 IF (MMATOM(I).LE.0) THEN 179 WRITE(LUPRI,'(/2X,A,I5,A,I2,A)') & 180 'Incorrect number of atoms ', MMATOM(I), & 181 'in MM subsystem', I, '!' 182 CALL QUIT('Corrupted POTENTIAL.INP') 183 END IF 184! Read MM subsystem data 185 DO J=1,MMATOM(I) 186 JOFF = ISTART + J 187 READ(LUQMNP,*) MMLBL(JOFF), MMMOL(JOFF), & 188 mm_cord(1,JOFF), mm_cord(2,JOFF), & 189 mm_cord(3,JOFF), MMFTYP(JOFF) 190 END DO 191 ISTART = ISTART + MMATOM(I) 192 TMMATM = TMMATM + MMATOM(I) 193 IF (TMMATM.GT.MXMMATM) THEN 194 WRITE(LUPRI,'(/2X,A)') & 195 'Maximum number of MM atoms exceeded!' 196 WRITE(LUPRI,'(2X,A,I3,2X,A,I3)') & 197 'Input TMMATM=', TMMATM, & 198 'Maximum allowed:', MXMMATM 199 CALL QUIT('Increase MXMMATM in qmnpmm.h') 200 END IF 201 END DO 202! Convert to atomic units if neeeded 203 IF (UNITS.EQ.'AA') THEN 204 CALL DSCAL(3*TMMATM,XFACT,mm_cord,1) 205 END IF 206! Read force field data 207 READ(LUQMNP,*) FFWORD, TMMFF 208 CALL UPCASE(FFWORD) 209! Check NP force field data consistency 210 IF (FFWORD.NE.'MMFF') THEN 211 WRITE(LUPRI,'(/2X,A,A)') 'Corrupted MM force field ', & 212 'header in POTENTIAL.INP file!' 213 CALL QUIT('Corrupted POTENTIAL.INP') 214 END IF 215! Read force field data 216 DO I=1,TMMFF 217 IF ((.NOT.DOMMCAP).AND.(.NOT.DOMMPOL)) THEN 218 READ(LUQMNP,*) MMFM0(I) 219 END IF 220 IF ((.NOT.DOMMCAP).AND.DOMMPOL) THEN 221 READ(LUQMNP,*) MMFM0(I),MMFPOL(I) 222 END IF 223 END DO 224! Check force field definitions 225 DO I=1,TMMATM 226 IF (MMFTYP(I).GT.TMMFF) THEN 227 WRITE(LUPRI,'(/2X,A,I4,1X,A)') & 228 'Unknown force field requested for atom', I, & 229 'in MM region!' 230 CALL QUIT('Corrupted POTENTIAL.INP') 231 END IF 232! Set polarization centers 233 IF (MMFPOL(MMFTYP(I)) .GT. 1.0D-6) THEN 234 TPOLATM = TPOLATM+1 235 MMSKIP(I) = 1 236! works only for one non-metallic MM region 237 END IF 238 END DO 239 IF (IPRTLVL.GE.5) THEN 240 CALL PRINT_MMREGION(MMLBL) 241 END IF 242 END IF 243! Close POTENTIAL.INP file 244 CALL GPCLOSE(LUQMNP,'KEEP') 245! 246! 247 end subroutine 248 SUBROUTINE PRINT_NPREGION(ATMLBL) 249! 250! Purpose: 251! prints detailed information about nanoparticles 252! in QM/NP/MM embedding 253! 254! Input: 255! ATMLBL - List of atom labels for NP region 256! 257! Last updated: 22/03/2013 by Z. Rinkevicius. 258! 259 260#include "priunit.h" 261#include "qmnpmm.h" 262! 263 CHARACTER ATMLBL(MXNPATM)*2 264 integer :: i, j, istart, joff 265! 266 IF (DOMMSUB) THEN 267 CALL TITLER('QM/NP/MM Embedding: Nanoparticle(s) data','*',103) 268 ELSE 269 CALL TITLER('QM/NP Embedding: Nanoparticle(s) data','*',103) 270 END IF 271! 272 IF (TNPBLK.GT.1) THEN 273 WRITE(LUPRI,'(/2X,A,I3,1X,A)') 'Nanoparticle region contains', & 274 TNPBLK, 'separate nanoparticles.' 275 ELSE 276 WRITE(LUPRI,'(/2X,A,A)') 'Nanoparticle region contains single', & 277 ' nanoparticle.' 278 END IF 279 ISTART = 0 280 DO I=1,TNPBLK 281 WRITE(LUPRI, '(/,2X,A,I3,1X,A,F6.3,A)') 'Nanoparticle',I, & 282 'with charge', NPCHRG(I), '. Coordinates in au.' 283 WRITE(LUPRI, '(2X,A)') & 284 '--------------------------------------------------' 285 WRITE(LUPRI, '(2X,A)') & 286 'Atom Coord. X Coord. Y Coord. Z FF Type' 287 WRITE(LUPRI, '(2X,A)') & 288 '==================================================' 289 DO J=1,NPATOM(I) 290 JOFF = ISTART+J 291 WRITE(LUPRI,'(3X,A,1X,F11.5,1X,F11.5,1X,F11.5,3X,I4)') & 292 ATMLBL(JOFF), NPCORD(1,JOFF), NPCORD(2,JOFF), & 293 NPCORD(3,JOFF), NPFTYP(JOFF) 294 END DO 295 WRITE(LUPRI, '(2X,A)') & 296 '==================================================' 297 ISTART = ISTART + NPATOM(I) 298 END DO 299 WRITE(LUPRI,'(/2X,A,A)') 'Static force field(s) data for ', & 300 'nanoparticle region.' 301 WRITE(LUPRI, '(2X,A)') & 302 '-------------------------------------' 303 WRITE(LUPRI, '(2X,A)') & 304 'FF Type Polarizability Capacitance' 305 WRITE(LUPRI, '(2X,A)') & 306 '=====================================' 307 DO I=1,TNPFF 308 WRITE(LUPRI,'(3X,I2,5X,F11.5,5X,F11.5)') I, NPFPOL(I), & 309 NPFCAP(I) 310 END DO 311 WRITE(LUPRI, '(2X,A)') & 312 '=====================================' 313 WRITE(LUPRI,'(/2X,A,A)') 'Dynamic force field(s) data for ', & 314 'nanoparticle region.' 315 WRITE(LUPRI, '(2X,A,A)') & 316 '---------------------------------------------------', & 317 '------' 318 WRITE(LUPRI, '(2X,A,A)') & 319 'FF Type Omega(1) Gamma(1) Omega(2) Gamma(2)', & 320 ' Factor' 321 WRITE(LUPRI, '(2X,A,A)') & 322 '===================================================', & 323 '======' 324 DO I=1,TNPFF 325 WRITE(LUPRI,'(3X,I2,3X,5(1X,F9.5))') I, NPFOMG1(I), NPFGAM1(I),& 326 NPFOMG2(I), NPFGAM2(I), NPFFAC(I) 327 END DO 328 WRITE(LUPRI, '(2X,A,A)') & 329 '===================================================', & 330 '======' 331! 332 end subroutine 333 SUBROUTINE PRINT_MMREGION(ATMLBL) 334! 335! Purpose: 336! prints detailed information about MM region 337! in QM/NP/MM embedding 338! 339! Input: 340! ATMLBL - List of atom labels for MM region 341! 342! Last updated: 22/03/2013 by Z. Rinkevicius. 343! 344#include "priunit.h" 345#include "qmnpmm.h" 346! 347 CHARACTER ATMLBL(MXNPATM)*2 348 integer :: i, istart, j, joff 349! 350 IF (DOMMSUB) THEN 351 CALL TITLER('QM/NP/MM Embedding: MM region(s) data','*',103) 352 END IF 353! 354 IF (TMMBLK.GT.1) THEN 355 WRITE(LUPRI,'(/2X,A,I3,1X,A)') 'MM region contains', & 356 TMMBLK, 'separate non-metallic subregions.' 357 END IF 358 ISTART = 0 359 DO I=1,TMMBLK 360 WRITE(LUPRI, '(/,2X,A,I3,1X,A)') 'MM subregion',I, & 361 '. Coordinates in au.' 362 WRITE(LUPRI, '(2X,A)') & 363 '--------------------------------------------------' 364 WRITE(LUPRI, '(2X,A)') & 365 'Atom Coord. X Coord. Y Coord. Z FF Type' 366 WRITE(LUPRI, '(2X,A)') & 367 '==================================================' 368 DO J=1,MMATOM(I) 369 JOFF = ISTART+J 370 WRITE(LUPRI,'(3X,A,1X,F11.5,1X,F11.5,1X,F11.5,3X,I4)') & 371 ATMLBL(JOFF), mm_cord(1,JOFF), mm_cord(2,JOFF), & 372 mm_cord(3,JOFF), MMFTYP(JOFF) 373 END DO 374 WRITE(LUPRI, '(2X,A)') & 375 '==================================================' 376! works only for one MM region 377 WRITE(LUPRI, '(2X,A,I4)') 'Number of polarizable centers : ',& 378 TPOLATM 379 WRITE(LUPRI, '(2X,A,I4)') 'Number of nonpolarizable centers: ',& 380 MMATOM(I)-TPOLATM 381 ISTART = ISTART + NPATOM(I) 382 END DO 383 WRITE(LUPRI,'(/2X,A)') 'Force field(s) data for MM region' 384! 385 IF ((.NOT.DOMMCAP).AND.(.NOT.DOMMPOL)) THEN 386 WRITE(LUPRI, '(2X,A)') '--------------------------------' 387 WRITE(LUPRI, '(2X,A)') 'FF Type MM Charge Polarizability' 388 WRITE(LUPRI, '(2X,A)') '================================' 389 DO I=1,TMMFF 390 WRITE(LUPRI,'(3X,I2,3X,1X,F9.5,4X,F9.5,1X,I2)') I, MMFM0(I),& 391 MMFPOL(I), MMSKIP(I) 392 END DO 393 WRITE(LUPRI, '(2X,A)') '================================' 394 END IF 395 end subroutine 396 397 398 pure function getdim_relmat(sqflag) 399! 400! Purpose: 401! determines Relay matrix dimensions. 402! 403! Input: 404! SQFLAG - Request total size of Relay matrix 405! Output: 406! IMATDIM - Size of Relay matrix or it's dimension 407! 408! Last updated: 22/03/2013 by Z. Rinkevicius. 409! 410 411 logical, intent(in) :: sqflag 412 integer :: getdim_relmat 413 414#include "qmnpmm.h" 415 416 getdim_relmat = 0 417! Add nanoparticle contribution 418 IF (DONPSUB) THEN 419 getdim_relmat = getdim_relmat + 3*TNPATM 420 IF (DONPCAP) THEN 421 getdim_relmat = getdim_relmat + TNPATM 422 END IF 423 END IF 424! Add Lagrangian term 425 IF (DONPCAP.OR.DOMMCAP) THEN 426 getdim_relmat = getdim_relmat + 1 427 END IF 428! Get requested dimmension 429 IF ((.NOT.MQITER).AND.SQFLAG) THEN 430 getdim_relmat = getdim_relmat*getdim_relmat 431 END IF 432 433 end function 434 435 436 pure subroutine getdim_mmmat(imatdim, sqflag) 437! 438! Purpose: 439! determines Relay matrix dimensions for MM region. 440! 441! Input: 442! SQFLAG - Request total size of Relay matrix 443! Output: 444! IMATDIM - Size of Relay matrix or it's dimension 445! 446! Last updated: 22/03/2013 by Z. Rinkevicius. 447! 448 449 integer, intent(out) :: imatdim 450 logical, intent(in) :: sqflag 451 452#include "qmnpmm.h" 453 454 IMATDIM = 0 455! Add molecular contribution 456 IF (DOMMSUB) THEN 457 IMATDIM = IMATDIM + 3*TPOLATM 458 END IF 459! Get requested dimmension 460 IF ((.NOT.MQITER).AND.SQFLAG) THEN 461 IMATDIM = IMATDIM*IMATDIM 462 END IF 463! 464 end subroutine 465 466 467 SUBROUTINE COMP_RELMAT(FMAT,WORK,LWORK) 468! 469! Purpose: 470! Computes Relay matrix and inverts it or estimates initial 471! MQ vector values for iterative MQ vector determination 472! algorithm. 473! 474! Input: 475! WORK - Temporary memory array 476! LWORK - Size of temporary memory array 477! Output: 478! if .not. MQITER 479! FMAT - Inverted Relay matrix 480! else 481! FMAT - Initial guess of real part of MQ vector 482! 483! Last updated: 22/03/2013 by Z. Rinkevicius. 484! 485 486 487#include "priunit.h" 488#include "qmnpmm.h" 489 490 real(8) :: fmat(:, :, :) 491 integer :: lwork 492 real(8) :: work(lwork) 493 integer, allocatable :: ipiv(:) 494 495 integer :: idimension, ierror 496! 497! Initialize arrays 498 idimension = getdim_relmat(.true.) 499 fmat = 0.0d0 500! Reset matrix dimension parameter 501 idimension = getdim_relmat(.false.) 502 IF (.NOT.MQITER) THEN 503 504! compute polarizabilty dependent terms 505 if (donppol .or. dommpol) then 506 call get_amat(fmat) 507 end if 508 509! compute capacitancy dependent term 510 if (donpcap .or. dommcap) then 511 !todo adapt these routines to handle complex case 512 call get_cmat(fmat) 513 call get_mmat(fmat) 514 call get_qlag(fmat) 515 end if 516 517 ELSE 518! Conjugated gradient method via paradiso solver 519 END IF 520! Print Relay matrix 521 IF ((IPRTLVL.GE.15).AND.(.NOT.MQITER)) THEN 522 WRITE(LUPRI,'(/,2X,A)') '*** Computed Relay matrix ***' 523 CALL OUTPUT(FMAT,1,idimension,1,idimension,idimension,idimension,1,LUPRI) 524 END IF 525! Invert Relay matrix 526 !todo adapt to complex case and invert that one 527 IF (.NOT.MQITER) THEN 528 allocate(ipiv(idimension)) 529 CALL DGETRF(idimension,idimension,FMAT,idimension,IPIV,IERROR) 530 IF (IERROR.NE.0) THEN 531 CALL QUIT('LU factorization failed in COMP_RELMAT') 532 END IF 533 CALL DGETRI(idimension,FMAT,idimension,IPIV,WORK,lwork, & 534 IERROR) 535 IF (IERROR.NE.0) THEN 536 CALL QUIT('Inversion failed in COMP_RELMAT') 537 END IF 538 deallocate(ipiv) 539 IF (IPRTLVL.GE.15) THEN 540 WRITE(LUPRI,'(/,2X,A)') '*** Inverted Relay matrix ***' 541 CALL OUTPUT(FMAT,1,idimension,1,idimension,idimension,idimension,1,LUPRI) 542 END IF 543 END IF 544! 545 end subroutine 546 547 548 SUBROUTINE COMP_MMRELMAT(FMAT,WORK,LWORK) 549! 550! Purpose: 551! Computes Relay matrix and inverts it or estimates initial 552! MQ vector values for iterative MQ vector determination 553! algorithm. (MM region) 554! 555! Input: 556! WORK - Temporary memory array 557! LWORK - Size of temporary memory array 558! Output: 559! if .not. MQITER 560! FMAT - Inverted Relay matrix 561! else 562! FMAT - Initial guess of real part of MQ vector 563! 564! Last updated: 22/03/2013 by Z. Rinkevicius. 565! 566 567#include "priunit.h" 568#include "qmnpmm.h" 569! 570 real(8) :: FMAT(*), WORK(LWORK) 571 integer, allocatable :: ipiv(:) 572 integer :: idimension, ierror, lwork 573! 574! Initialize arrays 575 idimension = getdim_relmat(.true.) 576 CALL DZERO(FMAT,idimension) 577! Reset matrix dimension parameter 578 idimension = getdim_relmat(.false.) 579 IF (.NOT.MQITER) THEN 580! Compute polarizabilty dependent terms 581 CALL GET_MM_AMAT(FMAT,idimension) 582 ELSE 583! Conjugated gradient method via paradiso solver 584 END IF 585! Print Relay matrix 586 IF ((IPRTLVL.GE.15).AND.(.NOT.MQITER)) THEN 587 WRITE(LUPRI,'(/,2X,A)') '*** Computed Relay (MM) matrix ***' 588 CALL OUTPUT(FMAT,1,idimension,1,idimension,idimension,idimension,1,LUPRI) 589 END IF 590! Invert Relay matrix 591 IF (.NOT.MQITER) THEN 592 allocate(ipiv(idimension)) 593 CALL DGETRF(idimension,idimension,FMAT,idimension,IPIV,IERROR) 594 IF (IERROR.NE.0) THEN 595 CALL QUIT('LU factorization failed in COMP_MMRELMAT') 596 END IF 597 CALL DGETRI(idimension,FMAT,idimension,IPIV,WORK,lwork, & 598 IERROR) 599 IF (IERROR.NE.0) THEN 600 CALL QUIT('Inversion failed in COMP_MMRELMAT') 601 END IF 602 deallocate(ipiv) 603 IF (IPRTLVL.GE.15) THEN 604 WRITE(LUPRI,'(/,2X,A)') '*** Inverted Relay(MM) matrix ***' 605 CALL OUTPUT(FMAT,1,idimension,1,idimension,idimension,idimension,1,LUPRI) 606 END IF 607 END IF 608! 609 end subroutine 610 611 612 subroutine get_amat(fmat) 613 ! computes a matrix component of relay matrix for np and mm regions 614 615 real(8), intent(inout) :: fmat(:, :, :) ! relay matrix with m matrix contribution added up 616 617#include "qmnpmm.h" 618 619 integer :: i, j, k, l, m 620 integer :: joff, koff, loff 621 real(8) :: rij(3), rpol, rad, rad2, rad51, rval 622 real(8) :: facta, factb 623 624! Set diagonal components 625 IF (DONPPOL) THEN 626 DO I=1,TNPATM 627 RPOL = 1.0d0/NPFPOL(NPFTYP(I)) 628 DO J=1,3 629 JOFF = (I-1)*3+J 630 FMAT(JOFF,JOFF,1) = RPOL 631 END DO 632 END DO 633 END IF 634! Set off-diagonal components 635 IF (DONPPOL) THEN 636 DO I=1,TNPATM 637 DO J=I+1,TNPATM 638! Compute distance dependent parameters 639 RIJ(1) = NPCORD(1,I)-NPCORD(1,J) 640 RIJ(2) = NPCORD(2,I)-NPCORD(2,J) 641 RIJ(3) = NPCORD(3,I)-NPCORD(3,J) 642 RAD = dsqrt(RIJ(1)*RIJ(1)+RIJ(2)*RIJ(2)+RIJ(3)*RIJ(3)) 643 RAD2 = RAD*RAD 644 RAD51 = 1.0d0/(RAD2*RAD2*RAD) 645 IF (NPMQGAU) CALL GET_GG_AFACT(FACTA,FACTB,I,J,RAD) 646! Distribute interaction tensor 647 DO K=1,3 648 KOFF = (I-1)*3+K 649 DO L=1,3 650 LOFF = (J-1)*3+L 651 RVAL = 3.0d0*RIJ(K)*RIJ(L) 652 IF (K.EQ.L) RVAL = RVAL-RAD2 653 RVAL = RVAL*RAD51 654 IF (NPMQGAU) THEN 655 RVAL = RVAL*FACTA-FACTB*RIJ(K)*RIJ(L) 656 END IF 657 FMAT(KOFF,LOFF,1) = -RVAL 658 FMAT(LOFF,KOFF,1) = -RVAL 659 END DO 660 END DO 661 END DO 662 END DO 663 END IF 664 end subroutine 665 666 667 pure subroutine get_mm_amat(fmat, idimension) 668 ! computes a matrix component of relay matrix for mm region 669 670 671 real(8), intent(inout) :: fmat(idimension, idimension) ! relay matrix with m matrix contribution added up 672 integer, intent(in) :: idimension 673 674#include "qmnpmm.h" 675 676 integer :: istart, i, j, k, l, m 677 integer :: ioff, joff, koff, loff 678 real(8) :: rij(3), rad, rad3, rval, fact 679 real(8) :: rpol, rad2, rad51 680 681! Set diagonal components 682 IOFF = 1 683 DO I=1,TMMATM 684 IF (MMSKIP(I) .EQ. 0) CYCLE 685 RPOL = 1.0d0/MMFPOL(MMFTYP(I)) 686 DO J=1,3 687 JOFF = (IOFF-1)*3+J 688 FMAT(JOFF,JOFF) = RPOL 689 END DO 690 IOFF = IOFF+1 691 END DO 692! Set off-diagonal components 693 IOFF = 1 694 DO I=1,TMMATM 695 IF (MMSKIP(I) .EQ. 0) CYCLE 696 JOFF = IOFF+1 697 DO J=I+1,TMMATM 698 IF (MMSKIP(J) .EQ. 0) CYCLE 699! Compute distance dependent parameters 700 RIJ(1) = mm_cord(1,I)-mm_cord(1,J) 701 RIJ(2) = mm_cord(2,I)-mm_cord(2,J) 702 RIJ(3) = mm_cord(3,I)-mm_cord(3,J) 703 RAD = dsqrt(RIJ(1)*RIJ(1)+RIJ(2)*RIJ(2)+RIJ(3)*RIJ(3)) 704 RAD2 = RAD*RAD 705 RAD51 = 1.0d0/(RAD2*RAD2*RAD) 706! Distribute interaction tensor 707 DO K=1,3 708 KOFF = (IOFF-1)*3+K 709 DO L=1,3 710 LOFF = (JOFF-1)*3+L 711 RVAL = 3.0d0*RIJ(K)*RIJ(L) 712 IF (K.EQ.L) RVAL = RVAL-RAD2 713 RVAL = RVAL*RAD51 714 FMAT(KOFF,LOFF) = -RVAL 715 FMAT(LOFF,KOFF) = -RVAL 716 END DO 717 END DO 718 JOFF = JOFF+1 719 END DO 720 IOFF = IOFF+1 721 END DO 722! 723 end subroutine 724 725 726 subroutine get_gg_afact(facta, factb, iatm, jatm, rad) 727! 728! Purpose: 729! Determines damping factors for T(2) operator for 730! Gaussian/Gaussian dipole model. 731! 732! Input: 733! IATM - I-th atom in NP region 734! JATM - J-th atom in NP region 735! RAD - Distance between I and J atoms 736! Output: 737! FACTA - Scalling factor 738! FACTB - Additional factor 739! 740! Last updated: 22/03/2013 by Z. Rinkevicius. 741! 742 743 real(8), intent(out) :: facta 744 real(8), intent(out) :: factb 745 integer, intent(in) :: iatm 746 integer, intent(in) :: jatm 747 real(8), intent(in) :: rad 748 749#include "pi.h" 750#include "qmnpmm.h" 751 752 real(8) :: ripol 753 real(8) :: rjpol 754 real(8) :: rdim 755 real(8) :: rim 756 real(8) :: rjm 757 real(8) :: rijm 758 real(8) :: rval 759 real(8),external :: erf 760! Get polarizabilities 761 RIPOL = NPFPOL(NPFTYP(IATM))/3.0d0 762 RJPOL = NPFPOL(NPFTYP(JATM))/3.0d0 763! Get damping radius 764 RDIM = dsqrt(2.0d0)/SQRTPI 765 RIM = (RDIM*RIPOL)**(1.0d0/3.0d0) 766 RJM = (RDIM*RJPOL)**(1.0d0/3.0d0) 767 RIJM = dsqrt(RIM*RIM+RJM*RJM) 768 RVAL = RAD/RIJM 769! Compute factors 770 FACTA = ERF(RVAL)-2.0d0*RVAL*DEXP(-RVAL*RVAL)/SQRTPI 771 FACTB = 4.0d0*DEXP(-RVAL*RVAL) 772 FACTB = FACTB/(RAD*RAD*RIJM*RIJM*RIJM*SQRTPI) 773! 774 end subroutine 775 776 777 subroutine get_cmat(fmat) 778 ! computes c matrix component of relay matrix for np and mm regions 779 780 real(8), intent(inout) :: fmat(:, :, :) ! relay matrix with m matrix contribution added up 781 782#include "qmnpmm.h" 783 784 integer :: istart, i, j, m 785 integer :: ioff, joff 786 real(8) :: rij(3), rad, rad3, rval, fact 787 788! Set diagonal components 789 IF (DONPCAP) THEN 790 ISTART = 0 791 IF (DONPPOL) ISTART = ISTART+3*TNPATM 792! Fix me: MM shift 793 DO I=1,TNPATM 794 IOFF = ISTART+I 795 FMAT(IOFF,IOFF,1) = -1.0d0/NPFCAP(NPFTYP(I)) 796 END DO 797 END IF 798! Set off-diagonal components 799 IF (DONPCAP) THEN 800 ISTART = 0 801 IF (DONPPOL) ISTART = ISTART+3*TNPATM 802! Fix me: MM shift 803 DO I=1,TNPATM 804 IOFF = ISTART+I 805 DO J=I+1,TNPATM 806 JOFF = ISTART+J 807! Compute distance dependent parameters 808 RIJ(1) = NPCORD(1,I)-NPCORD(1,J) 809 RIJ(2) = NPCORD(2,I)-NPCORD(2,J) 810 RIJ(3) = NPCORD(3,I)-NPCORD(3,J) 811 RAD = dsqrt(RIJ(1)*RIJ(1)+RIJ(2)*RIJ(2)+RIJ(3)*RIJ(3)) 812 RVAL = 1.0d0/RAD 813 IF (NPMQGAU) THEN 814 CALL GET_GG_CFACT(FACT,I,J,RAD) 815 RVAL = FACT*RVAL 816 END IF 817 FMAT(IOFF,JOFF,1) = -RVAL 818 FMAT(JOFF,IOFF,1) = -RVAL 819 END DO 820 END DO 821 END IF 822 end subroutine 823 824 825 subroutine get_gg_cfact(fact, iatm, jatm, rad) 826! 827! Purpose: 828! Determines damping factors for T(0) operator for 829! Gaussian/Gaussian dipole model. 830! 831! Input: 832! IATM - I-th atom in NP region 833! JATM - J-th atom in NP region 834! RAD - Distance between I and J atoms 835! Output: 836! FACT - Scalling factor 837! 838! Last updated: 22/03/2013 by Z. Rinkevicius. 839! 840 841 real(8), intent(out) :: fact 842 integer, intent(in) :: iatm 843 integer, intent(in) :: jatm 844 real(8), intent(in) :: rad 845 846#include "pi.h" 847#include "qmnpmm.h" 848 849 real(8) :: ricap 850 real(8) :: rjcap 851 real(8) :: rijm 852 real(8),external :: erf 853 854! Get capacitancies 855 RICAP = 1.41421356237309504880D0*NPFCAP(NPFTYP(IATM))/SQRTPI 856 RJCAP = 1.41421356237309504880D0*NPFCAP(NPFTYP(JATM))/SQRTPI 857 858! Get damping radius & scalling factor 859 RIJM = dsqrt(RICAP*RICAP+RJCAP*RJCAP) 860 FACT = ERF(RAD/RIJM) 861 862 end subroutine 863 864 865 subroutine get_mmat(fmat) 866 ! computes m matrix component of relay matrix for np and mm regions 867 868 real(8), intent(inout) :: fmat(:, :, :) ! relay matrix with m matrix contribution added up 869 870#include "qmnpmm.h" 871 872 integer :: istart, i, j, m 873 integer :: ioff, joff 874 real(8) :: rij(3), rad, rad3, rval, fact 875 876! Set off-diagonal components 877 IF (DONPPOL.AND.DONPCAP) THEN 878 ISTART = 0 879 IF (DONPPOL) ISTART = ISTART+3*TNPATM 880 DO I=1,TNPATM 881 DO J=1,TNPATM 882 IF (I.EQ.J) CYCLE 883! Compute distance dependent parameters 884 RIJ(1) = NPCORD(1,I)-NPCORD(1,J) 885 RIJ(2) = NPCORD(2,I)-NPCORD(2,J) 886 RIJ(3) = NPCORD(3,I)-NPCORD(3,J) 887 RAD = dsqrt(RIJ(1)*RIJ(1)+RIJ(2)*RIJ(2)+RIJ(3)*RIJ(3)) 888 RAD3 = 1.0d0/(RAD*RAD*RAD) 889! Get damping factor 890 IF (NPMQGAU) CALL GET_GG_MFACT(FACT,I,J,RAD) 891! Distribute contributions 892 DO M=1,3 893 IOFF = (I-1)*3+M 894 JOFF = ISTART+J 895 RVAL = RIJ(M)*RAD3 896 IF (NPMQGAU) RVAL = RVAL*FACT 897 FMAT(IOFF,JOFF,1) = -RVAL 898 FMAT(JOFF,IOFF,1) = -RVAL 899 END DO 900 END DO 901 END DO 902 END IF 903 904 end subroutine 905 906 907 subroutine get_gg_mfact(fact, iatm, jatm, distance_i_j) 908 ! determines damping factors for t(1) operator for 909 ! gaussian/gaussian dipole model 910 911 912 real(8), intent(inout) :: fact ! scaling factor 913 integer, intent(in) :: iatm ! i-th atom in np region 914 integer, intent(in) :: jatm ! j-th atom in np region 915 real(8), intent(in) :: distance_i_j 916 917! sqrtpi 918#include "pi.h" 919 920! npftyp, npfpol, npfcap 921#include "qmnpmm.h" 922 923 real(8) :: ripol 924 real(8) :: rdim 925 real(8) :: rim 926 real(8) :: rjcap 927 real(8) :: rijm 928 real(8) :: radx 929 real(8),external :: erf 930 931! get polarizabilities 932 ripol = npfpol(npftyp(iatm))/3.0d0 933! get damping radius 934 rdim = 1.41421356237309504880d0/sqrtpi 935 rim = (rdim*ripol)**(1.0d0/3.0d0) 936! get j-th capacitancy 937 rjcap = rdim*npfcap(npftyp(jatm)) 938! get damping radius & scalling factor 939 rijm = dsqrt(rim*rim+rjcap*rjcap) 940 radx = distance_i_j/rijm 941 fact = erf(radx)-2.0d0*radx*dexp(-radx*radx)/sqrtpi 942 943 end subroutine 944 945 946 pure subroutine get_qlag(fmat) 947 ! sets charge constrain in relay matrix for np and mm regions 948 949 real(8), intent(inout) :: fmat(:, :, :) 950 951! donpcap, dommcap, donppol 952#include "qmnpmm.h" 953 954 integer :: i, istart, ioff, idimension 955 idimension = size(fmat, 1) 956 957 ! set diagonal components 958 if (donpcap .or. dommcap) then 959 istart = 0 960 if (donppol) istart = istart + 3*tnpatm 961! fixme: mm shift 962 if (donpcap) then 963 do i = 1, tnpatm 964 ioff = istart + i 965 fmat(ioff, idimension, 1) = 1.0d0 966 fmat(idimension, ioff, 1) = 1.0d0 967 end do 968 end if 969 end if 970 971 end subroutine 972 973 974 975 SUBROUTINE WRITE_RELMAT(FMAT) 976! 977! Purpose: 978! Stores Relay matrix in binary file. 979! 980! Input: 981! FMAT - Inverted Relay matrix 982! 983! Last updated: 16/08/2013 by Z. Rinkevicius. 984! 985 986#include "qmnpmm.h" 987#include "dummy.h" 988#include "iratdef.h" 989#include "inftap.h" 990! 991 real(8) :: fmat(*) 992 integer :: idimension 993 integer :: luqmnp 994 995! determine dimensions 996 idimension = getdim_relmat(.true.) 997! write inverted Relay matrix 998 LUQMNP = -1 999 CALL GPOPEN(LUQMNP,'QMMMNP','UNKNOWN','SEQUENTIAL','UNFORMATTED', & 1000 IDUMMY,.FALSE.) 1001 REWIND(LUQMNP) 1002 CALL WRTIEF(FMAT, idimension, 'QQMNPMAT', LUQMNP) 1003 CALL GPCLOSE(LUQMNP,'KEEP') 1004! 1005 end subroutine 1006 SUBROUTINE READ_RELMAT(FMAT) 1007! 1008! Purpose: 1009! Reads Relay matrix in binary file. 1010! 1011! Input: 1012! FMAT - Inverted Relay matrix 1013! 1014! Last updated: 16/08/2013 by Z. Rinkevicius. 1015! 1016 1017#include "qmnpmm.h" 1018#include "dummy.h" 1019#include "iratdef.h" 1020#include "inftap.h" 1021! 1022 real(8) :: fmat(*) 1023 integer :: idimension 1024 integer :: luqmnp 1025 logical :: fndlab 1026 1027! determine dimensions 1028 idimension = getdim_relmat(.true.) 1029 1030! read inverted Relay matrix 1031 LUQMNP=-1 1032 CALL GPOPEN(LUQMNP,'QMMMNP','UNKNOWN','SEQUENTIAL','UNFORMATTED', & 1033 IDUMMY,.FALSE.) 1034 REWIND(LUQMNP) 1035 IF (FNDLAB('QQMNPMAT',LUQMNP)) THEN 1036 CALL READT(LUQMNP,idimension,FMAT) 1037 ELSE 1038 CALL QUIT('Problem reading the matrix from the QMMMNP file.') 1039 ENDIF 1040 CALL GPCLOSE(LUQMNP,'KEEP') 1041! 1042 end subroutine 1043 1044 1045 1046 1047 subroutine comp_dampvmat(fmat, mqvec) 1048 ! computes damped potential matrix in ao basis 1049 1050 1051 real(8), intent(in) :: fmat(*) ! inverted relay matrix 1052 real(8), intent(inout) :: mqvec(*) 1053 1054#include "priunit.h" 1055#include "qmnpmm.h" 1056#include "maxorb.h" 1057#include "aovec.h" 1058#include "shells.h" 1059#include "primit.h" 1060 1061#ifdef VAR_MPI 1062! qmcmm_work 1063#include "iprtyp.h" 1064#endif 1065 1066 real(8), allocatable :: rdvec(:) 1067 real(8), allocatable :: rqvec(:) 1068 integer :: idimension, i 1069 integer :: iprint 1070 1071 idimension = getdim_relmat(.false.) 1072 1073 ! allocate and set gaussian broadening paramenters 1074 1075 allocate(rdvec(tnpatm)) 1076 allocate(rqvec(tnpatm)) 1077 1078 call set_damparam(rdvec, rqvec) 1079 1080 if (iprtlvl > 14) then 1081 write(lupri, '(/,2x,a)') '*** Computed MQ vector start ***' 1082 do i = 1, idimension 1083 write(lupri, '(i8, f18.8)') i, mqvec(i) 1084 end do 1085 write(lupri, '(/,2x,a)') '*** Computed MQ vector end ***' 1086 end if 1087 1088#ifdef VAR_MPI 1089 call mpixbcast(QMCMM_WORK, 1, 'INTEGER', 0) 1090 iprint = 0 1091 call mpixbcast(iprint, 1, 'INTEGER', 0) 1092#endif 1093 1094#ifdef ENABLE_VPOTDAMP 1095 call vpotdamped(kmax, nhkt, nuco, nrco, jstrt, cent, priccf, priexp, & 1096 fmat, & 1097 npcord, & 1098 mqvec(3*tnpatm+1), rqvec, & 1099 mqvec , rdvec, & 1100 tnpatm) 1101#else 1102 call quit('VPOTDAMP not compiled in this version') 1103#endif 1104 1105 deallocate(rdvec) 1106 deallocate(rqvec) 1107 1108 end subroutine 1109 1110 1111 pure subroutine set_damparam(rdvec, rqvec) 1112 ! sets damping parameters vectors for dipoles and charges 1113 ! see Eqs 14 and 15 in J. Phys. Chem. C 112, 40 (2008) 1114 1115 1116 real(8), intent(inout) :: rdvec(*) 1117 real(8), intent(inout) :: rqvec(*) 1118 1119! tnpatm, npfpol, npfcap 1120#include "qmnpmm.h" 1121 1122! sqrtpi 1123#include "pi.h" 1124 1125 real(8) :: rdim 1126 real(8) :: ripol 1127 integer :: i 1128 1129 rdim = dsqrt(2.0d0)/sqrtpi 1130 do i = 1, tnpatm 1131 ripol = npfpol(npftyp(i))/3.0d0 1132 rdvec(i) = (rdim*ripol)**(1.0d0/3.0d0) 1133 rqvec(i) = rdim*npfcap(npftyp(i)) 1134 end do 1135 1136 end subroutine 1137 1138end module 1139