Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Util mixed_environment
      10             : !> \author Teodoro Laino [tlaino] - 02.2011
      11             : ! **************************************************************************************************
      12             : MODULE mixed_environment_utils
      13             : 
      14             :    USE cp_result_methods,               ONLY: cp_results_erase,&
      15             :                                               get_results,&
      16             :                                               put_results,&
      17             :                                               test_for_result
      18             :    USE cp_result_types,                 ONLY: cp_result_p_type,&
      19             :                                               cp_result_type
      20             :    USE input_section_types,             ONLY: section_vals_get,&
      21             :                                               section_vals_get_subs_vals,&
      22             :                                               section_vals_type,&
      23             :                                               section_vals_val_get
      24             :    USE kinds,                           ONLY: default_string_length,&
      25             :                                               dp
      26             :    USE mixed_energy_types,              ONLY: mixed_force_type
      27             :    USE particle_list_types,             ONLY: particle_list_type
      28             :    USE virial_types,                    ONLY: virial_p_type,&
      29             :                                               virial_type,&
      30             :                                               zero_virial
      31             : #include "./base/base_uses.f90"
      32             : 
      33             :    IMPLICIT NONE
      34             : 
      35             :    PRIVATE
      36             : 
      37             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mixed_environment_utils'
      38             : 
      39             :    PUBLIC :: mixed_map_forces, &
      40             :              get_subsys_map_index
      41             : 
      42             : CONTAINS
      43             : 
      44             : ! **************************************************************************************************
      45             : !> \brief Maps forces between the different force_eval sections/environments
      46             : !> \param particles_mix ...
      47             : !> \param virial_mix ...
      48             : !> \param results_mix ...
      49             : !> \param global_forces ...
      50             : !> \param virials ...
      51             : !> \param results ...
      52             : !> \param factor ...
      53             : !> \param iforce_eval ...
      54             : !> \param nforce_eval ...
      55             : !> \param map_index ...
      56             : !> \param mapping_section ...
      57             : !> \param overwrite ...
      58             : !> \author Teodoro Laino - University of Zurich [tlaino] - 05.2007
      59             : ! **************************************************************************************************
      60        1116 :    SUBROUTINE mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, &
      61             :                                virials, results, factor, iforce_eval, nforce_eval, map_index, &
      62             :                                mapping_section, overwrite)
      63             : 
      64             :       TYPE(particle_list_type), POINTER                  :: particles_mix
      65             :       TYPE(virial_type), POINTER                         :: virial_mix
      66             :       TYPE(cp_result_type), POINTER                      :: results_mix
      67             :       TYPE(mixed_force_type), DIMENSION(:), POINTER      :: global_forces
      68             :       TYPE(virial_p_type), DIMENSION(:), POINTER         :: virials
      69             :       TYPE(cp_result_p_type), DIMENSION(:), POINTER      :: results
      70             :       REAL(KIND=dp), INTENT(IN)                          :: factor
      71             :       INTEGER, INTENT(IN)                                :: iforce_eval, nforce_eval
      72             :       INTEGER, DIMENSION(:), POINTER                     :: map_index
      73             :       TYPE(section_vals_type), POINTER                   :: mapping_section
      74             :       LOGICAL, INTENT(IN)                                :: overwrite
      75             : 
      76             :       CHARACTER(LEN=default_string_length)               :: description
      77             :       INTEGER                                            :: iparticle, jparticle, natom, nres
      78             :       LOGICAL                                            :: dip_exists
      79             :       REAL(KIND=dp), DIMENSION(3)                        :: dip_mix, dip_tmp
      80             : 
      81             : ! Get Mapping index array
      82             : 
      83        1116 :       natom = SIZE(global_forces(iforce_eval)%forces, 2)
      84        1116 :       CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, map_index)
      85      878786 :       DO iparticle = 1, natom
      86      877670 :          jparticle = map_index(iparticle)
      87      878786 :          IF (overwrite) THEN
      88      405344 :             particles_mix%els(jparticle)%f(:) = factor*global_forces(iforce_eval)%forces(:, iparticle)
      89             :          ELSE
      90             :             particles_mix%els(jparticle)%f(:) = particles_mix%els(jparticle)%f(:) + &
      91     6616016 :                                                 factor*global_forces(iforce_eval)%forces(:, iparticle)
      92             :          END IF
      93             :       END DO
      94             :       ! Mixing Virial
      95        1116 :       IF (virial_mix%pv_availability) THEN
      96         180 :          IF (overwrite) CALL zero_virial(virial_mix, reset=.FALSE.)
      97        4680 :          virial_mix%pv_total = virial_mix%pv_total + factor*virials(iforce_eval)%virial%pv_total
      98        4680 :          virial_mix%pv_kinetic = virial_mix%pv_kinetic + factor*virials(iforce_eval)%virial%pv_kinetic
      99        4680 :          virial_mix%pv_virial = virial_mix%pv_virial + factor*virials(iforce_eval)%virial%pv_virial
     100        4680 :          virial_mix%pv_xc = virial_mix%pv_xc + factor*virials(iforce_eval)%virial%pv_xc
     101        4680 :          virial_mix%pv_fock_4c = virial_mix%pv_fock_4c + factor*virials(iforce_eval)%virial%pv_fock_4c
     102        4680 :          virial_mix%pv_constraint = virial_mix%pv_constraint + factor*virials(iforce_eval)%virial%pv_constraint
     103             :       END IF
     104             :       ! Deallocate map_index array
     105        1116 :       IF (ASSOCIATED(map_index)) THEN
     106        1116 :          DEALLOCATE (map_index)
     107             :       END IF
     108             : 
     109             :       ! Collect Requested Results info
     110        1116 :       description = '[DIPOLE]'
     111        1116 :       IF (overwrite) CALL cp_results_erase(results_mix)
     112             : 
     113        1116 :       dip_exists = test_for_result(results=results(iforce_eval)%results, description=description)
     114        1116 :       IF (dip_exists) THEN
     115         348 :          CALL get_results(results=results_mix, description=description, n_rep=nres)
     116         348 :          CPASSERT(nres <= 1)
     117         348 :          dip_mix = 0.0_dp
     118         348 :          IF (nres == 1) CALL get_results(results=results_mix, description=description, values=dip_mix)
     119         348 :          CALL get_results(results=results(iforce_eval)%results, description=description, n_rep=nres)
     120             :          CALL get_results(results=results(iforce_eval)%results, description=description, &
     121         348 :                           values=dip_tmp, nval=nres)
     122        1392 :          dip_mix = dip_mix + factor*dip_tmp
     123         348 :          CALL cp_results_erase(results=results_mix, description=description)
     124         348 :          CALL put_results(results=results_mix, description=description, values=dip_mix)
     125             :       END IF
     126             : 
     127        1116 :    END SUBROUTINE mixed_map_forces
     128             : 
     129             : ! **************************************************************************************************
     130             : !> \brief performs mapping of the subsystems of different force_eval
     131             : !> \param mapping_section ...
     132             : !> \param natom ...
     133             : !> \param iforce_eval ...
     134             : !> \param nforce_eval ...
     135             : !> \param map_index ...
     136             : !> \param force_eval_embed ...
     137             : !> \author Teodoro Laino - University of Zurich [tlaino] - 05.2007
     138             : ! **************************************************************************************************
     139        1898 :    SUBROUTINE get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, map_index, &
     140             :                                    force_eval_embed)
     141             : 
     142             :       TYPE(section_vals_type), POINTER                   :: mapping_section
     143             :       INTEGER, INTENT(IN)                                :: natom, iforce_eval, nforce_eval
     144             :       INTEGER, DIMENSION(:), POINTER                     :: map_index
     145             :       LOGICAL, OPTIONAL                                  :: force_eval_embed
     146             : 
     147             :       INTEGER                                            :: i, iatom, ival, j, jval, k, n_rep, &
     148             :                                                             n_rep_loc, n_rep_map, n_rep_sys, tmp
     149        1898 :       INTEGER, DIMENSION(:), POINTER                     :: index_glo, index_loc, list
     150             :       LOGICAL                                            :: check, explicit
     151             :       TYPE(section_vals_type), POINTER                   :: fragments_loc, fragments_sys, &
     152             :                                                             map_force_ev, map_full_sys
     153             : 
     154           0 :       CPASSERT(.NOT. ASSOCIATED(map_index))
     155        5694 :       ALLOCATE (map_index(natom))
     156        1898 :       CALL section_vals_get(mapping_section, explicit=explicit)
     157        1898 :       IF (.NOT. explicit) THEN
     158             :          ! Standard Mapping.. subsys are assumed to have the same structure
     159      284522 :          DO i = 1, natom
     160      284522 :             map_index(i) = i
     161             :          END DO
     162             :       ELSE
     163             :          ! Mapping systems with different structures
     164         784 :          IF (.NOT. PRESENT(force_eval_embed)) THEN
     165         684 :             map_full_sys => section_vals_get_subs_vals(mapping_section, "FORCE_EVAL_MIXED")
     166             :          ELSE
     167         100 :             map_full_sys => section_vals_get_subs_vals(mapping_section, "FORCE_EVAL_EMBED")
     168             :          END IF
     169         784 :          map_force_ev => section_vals_get_subs_vals(mapping_section, "FORCE_EVAL")
     170         784 :          CALL section_vals_get(map_full_sys, explicit=explicit)
     171         784 :          CPASSERT(explicit)
     172         784 :          CALL section_vals_get(map_force_ev, explicit=explicit, n_repetition=n_rep)
     173         784 :          CPASSERT(explicit)
     174         784 :          CPASSERT(n_rep == nforce_eval)
     175        1476 :          DO i = 1, n_rep
     176        1476 :             CALL section_vals_val_get(map_force_ev, "_SECTION_PARAMETERS_", i_rep_section=i, i_val=ival)
     177        1476 :             IF (ival == iforce_eval) EXIT
     178             :          END DO
     179         784 :          CPASSERT(i <= nforce_eval)
     180             :          MARK_USED(nforce_eval)
     181         784 :          fragments_sys => section_vals_get_subs_vals(map_full_sys, "FRAGMENT")
     182         784 :          fragments_loc => section_vals_get_subs_vals(map_force_ev, "FRAGMENT", i_rep_section=i)
     183             :          !Perform few check on the structure of the input mapping section. as provided by the user
     184         784 :          CALL section_vals_get(fragments_loc, n_repetition=n_rep_loc)
     185         784 :          CALL section_vals_get(fragments_sys, explicit=explicit, n_repetition=n_rep_sys)
     186         784 :          CPASSERT(explicit)
     187         784 :          CPASSERT(n_rep_sys >= n_rep_loc)
     188         784 :          IF (n_rep_loc == 0) THEN
     189         126 :             NULLIFY (list)
     190             :             ! We expect an easier syntax in this case..
     191         126 :             CALL section_vals_val_get(map_force_ev, "DEFINE_FRAGMENTS", i_rep_section=i, n_rep_val=n_rep_map)
     192         126 :             check = (n_rep_map /= 0)
     193         126 :             CPASSERT(check)
     194         126 :             CALL section_vals_val_get(map_force_ev, "DEFINE_FRAGMENTS", i_rep_section=i, i_vals=list)
     195         126 :             CPASSERT(SIZE(list) > 0)
     196         126 :             iatom = 0
     197         630 :             DO i = 1, SIZE(list)
     198         504 :                jval = list(i)
     199        1512 :                DO j = 1, n_rep_sys
     200        1512 :                   CALL section_vals_val_get(fragments_sys, "_SECTION_PARAMETERS_", i_rep_section=j, i_val=tmp)
     201        1512 :                   IF (tmp == jval) EXIT
     202             :                END DO
     203         504 :                CALL section_vals_val_get(fragments_sys, "_DEFAULT_KEYWORD_", i_rep_section=j, i_vals=index_glo)
     204      232533 :                DO k = 0, index_glo(2) - index_glo(1)
     205      231903 :                   iatom = iatom + 1
     206      231903 :                   CPASSERT(iatom <= natom)
     207      232407 :                   map_index(iatom) = index_glo(1) + k
     208             :                END DO
     209             :             END DO
     210         126 :             check = (iatom == natom)
     211         126 :             CPASSERT(check)
     212             :          ELSE
     213             :             ! General syntax..
     214             :             !Loop over the fragment of the force_eval
     215        2612 :             DO i = 1, n_rep_loc
     216        1954 :                CALL section_vals_val_get(fragments_loc, "_SECTION_PARAMETERS_", i_rep_section=i, i_val=ival)
     217        1954 :                CALL section_vals_val_get(fragments_loc, "MAP", i_rep_section=i, i_val=jval)
     218             :                ! Index corresponding to the mixed_force_eval fragment
     219        5520 :                DO j = 1, n_rep_sys
     220        5520 :                   CALL section_vals_val_get(fragments_sys, "_SECTION_PARAMETERS_", i_rep_section=j, i_val=tmp)
     221        5520 :                   IF (tmp == jval) EXIT
     222             :                END DO
     223        1954 :                CPASSERT(j <= n_rep_sys)
     224        1954 :                CALL section_vals_val_get(fragments_loc, "_DEFAULT_KEYWORD_", i_rep_section=i, i_vals=index_loc)
     225        1954 :                CALL section_vals_val_get(fragments_sys, "_DEFAULT_KEYWORD_", i_rep_section=j, i_vals=index_glo)
     226        1954 :                check = ((index_loc(2) - index_loc(1)) == (index_glo(2) - index_glo(1)))
     227        1954 :                CPASSERT(check)
     228             :                ! Now let's build the real mapping
     229      806342 :                DO k = 0, index_loc(2) - index_loc(1)
     230      803730 :                   map_index(index_loc(1) + k) = index_glo(1) + k
     231             :                END DO
     232             :             END DO
     233             :          END IF
     234             :       END IF
     235             : 
     236        3796 :    END SUBROUTINE get_subsys_map_index
     237             : 
     238             : END MODULE mixed_environment_utils