1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Set of routines to: 8!> Contract integrals over primitive Gaussians 9!> Decontract (density) matrices 10!> Trace matrices to get forces 11!> Block copy and add matrices 12!> \par History 13!> none 14!> \author JGH (01.07.2014) 15! ************************************************************************************************** 16MODULE ai_contraction 17 18 USE kinds, ONLY: dp 19#include "../base/base_uses.f90" 20 21 IMPLICIT NONE 22 23 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_contraction' 24 25 PRIVATE 26 27 PUBLIC :: contraction, decontraction, block_add, force_trace 28 29 INTERFACE contraction 30 MODULE PROCEDURE contraction_ab, contraction_abc 31 END INTERFACE 32 33 INTERFACE decontraction 34 MODULE PROCEDURE decontraction_ab 35 END INTERFACE 36 37 INTERFACE force_trace 38 MODULE PROCEDURE force_trace_ab 39 END INTERFACE 40 41 INTERFACE block_add 42 MODULE PROCEDURE block_add_ab 43 END INTERFACE 44 45! ************************************************************************************************** 46 47CONTAINS 48 49! ************************************************************************************************** 50!> \brief Applying the contraction coefficents to a set of two-center primitive 51!> integrals 52!> QAB <- CA(T) * SAB * CB 53!> QAB is optionally scaled with "fscale" 54!> Variable "trans" requests the output to be QAB(T) 55!> If only one of the transformation matrix is given, only a half 56!> transformation is done 57!> Active dimensions are: QAB(ma,mb), SAB(na,nb) 58!> \param sab Input matrix, dimension(:,:) 59!> \param qab Output matrix, dimension(:,:) 60!> \param ca Left transformation matrix, optional 61!> \param na First dimension of ca, optional 62!> \param ma Second dimension of ca, optional 63!> \param cb Right transformation matrix, optional 64!> \param nb First dimension of cb, optional 65!> \param mb Second dimension of cb, optional 66!> \param fscale Optional scaling of output 67!> \param trans Optional transposition of output 68! ************************************************************************************************** 69 SUBROUTINE contraction_ab(sab, qab, ca, na, ma, cb, nb, mb, fscale, trans) 70 71 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sab 72 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: qab 73 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 74 OPTIONAL :: ca 75 INTEGER, INTENT(IN), OPTIONAL :: na, ma 76 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 77 OPTIONAL :: cb 78 INTEGER, INTENT(IN), OPTIONAL :: nb, mb 79 REAL(KIND=dp), INTENT(IN), OPTIONAL :: fscale 80 LOGICAL, INTENT(IN), OPTIONAL :: trans 81 82 CHARACTER(len=*), PARAMETER :: routineN = 'contraction_ab', routineP = moduleN//':'//routineN 83 84 INTEGER :: lda, ldb, ldq, lds, ldw, mal, mbl, nal, & 85 nbl 86 LOGICAL :: my_trans 87 REAL(KIND=dp) :: fs 88 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work 89 90! Should output matrix be transposed? 91 92 IF (PRESENT(trans)) THEN 93 my_trans = trans 94 ELSE 95 my_trans = .FALSE. 96 END IF 97 98 ! Scaling of output matrix 99 IF (PRESENT(fscale)) THEN 100 fs = fscale 101 ELSE 102 fs = 1.0_dp 103 END IF 104 105 ! Active matrix size 106 IF (PRESENT(ca)) THEN 107 IF (PRESENT(na)) THEN 108 nal = na 109 ELSE 110 nal = SIZE(ca, 1) 111 END IF 112 IF (PRESENT(ma)) THEN 113 mal = ma 114 ELSE 115 mal = SIZE(ca, 2) 116 END IF 117 lda = SIZE(ca, 1) 118 END IF 119 IF (PRESENT(cb)) THEN 120 IF (PRESENT(nb)) THEN 121 nbl = nb 122 ELSE 123 nbl = SIZE(cb, 1) 124 END IF 125 IF (PRESENT(mb)) THEN 126 mbl = mb 127 ELSE 128 mbl = SIZE(cb, 2) 129 END IF 130 ldb = SIZE(cb, 1) 131 END IF 132 133 lds = SIZE(sab, 1) 134 ldq = SIZE(qab, 1) 135 136 IF (PRESENT(ca) .AND. PRESENT(cb)) THEN 137 ! Full transform 138 ALLOCATE (work(nal, mbl)) 139 ldw = nal 140 CALL dgemm("N", "N", nal, mbl, nbl, 1.0_dp, sab(1, 1), lds, cb(1, 1), ldb, 0.0_dp, work(1, 1), ldw) 141 IF (my_trans) THEN 142 CALL dgemm("T", "N", mbl, mal, nal, fs, work(1, 1), ldw, ca(1, 1), lda, 0.0_dp, qab(1, 1), ldq) 143 ELSE 144 CALL dgemm("T", "N", mal, mbl, nal, fs, ca(1, 1), lda, work(1, 1), ldw, 0.0_dp, qab(1, 1), ldq) 145 END IF 146 DEALLOCATE (work) 147 ELSE IF (PRESENT(ca)) THEN 148 IF (PRESENT(nb)) THEN 149 nbl = nb 150 ELSE 151 nbl = SIZE(sab, 2) 152 END IF 153 IF (my_trans) THEN 154 CALL dgemm("T", "N", nbl, mal, nal, fs, sab(1, 1), lds, ca(1, 1), lda, 0.0_dp, qab(1, 1), ldq) 155 ELSE 156 CALL dgemm("T", "N", mal, nbl, nal, fs, ca(1, 1), lda, sab(1, 1), lds, 0.0_dp, qab(1, 1), ldq) 157 END IF 158 ELSE IF (PRESENT(cb)) THEN 159 IF (PRESENT(na)) THEN 160 nal = na 161 ELSE 162 nal = SIZE(sab, 1) 163 END IF 164 IF (my_trans) THEN 165 CALL dgemm("N", "N", nal, mbl, nbl, fs, sab(1, 1), lds, cb(1, 1), ldb, 0.0_dp, qab, ldq) 166 ELSE 167 CALL dgemm("T", "T", mbl, nal, nbl, fs, cb(1, 1), ldb, sab(1, 1), lds, 0.0_dp, qab, ldq) 168 END IF 169 ELSE 170 ! Copy of arrays is not covered here 171 CPABORT("Copy of arrays is not covered here") 172 END IF 173 174 END SUBROUTINE contraction_ab 175 176! ************************************************************************************************** 177!> \brief Applying the contraction coefficents to a tripple set integrals 178!> QABC <- CA(T) * SABC * CB * CC 179!> If only one or two of the transformation matrices are given, only a 180!> part transformation is done 181!> \param sabc Input matrix, dimension(:,:) 182!> \param qabc Output matrix, dimension(:,:) 183!> \param ca Transformation matrix (index 1), optional 184!> \param na First dimension of ca, optional 185!> \param ma Second dimension of ca, optional 186!> \param cb Transformation matrix (index 2), optional 187!> \param nb First dimension of cb, optional 188!> \param mb Second dimension of cb, optional 189!> \param cc Transformation matrix (index 3), optional 190!> \param nc First dimension of cc, optional 191!> \param mc Second dimension of cc, optional 192! ************************************************************************************************** 193 SUBROUTINE contraction_abc(sabc, qabc, ca, na, ma, cb, nb, mb, cc, nc, mc) 194 195 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: sabc 196 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: qabc 197 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 198 OPTIONAL :: ca 199 INTEGER, INTENT(IN), OPTIONAL :: na, ma 200 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 201 OPTIONAL :: cb 202 INTEGER, INTENT(IN), OPTIONAL :: nb, mb 203 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 204 OPTIONAL :: cc 205 INTEGER, INTENT(IN), OPTIONAL :: nc, mc 206 207 CHARACTER(len=*), PARAMETER :: routineN = 'contraction_abc', & 208 routineP = moduleN//':'//routineN 209 210 INTEGER :: lda, ldb, ldc, mal, mbl, mcl, nal, nbl, & 211 ncl 212 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work1, work2, work3, work4 213 214! Active matrix size 215 216 IF (PRESENT(ca)) THEN 217 IF (PRESENT(na)) THEN 218 nal = na 219 ELSE 220 nal = SIZE(ca, 1) 221 END IF 222 IF (PRESENT(ma)) THEN 223 mal = ma 224 ELSE 225 mal = SIZE(ca, 2) 226 END IF 227 lda = SIZE(ca, 1) 228 END IF 229 IF (PRESENT(cb)) THEN 230 IF (PRESENT(nb)) THEN 231 nbl = nb 232 ELSE 233 nbl = SIZE(cb, 1) 234 END IF 235 IF (PRESENT(mb)) THEN 236 mbl = mb 237 ELSE 238 mbl = SIZE(cb, 2) 239 END IF 240 ldb = SIZE(cb, 1) 241 END IF 242 IF (PRESENT(cc)) THEN 243 IF (PRESENT(nc)) THEN 244 ncl = nc 245 ELSE 246 ncl = SIZE(cc, 1) 247 END IF 248 IF (PRESENT(mc)) THEN 249 mcl = mc 250 ELSE 251 mcl = SIZE(cc, 2) 252 END IF 253 ldc = SIZE(cc, 1) 254 END IF 255 256 IF (PRESENT(ca) .AND. PRESENT(cb) .AND. PRESENT(cc)) THEN 257 ! Full transform 258 ALLOCATE (work1(nal, nbl, ncl)) 259 ! make sure that we have contigous memory, needed for transpose algorithm 260 work1(1:nal, 1:nbl, 1:ncl) = sabc(1:nal, 1:nbl, 1:ncl) 261 ! 262 ALLOCATE (work2(nbl, ncl, mal)) 263 CALL dgemm("T", "N", nbl*ncl, mal, nal, 1.0_dp, work1(1, 1, 1), nal, ca(1, 1), lda, 0.0_dp, work2(1, 1, 1), nbl*ncl) 264 ! 265 ALLOCATE (work3(ncl, mal, mbl)) 266 CALL dgemm("T", "N", ncl*mal, mbl, nbl, 1.0_dp, work2(1, 1, 1), nbl, cb(1, 1), ldb, 0.0_dp, work3(1, 1, 1), ncl*mal) 267 ! 268 ALLOCATE (work4(mal, mbl, mcl)) 269 CALL dgemm("T", "N", mal*mbl, mcl, ncl, 1.0_dp, work3(1, 1, 1), ncl, cc(1, 1), ldc, 0.0_dp, work4(1, 1, 1), mal*mbl) 270 ! 271 work4(1:mal, 1:mbl, 1:mcl) = qabc(1:mal, 1:mbl, 1:mcl) 272 ! 273 DEALLOCATE (work1, work2, work3, work4) 274 ! 275 ELSE IF (PRESENT(ca) .AND. PRESENT(cb)) THEN 276 CPABORT("Not implemented") 277 ELSE IF (PRESENT(ca) .AND. PRESENT(cc)) THEN 278 CPABORT("Not implemented") 279 ELSE IF (PRESENT(cb) .AND. PRESENT(cc)) THEN 280 CPABORT("Not implemented") 281 ELSE IF (PRESENT(ca)) THEN 282 CPABORT("Not implemented") 283 ELSE IF (PRESENT(cb)) THEN 284 CPABORT("Not implemented") 285 ELSE IF (PRESENT(cc)) THEN 286 CPABORT("Not implemented") 287 ELSE 288 ! Copy of arrays is not covered here 289 CPABORT("Copy of arrays is not covered here") 290 END IF 291 292 END SUBROUTINE contraction_abc 293 294! ************************************************************************************************** 295!> \brief Applying the de-contraction coefficents to a matrix 296!> QAB <- CA * SAB * CB(T) 297!> Variable "trans" requests the input matrix to be SAB(T) 298!> Active dimensions are: QAB(na,nb), SAB(ma,mb) 299!> \param sab Input matrix, dimension(:,:) 300!> \param qab Output matrix, dimension(:,:) 301!> \param ca Left transformation matrix 302!> \param na First dimension of ca 303!> \param ma Second dimension of ca 304!> \param cb Right transformation matrix 305!> \param nb First dimension of cb 306!> \param mb Second dimension of cb 307!> \param trans Optional transposition of input matrix 308! ************************************************************************************************** 309 SUBROUTINE decontraction_ab(sab, qab, ca, na, ma, cb, nb, mb, trans) 310 311 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sab 312 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: qab 313 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ca 314 INTEGER, INTENT(IN) :: na, ma 315 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cb 316 INTEGER, INTENT(IN) :: nb, mb 317 LOGICAL, INTENT(IN), OPTIONAL :: trans 318 319 CHARACTER(len=*), PARAMETER :: routineN = 'decontraction_ab', & 320 routineP = moduleN//':'//routineN 321 322 INTEGER :: lda, ldb, ldq, lds, ldw 323 LOGICAL :: my_trans 324 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work 325 326! Should input matrix be transposed? 327 328 IF (PRESENT(trans)) THEN 329 my_trans = trans 330 ELSE 331 my_trans = .FALSE. 332 END IF 333 334 lds = SIZE(sab, 1) 335 ldq = SIZE(qab, 1) 336 lda = SIZE(ca, 1) 337 ldb = SIZE(cb, 1) 338 339 ALLOCATE (work(na, mb)) 340 ldw = na 341 342 IF (my_trans) THEN 343 CALL dgemm("N", "T", na, mb, ma, 1.0_dp, ca, lda, sab, lds, 0.0_dp, work, ldw) 344 ELSE 345 CALL dgemm("N", "N", na, mb, ma, 1.0_dp, ca, lda, sab, lds, 0.0_dp, work, ldw) 346 END IF 347 CALL dgemm("N", "T", na, nb, mb, 1.0_dp, work, ldw, cb, ldb, 0.0_dp, qab, ldq) 348 349 DEALLOCATE (work) 350 351 END SUBROUTINE decontraction_ab 352 353! ************************************************************************************************** 354!> \brief Routine to trace a series of matrices with another matrix 355!> Calculate forces of type f(:) = Trace(Pab*Sab(:)) 356!> \param force Vector to hold output forces 357!> \param sab Input vector of matrices, dimension (:,:,:) 358!> \param pab Input matrix 359!> \param na Active first dimension 360!> \param nb Active second dimension 361!> \param m Number of matrices to be traced 362!> \param trans Matrices are transposed (Sab and Pab) 363! ************************************************************************************************** 364 SUBROUTINE force_trace_ab(force, sab, pab, na, nb, m, trans) 365 366 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: force 367 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: sab 368 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pab 369 INTEGER, INTENT(IN) :: na, nb, m 370 LOGICAL, INTENT(IN), OPTIONAL :: trans 371 372 CHARACTER(len=*), PARAMETER :: routineN = 'force_trace_ab', routineP = moduleN//':'//routineN 373 374 INTEGER :: i 375 LOGICAL :: my_trans 376 377 CPASSERT(m <= SIZE(SAB, 3)) 378 CPASSERT(m <= SIZE(force, 1)) 379 380 ! are matrices transposed? 381 IF (PRESENT(trans)) THEN 382 my_trans = trans 383 ELSE 384 my_trans = .FALSE. 385 END IF 386 387 DO i = 1, m 388 IF (my_trans) THEN 389 force(i) = SUM(sab(1:nb, 1:na, i)*pab(1:nb, 1:na)) 390 ELSE 391 force(i) = SUM(sab(1:na, 1:nb, i)*pab(1:na, 1:nb)) 392 END IF 393 END DO 394 395 END SUBROUTINE force_trace_ab 396 397! ************************************************************************************************** 398!> \brief Copy a block out of a matrix and add it to another matrix 399!> SAB = SAB + QAB or QAB = QAB + SAB 400!> QAB(ia:,ib:) and SAB(1:,1:) 401!> \param dir "IN" and "OUT" defines direction of copy 402!> \param sab Matrix input for "IN", output for "OUT" 403!> \param na first dimension of matrix to copy 404!> \param nb second dimension of matrix to copy 405!> \param qab Matrix output for "IN", input for "OUT" 406!> Use subblock of this matrix 407!> \param ia Starting index in qab first dimension 408!> \param ib Starting index in qab second dimension 409!> \param trans Matrices (qab and sab) are transposed 410! ************************************************************************************************** 411 SUBROUTINE block_add_ab(dir, sab, na, nb, qab, ia, ib, trans) 412 413 CHARACTER(LEN=*), INTENT(IN) :: dir 414 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab 415 INTEGER, INTENT(IN) :: na, nb 416 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: qab 417 INTEGER, INTENT(IN) :: ia, ib 418 LOGICAL, INTENT(IN), OPTIONAL :: trans 419 420 CHARACTER(len=*), PARAMETER :: routineN = 'block_add_ab', routineP = moduleN//':'//routineN 421 422 INTEGER :: ja, jb 423 LOGICAL :: my_trans 424 425 IF (PRESENT(trans)) THEN 426 my_trans = trans 427 ELSE 428 my_trans = .FALSE. 429 END IF 430 431 IF (dir == "IN" .OR. dir == "in") THEN 432 ! QAB(block) <= SAB 433 ja = ia + na - 1 434 jb = ib + nb - 1 435 IF (my_trans) THEN 436 qab(ib:jb, ia:ja) = qab(ib:jb, ia:ja) + sab(1:nb, 1:na) 437 ELSE 438 qab(ia:ja, ib:jb) = qab(ia:ja, ib:jb) + sab(1:na, 1:nb) 439 END IF 440 ELSEIF (dir == "OUT" .OR. dir == "out") THEN 441 ! SAB <= QAB(block) 442 ja = ia + na - 1 443 jb = ib + nb - 1 444 IF (my_trans) THEN 445 sab(1:nb, 1:na) = sab(1:nb, 1:na) + qab(ib:jb, ia:ja) 446 ELSE 447 sab(1:na, 1:nb) = sab(1:na, 1:nb) + qab(ia:ja, ib:jb) 448 END IF 449 ELSE 450 CPABORT("") 451 END IF 452 453 END SUBROUTINE block_add_ab 454! ************************************************************************************************** 455 456END MODULE ai_contraction 457