LCOV - code coverage report
Current view: top level - src - pao_ml_descriptor.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.8 % 184 180
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 4 4

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Feature vectors for describing chemical environments in a rotationally invariant fashion.
      10              : !> \author Ole Schuett
      11              : ! **************************************************************************************************
      12              : MODULE pao_ml_descriptor
      13              :    USE atomic_kind_types,               ONLY: get_atomic_kind
      14              :    USE basis_set_types,                 ONLY: gto_basis_set_type
      15              :    USE cell_types,                      ONLY: cell_type,&
      16              :                                               pbc
      17              :    USE kinds,                           ONLY: dp
      18              :    USE mathconstants,                   ONLY: fourpi,&
      19              :                                               rootpi
      20              :    USE mathlib,                         ONLY: diamat_all
      21              :    USE pao_input,                       ONLY: pao_ml_desc_overlap,&
      22              :                                               pao_ml_desc_pot,&
      23              :                                               pao_ml_desc_r12
      24              :    USE pao_potentials,                  ONLY: pao_calc_gaussian
      25              :    USE pao_types,                       ONLY: pao_env_type
      26              :    USE particle_types,                  ONLY: particle_type
      27              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      28              :                                               pao_descriptor_type,&
      29              :                                               qs_kind_type
      30              :    USE util,                            ONLY: sort
      31              : #include "./base/base_uses.f90"
      32              : 
      33              :    IMPLICIT NONE
      34              : 
      35              :    PRIVATE
      36              : 
      37              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_ml_descriptor'
      38              : 
      39              :    PUBLIC :: pao_ml_calc_descriptor
      40              : 
      41              : CONTAINS
      42              : 
      43              : ! **************************************************************************************************
      44              : !> \brief Calculates a descriptor for chemical environment of given atom
      45              : !> \param pao ...
      46              : !> \param particle_set ...
      47              : !> \param qs_kind_set ...
      48              : !> \param cell ...
      49              : !> \param iatom ...
      50              : !> \param descriptor ...
      51              : !> \param descr_grad ...
      52              : !> \param forces ...
      53              : ! **************************************************************************************************
      54          208 :    SUBROUTINE pao_ml_calc_descriptor(pao, particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
      55              :       TYPE(pao_env_type), POINTER                        :: pao
      56              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      57              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      58              :       TYPE(cell_type), POINTER                           :: cell
      59              :       INTEGER, INTENT(IN)                                :: iatom
      60              :       REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL      :: descriptor
      61              :       REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL       :: descr_grad
      62              :       REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
      63              : 
      64              :       CHARACTER(len=*), PARAMETER :: routineN = 'pao_ml_calc_descriptor'
      65              : 
      66              :       INTEGER                                            :: handle
      67              : 
      68          208 :       CALL timeset(routineN, handle)
      69              : 
      70          208 :       CPASSERT(PRESENT(forces) .EQV. PRESENT(descr_grad))
      71              : 
      72          208 :       SELECT CASE (pao%ml_descriptor)
      73              :       CASE (pao_ml_desc_pot)
      74          150 :          CALL calc_descriptor_pot(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
      75              :       CASE (pao_ml_desc_overlap)
      76          118 :          CALL calc_descriptor_overlap(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
      77              :       CASE (pao_ml_desc_r12)
      78          320 :          CALL calc_descriptor_r12(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
      79              :       CASE DEFAULT
      80          208 :          CPABORT("PAO: unknown descriptor")
      81              :       END SELECT
      82              : 
      83          208 :       CALL timestop(handle)
      84          208 :    END SUBROUTINE pao_ml_calc_descriptor
      85              : 
      86              : ! **************************************************************************************************
      87              : !> \brief Calculates a descriptor based on the eigenvalues of V_neighbors
      88              : !> \param particle_set ...
      89              : !> \param qs_kind_set ...
      90              : !> \param cell ...
      91              : !> \param iatom ...
      92              : !> \param descriptor ...
      93              : !> \param descr_grad ...
      94              : !> \param forces ...
      95              : ! **************************************************************************************************
      96           54 :    SUBROUTINE calc_descriptor_pot(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
      97              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      98              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      99              :       TYPE(cell_type), POINTER                           :: cell
     100              :       INTEGER, INTENT(IN)                                :: iatom
     101              :       REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL      :: descriptor
     102              :       REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL       :: descr_grad
     103              :       REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
     104              : 
     105              :       CHARACTER(len=*), PARAMETER :: routineN = 'calc_descriptor_pot'
     106              : 
     107              :       INTEGER                                            :: handle, i, idesc, ikind, jatom, jkind, &
     108              :                                                             k, N, natoms, ndesc
     109              :       REAL(dp)                                           :: beta, w, weight
     110              :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: V_evals
     111           54 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: block_M, block_V, V_evecs
     112           54 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :, :)          :: block_D
     113              :       REAL(dp), DIMENSION(3)                             :: Ra, Rab, Rb
     114              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
     115           54 :       TYPE(pao_descriptor_type), DIMENSION(:), POINTER   :: pao_descriptors
     116              : 
     117           54 :       CALL timeset(routineN, handle)
     118              : 
     119           54 :       CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     120           54 :       CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_descriptors=pao_descriptors)
     121           54 :       N = basis_set%nsgf
     122           54 :       natoms = SIZE(particle_set)
     123           54 :       ndesc = SIZE(pao_descriptors)
     124           54 :       IF (ndesc == 0) CPABORT("No PAO_DESCRIPTOR section found")
     125              : 
     126          432 :       ALLOCATE (block_V(N, N), V_evecs(N, N), V_evals(N))
     127          150 :       IF (PRESENT(descriptor)) ALLOCATE (descriptor(N*ndesc))
     128           90 :       IF (PRESENT(forces)) ALLOCATE (block_D(N, N, 3), block_M(N, N))
     129              : 
     130          152 :       DO idesc = 1, ndesc
     131              : 
     132              :          ! construct matrix V_block from neighboring atoms
     133         3038 :          block_V = 0.0_dp
     134          294 :          DO jatom = 1, natoms
     135          196 :             IF (jatom == iatom) CYCLE
     136          392 :             Ra = particle_set(iatom)%r
     137          392 :             Rb = particle_set(jatom)%r
     138           98 :             Rab = pbc(ra, rb, cell)
     139           98 :             CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
     140           98 :             CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=pao_descriptors)
     141           98 :             IF (SIZE(pao_descriptors) /= ndesc) &
     142            0 :                CPABORT("Not all KINDs have the same number of PAO_DESCRIPTOR sections")
     143           98 :             weight = pao_descriptors(idesc)%weight
     144           98 :             beta = pao_descriptors(idesc)%beta
     145          392 :             CALL pao_calc_gaussian(basis_set, block_V=block_V, Rab=Rab, lpot=0, beta=beta, weight=weight)
     146              :          END DO
     147              : 
     148              :          ! diagonalize block_V
     149         3038 :          V_evecs(:, :) = block_V(:, :)
     150           98 :          CALL diamat_all(V_evecs, V_evals)
     151              : 
     152              :          ! use eigenvalues of V_block as descriptor
     153           98 :          IF (PRESENT(descriptor)) &
     154          528 :             descriptor((idesc - 1)*N + 1:idesc*N) = V_evals(:)
     155              : 
     156              :          ! FORCES ----------------------------------------------------------------------------------
     157          152 :          IF (PRESENT(forces)) THEN
     158           10 :             CPASSERT(PRESENT(descr_grad))
     159          310 :             block_M = 0.0_dp
     160           60 :             DO k = 1, N
     161           50 :                w = descr_grad((idesc - 1)*N + k)
     162         5060 :                block_M(:, :) = block_M(:, :) + w*MATMUL(V_evecs(:, k:k), TRANSPOSE(V_evecs(:, k:k)))
     163              :             END DO
     164           30 :             DO jatom = 1, natoms
     165           20 :                IF (jatom == iatom) CYCLE
     166           40 :                Ra = particle_set(iatom)%r
     167           40 :                Rb = particle_set(jatom)%r
     168           10 :                Rab = pbc(ra, rb, cell)
     169           10 :                CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
     170           10 :                CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=pao_descriptors)
     171           10 :                weight = pao_descriptors(idesc)%weight
     172           10 :                beta = pao_descriptors(idesc)%beta
     173          940 :                block_D = 0.0_dp
     174           10 :                CALL pao_calc_gaussian(basis_set, block_D=block_D, Rab=Rab, lpot=0, beta=beta, weight=weight)
     175           60 :                DO i = 1, 3
     176          930 :                   forces(iatom, i) = forces(iatom, i) - SUM(block_M*block_D(:, :, i))
     177          950 :                   forces(jatom, i) = forces(jatom, i) + SUM(block_M*block_D(:, :, i))
     178              :                END DO
     179              :             END DO
     180              :          END IF
     181              : 
     182              :       END DO
     183              : 
     184           54 :       CALL timestop(handle)
     185          108 :    END SUBROUTINE calc_descriptor_pot
     186              : 
     187              : ! **************************************************************************************************
     188              : !> \brief Calculates a descriptor based on the eigenvalues of local overlap matrix
     189              : !> \param particle_set ...
     190              : !> \param qs_kind_set ...
     191              : !> \param cell ...
     192              : !> \param iatom ...
     193              : !> \param descriptor ...
     194              : !> \param descr_grad ...
     195              : !> \param forces ...
     196              : ! **************************************************************************************************
     197           42 :    SUBROUTINE calc_descriptor_overlap(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
     198              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     199              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     200              :       TYPE(cell_type), POINTER                           :: cell
     201              :       INTEGER, INTENT(IN)                                :: iatom
     202              :       REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL      :: descriptor
     203              :       REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL       :: descr_grad
     204              :       REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
     205              : 
     206              :       CHARACTER(len=*), PARAMETER :: routineN = 'calc_descriptor_overlap'
     207              : 
     208              :       INTEGER                                            :: handle, idesc, ikind, j, jatom, jkind, &
     209              :                                                             k, katom, kkind, N, natoms, ndesc
     210           42 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: neighbor_order
     211              :       REAL(dp) :: beta_sum, deriv, exponent, integral, jbeta, jweight, kbeta, kweight, &
     212              :          normalization, Rij2, Rik2, Rjk2, sbeta, screening_radius, screening_volume, w
     213              :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: S_evals
     214           42 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: block_M, block_S, S_evecs
     215              :       REAL(dp), DIMENSION(3)                             :: Ri, Rij, Rik, Rj, Rjk, Rk
     216           42 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: neighbor_dist
     217           42 :       TYPE(pao_descriptor_type), DIMENSION(:), POINTER   :: ipao_descriptors, jpao_descriptors, &
     218           42 :                                                             kpao_descriptors
     219              : 
     220           42 :       CALL timeset(routineN, handle)
     221              : 
     222           42 :       CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     223           42 :       CALL get_qs_kind(qs_kind_set(ikind), pao_descriptors=ipao_descriptors)
     224              : 
     225           42 :       natoms = SIZE(particle_set)
     226           42 :       ndesc = SIZE(ipao_descriptors)
     227           42 :       IF (ndesc == 0) CPABORT("No PAO_DESCRIPTOR section found")
     228              : 
     229              :       ! determine largest screening radius
     230           42 :       screening_radius = 0.0_dp
     231          118 :       DO idesc = 1, ndesc
     232          118 :          screening_radius = MAX(screening_radius, ipao_descriptors(idesc)%screening_radius)
     233              :       END DO
     234              : 
     235              :       ! estimate maximum number of neighbors within screening
     236           42 :       screening_volume = fourpi/3.0_dp*screening_radius**3
     237           42 :       N = INT(screening_volume/35.0_dp) ! rule of thumb
     238              : 
     239          336 :       ALLOCATE (block_S(N, N), S_evals(N), S_evecs(N, N))
     240          118 :       IF (PRESENT(descriptor)) ALLOCATE (descriptor(N*ndesc))
     241           50 :       IF (PRESENT(forces)) ALLOCATE (block_M(N, N))
     242              : 
     243              :       !find neighbors
     244              :       !TODO: this is a quadratic algorithm, use a neighbor-list instead
     245          210 :       ALLOCATE (neighbor_dist(natoms), neighbor_order(natoms))
     246          168 :       Ri = particle_set(iatom)%r
     247          132 :       DO jatom = 1, natoms
     248          360 :          Rj = particle_set(jatom)%r
     249           90 :          Rij = pbc(Ri, Rj, cell)
     250          402 :          neighbor_dist(jatom) = SQRT(SUM(Rij**2))
     251              :       END DO
     252           42 :       CALL sort(neighbor_dist, natoms, neighbor_order)
     253           42 :       CPASSERT(neighbor_order(1) == iatom) !central atom should be closesd to itself
     254              : 
     255              :       ! check if N was chosen large enough
     256           42 :       IF (natoms > N) THEN
     257            0 :          IF (neighbor_dist(N + 1) < screening_radius) &
     258            0 :             CPABORT("PAO heuristic for descriptor size broke down")
     259              :       END IF
     260              : 
     261          118 :       DO idesc = 1, ndesc
     262           76 :          sbeta = ipao_descriptors(idesc)%screening
     263              : 
     264              :          ! construct matrix S_block from neighboring atoms
     265      1653532 :          block_S = 0.0_dp
     266          234 :          DO j = 1, MIN(natoms, N)
     267          568 :          DO k = 1, MIN(natoms, N)
     268          334 :             jatom = neighbor_order(j)
     269          334 :             katom = neighbor_order(k)
     270              : 
     271              :             ! get weigths and betas
     272          334 :             CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
     273          334 :             CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=jpao_descriptors)
     274          334 :             CALL get_atomic_kind(particle_set(katom)%atomic_kind, kind_number=kkind)
     275          334 :             CALL get_qs_kind(qs_kind_set(kkind), pao_descriptors=kpao_descriptors)
     276          334 :             IF (SIZE(jpao_descriptors) /= ndesc .OR. SIZE(kpao_descriptors) /= ndesc) &
     277            0 :                CPABORT("Not all KINDs have the same number of PAO_DESCRIPTOR sections")
     278          334 :             jweight = jpao_descriptors(idesc)%weight
     279          334 :             jbeta = jpao_descriptors(idesc)%beta
     280          334 :             kweight = kpao_descriptors(idesc)%weight
     281          334 :             kbeta = kpao_descriptors(idesc)%beta
     282          334 :             beta_sum = sbeta + jbeta + kbeta
     283              : 
     284              :             ! get distances
     285         1336 :             Rj = particle_set(jatom)%r
     286         1336 :             Rk = particle_set(katom)%r
     287          334 :             Rij = pbc(Ri, Rj, cell)
     288          334 :             Rik = pbc(Ri, Rk, cell)
     289          334 :             Rjk = pbc(Rj, Rk, cell)
     290         1336 :             Rij2 = SUM(Rij**2)
     291         1336 :             Rik2 = SUM(Rik**2)
     292         1336 :             Rjk2 = SUM(Rjk**2)
     293              : 
     294              :             ! calculate integral over three Gaussians
     295          334 :             exponent = -(sbeta*jbeta*Rij2 + sbeta*kbeta*Rik2 + jbeta*kbeta*Rjk2)/beta_sum
     296          334 :             integral = EXP(exponent)*rootpi/SQRT(beta_sum)
     297          334 :             normalization = SQRT(jbeta*kbeta)/rootpi**2
     298         1160 :             block_S(j, k) = jweight*kweight*normalization*integral
     299              :          END DO
     300              :          END DO
     301              : 
     302              :          ! diagonalize V_block
     303      1653532 :          S_evecs(:, :) = block_S(:, :)
     304           76 :          CALL diamat_all(S_evecs, S_evals)
     305              : 
     306              :          ! use eigenvalues of S_block as descriptor
     307           76 :          IF (PRESENT(descriptor)) &
     308        10360 :             descriptor((idesc - 1)*N + 1:idesc*N) = S_evals(:)
     309              : 
     310              :          ! FORCES ----------------------------------------------------------------------------------
     311          118 :          IF (PRESENT(forces)) THEN
     312            6 :             CPASSERT(PRESENT(descr_grad))
     313       130542 :             block_M = 0.0_dp
     314          888 :             DO k = 1, N
     315          882 :                w = descr_grad((idesc - 1)*N + k)
     316     57699564 :                block_M(:, :) = block_M(:, :) + w*MATMUL(S_evecs(:, k:k), TRANSPOSE(S_evecs(:, k:k)))
     317              :             END DO
     318              : 
     319           20 :             DO j = 1, MIN(natoms, N)
     320           54 :             DO k = 1, MIN(natoms, N)
     321           34 :                jatom = neighbor_order(j)
     322           34 :                katom = neighbor_order(k)
     323              : 
     324              :                ! get weigths and betas
     325           34 :                CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
     326           34 :                CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=jpao_descriptors)
     327           34 :                CALL get_atomic_kind(particle_set(katom)%atomic_kind, kind_number=kkind)
     328           34 :                CALL get_qs_kind(qs_kind_set(kkind), pao_descriptors=kpao_descriptors)
     329           34 :                jweight = jpao_descriptors(idesc)%weight
     330           34 :                jbeta = jpao_descriptors(idesc)%beta
     331           34 :                kweight = kpao_descriptors(idesc)%weight
     332           34 :                kbeta = kpao_descriptors(idesc)%beta
     333           34 :                beta_sum = sbeta + jbeta + kbeta
     334              : 
     335              :                ! get distances
     336          136 :                Rj = particle_set(jatom)%r
     337          136 :                Rk = particle_set(katom)%r
     338           34 :                Rij = pbc(Ri, Rj, cell)
     339           34 :                Rik = pbc(Ri, Rk, cell)
     340           34 :                Rjk = pbc(Rj, Rk, cell)
     341          136 :                Rij2 = SUM(Rij**2)
     342          136 :                Rik2 = SUM(Rik**2)
     343          136 :                Rjk2 = SUM(Rjk**2)
     344              : 
     345              :                ! calculate integral over three Gaussians
     346           34 :                exponent = -(sbeta*jbeta*Rij2 + sbeta*kbeta*Rik2 + jbeta*kbeta*Rjk2)/beta_sum
     347           34 :                integral = EXP(exponent)*rootpi/SQRT(beta_sum)
     348           34 :                normalization = SQRT(jbeta*kbeta)/rootpi**2
     349           34 :                deriv = 2.0_dp/beta_sum*block_M(j, k)
     350           34 :                w = jweight*kweight*normalization*integral*deriv
     351          136 :                forces(iatom, :) = forces(iatom, :) - sbeta*jbeta*Rij*w
     352          136 :                forces(jatom, :) = forces(jatom, :) + sbeta*jbeta*Rij*w
     353          136 :                forces(iatom, :) = forces(iatom, :) - sbeta*kbeta*Rik*w
     354          136 :                forces(katom, :) = forces(katom, :) + sbeta*kbeta*Rik*w
     355          136 :                forces(jatom, :) = forces(jatom, :) - jbeta*kbeta*Rjk*w
     356          218 :                forces(katom, :) = forces(katom, :) + jbeta*kbeta*Rjk*w
     357              :             END DO
     358              :             END DO
     359              :          END IF
     360              :       END DO
     361              : 
     362           42 :       CALL timestop(handle)
     363           84 :    END SUBROUTINE calc_descriptor_overlap
     364              : 
     365              : ! **************************************************************************************************
     366              : !> \brief Calculates a descriptor based on distance between two atoms
     367              : !> \param particle_set ...
     368              : !> \param qs_kind_set ...
     369              : !> \param cell ...
     370              : !> \param iatom ...
     371              : !> \param descriptor ...
     372              : !> \param descr_grad ...
     373              : !> \param forces ...
     374              : ! **************************************************************************************************
     375          112 :    SUBROUTINE calc_descriptor_r12(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
     376              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     377              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     378              :       TYPE(cell_type), POINTER                           :: cell
     379              :       INTEGER, INTENT(IN)                                :: iatom
     380              :       REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL      :: descriptor
     381              :       REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL       :: descr_grad
     382              :       REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
     383              : 
     384              :       REAL(dp), DIMENSION(3)                             :: G, R1, R12, R2
     385              : 
     386          112 :       CPASSERT(SIZE(particle_set) == 2)
     387              : 
     388              :       MARK_USED(qs_kind_set)
     389              :       MARK_USED(iatom)
     390              :       MARK_USED(cell)
     391              : 
     392          448 :       R1 = particle_set(1)%r
     393          448 :       R2 = particle_set(2)%r
     394          112 :       R12 = pbc(R1, R2, cell)
     395              : 
     396          112 :       IF (PRESENT(descriptor)) THEN
     397          104 :          ALLOCATE (descriptor(1))
     398          416 :          descriptor(1) = SQRT(SUM(R12**2))
     399              :       END IF
     400              : 
     401          112 :       IF (PRESENT(forces)) THEN
     402            8 :          CPASSERT(PRESENT(descr_grad))
     403           56 :          G = R12/SQRT(SUM(R12**2))*descr_grad(1)
     404           32 :          forces(1, :) = forces(1, :) + G
     405           32 :          forces(2, :) = forces(2, :) - G
     406              :       END IF
     407          112 :    END SUBROUTINE calc_descriptor_r12
     408              : 
     409          932 : END MODULE pao_ml_descriptor
        

Generated by: LCOV version 2.0-1