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             : !> \par History
      10             : !>      HAF (16-Apr-2025) : Import into CP2K
      11             : !> \author HAF and yury-lysogorskiy and ralf-drautz
      12             : ! **************************************************************************************************
      13             : 
      14             : MODULE manybody_ace
      15             : 
      16             :    USE ISO_C_BINDING,                   ONLY: C_ASSOCIATED
      17             :    USE ace_nlist,                       ONLY: ace_interface
      18             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      19             :                                               get_atomic_kind
      20             :    USE bibliography,                    ONLY: Bochkarev2024,&
      21             :                                               Drautz2019,&
      22             :                                               Lysogorskiy2021,&
      23             :                                               cite_reference
      24             :    USE cell_types,                      ONLY: cell_type
      25             :    USE fist_nonbond_env_types,          ONLY: ace_data_type,&
      26             :                                               fist_nonbond_env_get,&
      27             :                                               fist_nonbond_env_set,&
      28             :                                               fist_nonbond_env_type
      29             :    USE kinds,                           ONLY: dp
      30             :    USE pair_potential_types,            ONLY: ace_type,&
      31             :                                               pair_potential_pp_type,&
      32             :                                               pair_potential_single_type
      33             :    USE particle_types,                  ONLY: particle_type
      34             :    USE physcon,                         ONLY: angstrom,&
      35             :                                               evolt
      36             : #include "./base/base_uses.f90"
      37             : 
      38             :    IMPLICIT NONE
      39             : 
      40             :    PRIVATE
      41             :    PUBLIC ace_energy_store_force_virial, ace_add_force_virial
      42             : 
      43             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_ace'
      44             : 
      45             : CONTAINS
      46             : 
      47             : ! **************************************************************************************************
      48             : !> \brief ...
      49             : !> \param particle_set ...
      50             : !> \param atomic_kind_set ...
      51             : !> \param potparm ...
      52             : !> \param ace_data ...
      53             : ! **************************************************************************************************
      54           6 :    SUBROUTINE init_ace_data(particle_set, atomic_kind_set, potparm, &
      55             :                             ace_data)
      56             : 
      57             :       TYPE(particle_type), POINTER                       :: particle_set(:)
      58             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
      59             :       TYPE(pair_potential_pp_type), POINTER              :: potparm
      60             :       TYPE(ace_data_type), POINTER                       :: ace_data
      61             : 
      62             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'init_ace_data'
      63             : 
      64             :       CHARACTER(2)                                       :: element_symbol
      65             :       INTEGER                                            :: ace_natom, handle, i, iat, iat_use, &
      66             :                                                             ikind, jkind, lkind, n_atoms
      67           6 :       INTEGER, ALLOCATABLE                               :: use_atom_type(:)
      68           6 :       INTEGER, DIMENSION(:), POINTER                     :: ak_alist
      69           6 :       LOGICAL, ALLOCATABLE                               :: use_atom(:)
      70             :       TYPE(pair_potential_single_type), POINTER          :: pot
      71             : 
      72           6 :       CALL timeset(routineN, handle)
      73             : 
      74             :       ! init ace_data
      75           6 :       IF (.NOT. ASSOCIATED(ace_data)) THEN
      76           0 :          ALLOCATE (ace_data)
      77             :       END IF
      78             : 
      79           6 :       n_atoms = SIZE(particle_set)
      80          18 :       ALLOCATE (use_atom(n_atoms))
      81          12 :       ALLOCATE (use_atom_type(n_atoms))
      82         996 :       use_atom = .FALSE.
      83         996 :       use_atom_type = 0
      84             : 
      85          18 :       DO ikind = 1, SIZE(atomic_kind_set)
      86          12 :          pot => potparm%pot(ikind, ikind)%pot
      87          30 :          DO i = 1, SIZE(pot%type)
      88          12 :             IF (pot%type(i) /= ace_type) CYCLE
      89             :             CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), &
      90             :                                  element_symbol=element_symbol, &
      91          12 :                                  natom=lkind, atom_list=ak_alist)
      92          12 :             IF (lkind < 1) CYCLE
      93          12 :             ace_data%model = pot%set(i)%ace%model
      94          12 :             jkind = 0
      95          18 :             DO iat = 1, SIZE(ace_data%model%symbolc)
      96          18 :                IF (element_symbol == ace_data%model%symbolc(iat)) THEN
      97             :                   jkind = iat
      98             :                   EXIT
      99             :                END IF
     100             :             END DO
     101          12 :             CPASSERT(jkind > 0)
     102        1026 :             DO iat = 1, lkind
     103         990 :                use_atom_type(ak_alist(iat)) = jkind
     104        1002 :                use_atom(ak_alist(iat)) = .TRUE.
     105             :             END DO
     106             :          END DO ! i
     107             :       END DO ! ikind
     108           6 :       CPASSERT(C_ASSOCIATED(ace_data%model%c_ptr))
     109             : 
     110         996 :       ace_natom = COUNT(use_atom)
     111             : 
     112           6 :       IF (.NOT. ALLOCATED(ace_data%uctype)) THEN
     113          18 :          ALLOCATE (ace_data%uctype(1:ace_natom))
     114             :       END IF
     115             : 
     116             :       iat_use = 0
     117         996 :       DO iat = 1, n_atoms
     118         990 :          IF (.NOT. use_atom(iat)) CYCLE
     119         990 :          iat_use = iat_use + 1
     120         996 :          ace_data%uctype(iat_use) = use_atom_type(iat)
     121             :       END DO
     122             : 
     123           6 :       IF (iat_use > 0) THEN
     124           6 :          CALL cite_reference(Drautz2019)
     125           6 :          CALL cite_reference(Lysogorskiy2021)
     126           6 :          CALL cite_reference(Bochkarev2024)
     127             :       END IF
     128             : 
     129           6 :       IF (.NOT. ALLOCATED(ace_data%force)) THEN
     130          18 :          ALLOCATE (ace_data%force(3, ace_natom))
     131          18 :          ALLOCATE (ace_data%use_indices(ace_natom))
     132          12 :          ALLOCATE (ace_data%inverse_index_map(n_atoms))
     133             :       END IF
     134           6 :       CPASSERT(SIZE(ace_data%force, 2) == ace_natom)
     135             : 
     136           6 :       iat_use = 0
     137         996 :       ace_data%inverse_index_map(:) = 0
     138         996 :       DO iat = 1, n_atoms
     139         996 :          IF (use_atom(iat)) THEN
     140         990 :             iat_use = iat_use + 1
     141         990 :             ace_data%use_indices(iat_use) = iat
     142         990 :             ace_data%inverse_index_map(iat) = iat_use
     143             :          END IF
     144             :       END DO
     145           6 :       ace_data%natom = ace_natom
     146           6 :       DEALLOCATE (use_atom, use_atom_type)
     147             : 
     148           6 :       CALL timestop(handle)
     149             : 
     150           6 :    END SUBROUTINE init_ace_data
     151             : 
     152             : ! **************************************************************************************************
     153             : !> \brief ...
     154             : !>     > \param particle_set ...
     155             : !> \param particle_set ...
     156             : !> \param cell ...
     157             : !> \param atomic_kind_set ...
     158             : !> \param potparm ...
     159             : !> \param fist_nonbond_env ...
     160             : !> \param pot_ace ...
     161             : ! **************************************************************************************************
     162         412 :    SUBROUTINE ace_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, &
     163             :                                             fist_nonbond_env, pot_ace)
     164             : 
     165             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     166             :       TYPE(cell_type), POINTER                           :: cell
     167             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
     168             :       TYPE(pair_potential_pp_type), POINTER              :: potparm
     169             :       TYPE(fist_nonbond_env_type), POINTER               :: fist_nonbond_env
     170             :       REAL(kind=dp), INTENT(OUT)                         :: pot_ace
     171             : 
     172             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'ace_energy_store_force_virial'
     173             : 
     174             :       INTEGER                                            :: handle
     175             :       REAL(kind=dp)                                      :: ace_virial(1:6)
     176             :       TYPE(ace_data_type), POINTER                       :: ace_data
     177             : 
     178         206 :       CALL timeset(routineN, handle)
     179             : 
     180             :       ! get ace_data to save force, virial info
     181         206 :       CALL fist_nonbond_env_get(fist_nonbond_env, ace_data=ace_data)
     182         206 :       IF (.NOT. ASSOCIATED(ace_data)) THEN
     183          84 :          ALLOCATE (ace_data)
     184             :          !initialize ace_data:
     185           6 :          CALL init_ace_data(particle_set, atomic_kind_set, potparm, ace_data)
     186           6 :          CALL fist_nonbond_env_set(fist_nonbond_env, ace_data=ace_data)
     187             :       END IF
     188             : 
     189             :       CALL ace_interface(ace_data%natom, ace_data%uctype, &
     190             :                          pot_ace, ace_data%force, ace_virial, &
     191         206 :                          fist_nonbond_env, cell, ace_data)
     192             : 
     193             :       ! convert units
     194         206 :       pot_ace = pot_ace/evolt
     195      157766 :       ace_data%force = ace_data%force/(evolt/angstrom)
     196        1442 :       ace_virial = ace_virial/evolt
     197             : 
     198             :       ! minus sign due to CP2K conventions
     199         206 :       ace_data%virial(1, 1) = -ace_virial(1)
     200         206 :       ace_data%virial(2, 2) = -ace_virial(2)
     201         206 :       ace_data%virial(3, 3) = -ace_virial(3)
     202         206 :       ace_data%virial(1, 2) = -ace_virial(4)
     203         206 :       ace_data%virial(2, 1) = -ace_virial(4)
     204         206 :       ace_data%virial(1, 3) = -ace_virial(5)
     205         206 :       ace_data%virial(3, 1) = -ace_virial(5)
     206         206 :       ace_data%virial(2, 3) = -ace_virial(6)
     207         206 :       ace_data%virial(3, 2) = -ace_virial(6)
     208             : 
     209         206 :       CALL timestop(handle)
     210         206 :    END SUBROUTINE ace_energy_store_force_virial
     211             : 
     212             : ! **************************************************************************************************
     213             : !> \brief ...
     214             : !> \param fist_nonbond_env ...
     215             : !> \param force ...
     216             : !> \param pv_nonbond ...
     217             : !> \param use_virial ...
     218             : ! **************************************************************************************************
     219         206 :    SUBROUTINE ace_add_force_virial(fist_nonbond_env, force, pv_nonbond, use_virial)
     220             :       TYPE(fist_nonbond_env_type), POINTER               :: fist_nonbond_env
     221             :       REAL(KIND=dp), INTENT(INOUT)                       :: force(:, :), pv_nonbond(3, 3)
     222             :       LOGICAL, OPTIONAL                                  :: use_virial
     223             : 
     224             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'ace_add_force_virial'
     225             : 
     226             :       INTEGER                                            :: handle, iat, iat_use
     227             :       TYPE(ace_data_type), POINTER                       :: ace_data
     228             : 
     229         206 :       CALL timeset(routineN, handle)
     230             : 
     231         206 :       CALL fist_nonbond_env_get(fist_nonbond_env, ace_data=ace_data)
     232             : 
     233         206 :       IF (.NOT. ASSOCIATED(ace_data)) RETURN
     234             : 
     235       39596 :       DO iat_use = 1, SIZE(ace_data%use_indices)
     236       39390 :          iat = ace_data%use_indices(iat_use)
     237       39390 :          CPASSERT(iat >= 1 .AND. iat <= SIZE(force, 2))
     238      157766 :          force(1:3, iat) = force(1:3, iat) + ace_data%force(1:3, iat_use)
     239             :       END DO
     240             : 
     241         206 :       IF (use_virial) THEN
     242          26 :          pv_nonbond = pv_nonbond + ace_data%virial
     243             :       END IF
     244             : 
     245         206 :       CALL timestop(handle)
     246             :    END SUBROUTINE ace_add_force_virial
     247             : 
     248             : END MODULE manybody_ace