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 Initialize a small environment for a particular calculation
      10             : !> \par History
      11             : !>      5.2004 created [fawzi]
      12             : !>      9.2007 cleaned [tlaino] - University of Zurich
      13             : !> \author Teodoro Laino
      14             : ! **************************************************************************************************
      15             : MODULE cp_subsys_methods
      16             :    USE atomic_kind_list_types,          ONLY: atomic_kind_list_create,&
      17             :                                               atomic_kind_list_release,&
      18             :                                               atomic_kind_list_type
      19             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      20             :    USE atprop_types,                    ONLY: atprop_create
      21             :    USE cell_types,                      ONLY: cell_retain,&
      22             :                                               cell_type
      23             :    USE colvar_methods,                  ONLY: colvar_read
      24             :    USE cp_result_types,                 ONLY: cp_result_create
      25             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      26             :                                               cp_subsys_set,&
      27             :                                               cp_subsys_type
      28             :    USE exclusion_types,                 ONLY: exclusion_type
      29             :    USE input_constants,                 ONLY: do_conn_off,&
      30             :                                               do_stress_analytical,&
      31             :                                               do_stress_diagonal_anal,&
      32             :                                               do_stress_diagonal_numer,&
      33             :                                               do_stress_none,&
      34             :                                               do_stress_numerical
      35             :    USE input_section_types,             ONLY: section_vals_get,&
      36             :                                               section_vals_get_subs_vals,&
      37             :                                               section_vals_type,&
      38             :                                               section_vals_val_get
      39             :    USE kinds,                           ONLY: default_string_length,&
      40             :                                               dp
      41             :    USE message_passing,                 ONLY: mp_para_env_type
      42             :    USE molecule_kind_list_types,        ONLY: molecule_kind_list_create,&
      43             :                                               molecule_kind_list_release,&
      44             :                                               molecule_kind_list_type
      45             :    USE molecule_kind_types,             ONLY: molecule_kind_type
      46             :    USE molecule_list_types,             ONLY: molecule_list_create,&
      47             :                                               molecule_list_release,&
      48             :                                               molecule_list_type
      49             :    USE molecule_types,                  ONLY: molecule_type
      50             :    USE particle_list_types,             ONLY: particle_list_create,&
      51             :                                               particle_list_release,&
      52             :                                               particle_list_type
      53             :    USE particle_types,                  ONLY: particle_type
      54             :    USE qmmm_types_low,                  ONLY: qmmm_env_mm_type
      55             :    USE string_table,                    ONLY: id2str,&
      56             :                                               s2s,&
      57             :                                               str2id
      58             :    USE topology,                        ONLY: connectivity_control,&
      59             :                                               topology_control
      60             :    USE topology_connectivity_util,      ONLY: topology_connectivity_pack
      61             :    USE topology_coordinate_util,        ONLY: topology_coordinate_pack
      62             :    USE topology_types,                  ONLY: deallocate_topology,&
      63             :                                               init_topology,&
      64             :                                               topology_parameters_type
      65             :    USE topology_util,                   ONLY: check_subsys_element
      66             :    USE virial_types,                    ONLY: virial_set
      67             : #include "./base/base_uses.f90"
      68             : 
      69             :    IMPLICIT NONE
      70             :    PRIVATE
      71             : 
      72             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
      73             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_subsys_methods'
      74             : 
      75             :    PUBLIC :: create_small_subsys, cp_subsys_create
      76             : 
      77             : CONTAINS
      78             : 
      79             : ! **************************************************************************************************
      80             : !> \brief Creates allocates and fills subsys from given input.
      81             : !> \param subsys ...
      82             : !> \param para_env ...
      83             : !> \param root_section ...
      84             : !> \param force_env_section ...
      85             : !> \param subsys_section ...
      86             : !> \param use_motion_section ...
      87             : !> \param qmmm ...
      88             : !> \param qmmm_env ...
      89             : !> \param exclusions ...
      90             : !> \param elkind ...
      91             : !> \author Ole Schuett
      92             : ! **************************************************************************************************
      93       28791 :    SUBROUTINE cp_subsys_create(subsys, para_env, &
      94             :                                root_section, force_env_section, subsys_section, &
      95             :                                use_motion_section, qmmm, qmmm_env, exclusions, elkind)
      96             :       TYPE(cp_subsys_type), POINTER                      :: subsys
      97             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      98             :       TYPE(section_vals_type), POINTER                   :: root_section
      99             :       TYPE(section_vals_type), OPTIONAL, POINTER         :: force_env_section, subsys_section
     100             :       LOGICAL, INTENT(IN), OPTIONAL                      :: use_motion_section
     101             :       LOGICAL, OPTIONAL                                  :: qmmm
     102             :       TYPE(qmmm_env_mm_type), OPTIONAL, POINTER          :: qmmm_env
     103             :       TYPE(exclusion_type), DIMENSION(:), OPTIONAL, &
     104             :          POINTER                                         :: exclusions
     105             :       LOGICAL, INTENT(IN), OPTIONAL                      :: elkind
     106             : 
     107             :       INTEGER                                            :: stress_tensor
     108        9597 :       INTEGER, DIMENSION(:), POINTER                     :: seed_vals
     109             :       LOGICAL                                            :: atomic_energy, my_use_motion_section, &
     110             :                                                             pv_availability, pv_diagonal, &
     111             :                                                             pv_numerical
     112             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     113        9597 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     114             :       TYPE(molecule_kind_list_type), POINTER             :: mol_kinds
     115        9597 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     116             :       TYPE(molecule_list_type), POINTER                  :: mols
     117        9597 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     118             :       TYPE(particle_list_type), POINTER                  :: particles
     119        9597 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     120             :       TYPE(section_vals_type), POINTER                   :: colvar_section, my_force_env_section, &
     121             :                                                             my_subsys_section
     122             : 
     123           0 :       CPASSERT(.NOT. ASSOCIATED(subsys))
     124       95970 :       ALLOCATE (subsys)
     125             : 
     126        9597 :       CALL para_env%retain()
     127        9597 :       subsys%para_env => para_env
     128             : 
     129        9597 :       my_use_motion_section = .FALSE.
     130        9597 :       IF (PRESENT(use_motion_section)) &
     131        9595 :          my_use_motion_section = use_motion_section
     132             : 
     133        9597 :       my_force_env_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
     134        9597 :       IF (PRESENT(force_env_section)) &
     135        9595 :          my_force_env_section => force_env_section
     136             : 
     137        9597 :       my_subsys_section => section_vals_get_subs_vals(my_force_env_section, "SUBSYS")
     138        9597 :       IF (PRESENT(subsys_section)) &
     139        9441 :          my_subsys_section => subsys_section
     140             : 
     141        9597 :       CALL section_vals_val_get(my_subsys_section, "SEED", i_vals=seed_vals)
     142        9597 :       IF (SIZE(seed_vals) == 1) THEN
     143       86337 :          subsys%seed(:, :) = REAL(seed_vals(1), KIND=dp)
     144           4 :       ELSE IF (SIZE(seed_vals) == 6) THEN
     145          60 :          subsys%seed(1:3, 1:2) = RESHAPE(REAL(seed_vals(:), KIND=dp), (/3, 2/))
     146             :       ELSE
     147           0 :          CPABORT("Supply exactly 1 or 6 arguments for SEED in &SUBSYS only!")
     148             :       END IF
     149             : 
     150        9597 :       colvar_section => section_vals_get_subs_vals(my_subsys_section, "COLVAR")
     151             : 
     152        9597 :       CALL cp_subsys_read_colvar(subsys, colvar_section)
     153             : 
     154             :       !   *** Read the particle coordinates and allocate the atomic kind, ***
     155             :       !   *** the molecule kind, and the molecule data structures         ***
     156             :       CALL topology_control(atomic_kind_set, particle_set, molecule_kind_set, molecule_set, &
     157             :                             subsys%colvar_p, subsys%gci, root_section, para_env, &
     158             :                             force_env_section=my_force_env_section, &
     159             :                             subsys_section=my_subsys_section, use_motion_section=my_use_motion_section, &
     160        9597 :                             qmmm=qmmm, qmmm_env=qmmm_env, exclusions=exclusions, elkind=elkind)
     161             : 
     162        9597 :       CALL particle_list_create(particles, els_ptr=particle_set)
     163        9597 :       CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
     164        9597 :       CALL molecule_list_create(mols, els_ptr=molecule_set)
     165        9597 :       CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
     166             : 
     167             :       CALL cp_subsys_set(subsys, particles=particles, atomic_kinds=atomic_kinds, &
     168        9597 :                          molecules=mols, molecule_kinds=mol_kinds)
     169             : 
     170        9597 :       CALL particle_list_release(particles)
     171        9597 :       CALL atomic_kind_list_release(atomic_kinds)
     172        9597 :       CALL molecule_list_release(mols)
     173        9597 :       CALL molecule_kind_list_release(mol_kinds)
     174             : 
     175             :       ! Should we compute the virial?
     176        9597 :       CALL section_vals_val_get(my_force_env_section, "STRESS_TENSOR", i_val=stress_tensor)
     177        8809 :       SELECT CASE (stress_tensor)
     178             :       CASE (do_stress_none)
     179        8809 :          pv_availability = .FALSE.
     180        8809 :          pv_numerical = .FALSE.
     181        8809 :          pv_diagonal = .FALSE.
     182             :       CASE (do_stress_analytical)
     183         776 :          pv_availability = .TRUE.
     184         776 :          pv_numerical = .FALSE.
     185         776 :          pv_diagonal = .FALSE.
     186             :       CASE (do_stress_numerical)
     187           2 :          pv_availability = .TRUE.
     188           2 :          pv_numerical = .TRUE.
     189           2 :          pv_diagonal = .FALSE.
     190             :       CASE (do_stress_diagonal_anal)
     191           0 :          pv_availability = .TRUE.
     192           0 :          pv_numerical = .FALSE.
     193           0 :          pv_diagonal = .TRUE.
     194             :       CASE (do_stress_diagonal_numer)
     195          10 :          pv_availability = .TRUE.
     196          10 :          pv_numerical = .TRUE.
     197        9597 :          pv_diagonal = .TRUE.
     198             :       END SELECT
     199             : 
     200     2380056 :       ALLOCATE (subsys%virial)
     201             :       CALL virial_set(virial=subsys%virial, &
     202             :                       pv_availability=pv_availability, &
     203             :                       pv_numer=pv_numerical, &
     204        9597 :                       pv_diagonal=pv_diagonal)
     205             : 
     206             :       ! Should we compute atomic properties?
     207        9597 :       CALL atprop_create(subsys%atprop)
     208        9597 :       CALL section_vals_val_get(my_force_env_section, "PROPERTIES%ATOMIC%ENERGY", l_val=atomic_energy)
     209        9597 :       subsys%atprop%energy = atomic_energy
     210             : 
     211        9597 :       CALL cp_result_create(subsys%results)
     212        9597 :    END SUBROUTINE cp_subsys_create
     213             : 
     214             : ! **************************************************************************************************
     215             : !> \brief reads the colvar section of the colvar
     216             : !> \param subsys ...
     217             : !> \param colvar_section ...
     218             : !> \par History
     219             : !>      2006.01 Joost VandeVondele
     220             : ! **************************************************************************************************
     221        9597 :    SUBROUTINE cp_subsys_read_colvar(subsys, colvar_section)
     222             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     223             :       TYPE(section_vals_type), POINTER                   :: colvar_section
     224             : 
     225             :       INTEGER                                            :: ig, ncol
     226             : 
     227        9597 :       CALL section_vals_get(colvar_section, n_repetition=ncol)
     228       19928 :       ALLOCATE (subsys%colvar_p(ncol))
     229       10059 :       DO ig = 1, ncol
     230         462 :          NULLIFY (subsys%colvar_p(ig)%colvar)
     231       10059 :          CALL colvar_read(subsys%colvar_p(ig)%colvar, ig, colvar_section, subsys%para_env)
     232             :       END DO
     233        9597 :    END SUBROUTINE cp_subsys_read_colvar
     234             : 
     235             : ! **************************************************************************************************
     236             : !> \brief updates the molecule information of the given subsys
     237             : !> \param small_subsys the subsys to create
     238             : !> \param big_subsys the superset of small_subsys
     239             : !> \param small_cell ...
     240             : !> \param small_para_env the parallel environment for the new (small)
     241             : !>        subsys
     242             : !> \param sub_atom_index indexes of the atoms that should be in small_subsys
     243             : !> \param sub_atom_kind_name ...
     244             : !> \param para_env ...
     245             : !> \param force_env_section ...
     246             : !> \param subsys_section ...
     247             : !> \param ignore_outside_box ...
     248             : !> \par History
     249             : !>      05.2004 created [fawzi]
     250             : !> \author Fawzi Mohamed, Teodoro Laino
     251             : !> \note
     252             : !>      not really ready to be used with different para_envs for the small
     253             : !>      and big part
     254             : ! **************************************************************************************************
     255         538 :    SUBROUTINE create_small_subsys(small_subsys, big_subsys, small_cell, &
     256         538 :                                   small_para_env, sub_atom_index, sub_atom_kind_name, &
     257             :                                   para_env, force_env_section, subsys_section, ignore_outside_box)
     258             : 
     259             :       TYPE(cp_subsys_type), POINTER                      :: small_subsys, big_subsys
     260             :       TYPE(cell_type), POINTER                           :: small_cell
     261             :       TYPE(mp_para_env_type), POINTER                    :: small_para_env
     262             :       INTEGER, DIMENSION(:), INTENT(in)                  :: sub_atom_index
     263             :       CHARACTER(len=default_string_length), &
     264             :          DIMENSION(:), INTENT(in)                        :: sub_atom_kind_name
     265             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     266             :       TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
     267             :       LOGICAL, INTENT(in), OPTIONAL                      :: ignore_outside_box
     268             : 
     269             :       CHARACTER(len=default_string_length)               :: my_element, strtmp1
     270             :       INTEGER                                            :: iat, id_, nat
     271             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     272         538 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     273             :       TYPE(molecule_kind_list_type), POINTER             :: mol_kinds
     274         538 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     275             :       TYPE(molecule_list_type), POINTER                  :: mols
     276         538 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     277             :       TYPE(particle_list_type), POINTER                  :: particles
     278         538 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     279             :       TYPE(topology_parameters_type)                     :: topology
     280             : 
     281         538 :       NULLIFY (mol_kinds, mols, particles, atomic_kinds, atomic_kind_set, particle_set, &
     282         538 :                molecule_kind_set, molecule_set, particles, atomic_kinds)
     283             : 
     284           0 :       CPASSERT(.NOT. ASSOCIATED(small_subsys))
     285         538 :       CPASSERT(ASSOCIATED(big_subsys))
     286         538 :       IF (big_subsys%para_env /= small_para_env) &
     287           0 :          CPABORT("big_subsys%para_env==small_para_env")
     288             : 
     289             :       !-----------------------------------------------------------------------------
     290             :       !-----------------------------------------------------------------------------
     291             :       ! 1. Initialize the topology structure type
     292             :       !-----------------------------------------------------------------------------
     293         538 :       CALL init_topology(topology)
     294             : 
     295             :       !-----------------------------------------------------------------------------
     296             :       !-----------------------------------------------------------------------------
     297             :       ! 2. Get the cell info
     298             :       !-----------------------------------------------------------------------------
     299         538 :       topology%cell => small_cell
     300         538 :       CALL cell_retain(small_cell)
     301             : 
     302             :       !-----------------------------------------------------------------------------
     303             :       !-----------------------------------------------------------------------------
     304             :       ! 3. Initialize atom coords from the bigger system
     305             :       !-----------------------------------------------------------------------------
     306         538 :       nat = SIZE(sub_atom_index)
     307         538 :       topology%natoms = nat
     308         538 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%r))
     309         538 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_atmname))
     310         538 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_molname))
     311         538 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_resname))
     312         538 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_mass))
     313         538 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_charge))
     314             :       ALLOCATE (topology%atom_info%r(3, nat), topology%atom_info%id_atmname(nat), &
     315             :                 topology%atom_info%id_molname(nat), topology%atom_info%id_resname(nat), &
     316             :                 topology%atom_info%id_element(nat), topology%atom_info%atm_mass(nat), &
     317        5918 :                 topology%atom_info%atm_charge(nat))
     318             : 
     319         538 :       CALL cp_subsys_get(big_subsys, particles=particles)
     320        4250 :       DO iat = 1, nat
     321       14848 :          topology%atom_info%r(:, iat) = particles%els(sub_atom_index(iat))%r
     322        3712 :          topology%atom_info%id_atmname(iat) = str2id(s2s(sub_atom_kind_name(iat)))
     323        3712 :          topology%atom_info%id_molname(iat) = topology%atom_info%id_atmname(iat)
     324        3712 :          topology%atom_info%id_resname(iat) = topology%atom_info%id_atmname(iat)
     325             :          !
     326             :          ! Defining element
     327             :          !
     328        3712 :          id_ = INDEX(id2str(topology%atom_info%id_atmname(iat)), "_") - 1
     329        3712 :          IF (id_ == -1) id_ = LEN_TRIM(id2str(topology%atom_info%id_atmname(iat)))
     330        3712 :          strtmp1 = id2str(topology%atom_info%id_atmname(iat))
     331        3712 :          strtmp1 = strtmp1(1:id_)
     332             :          CALL check_subsys_element(strtmp1, strtmp1, my_element, &
     333        3712 :                                    subsys_section, use_mm_map_first=.FALSE.)
     334        3712 :          topology%atom_info%id_element(iat) = str2id(s2s(my_element))
     335        3712 :          topology%atom_info%atm_mass(iat) = 0._dp
     336        4250 :          topology%atom_info%atm_charge(iat) = 0._dp
     337             :       END DO
     338         538 :       topology%conn_type = do_conn_off
     339             : 
     340             :       !-----------------------------------------------------------------------------
     341             :       !-----------------------------------------------------------------------------
     342             :       ! 4. Read in or generate the molecular connectivity
     343             :       !-----------------------------------------------------------------------------
     344             :       CALL connectivity_control(topology, para_env, subsys_section=subsys_section, &
     345         538 :                                 force_env_section=force_env_section)
     346             : 
     347             :       !-----------------------------------------------------------------------------
     348             :       !-----------------------------------------------------------------------------
     349             :       ! 5. Pack everything into the molecular types
     350             :       !-----------------------------------------------------------------------------
     351             :       CALL topology_connectivity_pack(molecule_kind_set, molecule_set, &
     352         538 :                                       topology, subsys_section=subsys_section)
     353             : 
     354             :       !-----------------------------------------------------------------------------
     355             :       !-----------------------------------------------------------------------------
     356             :       ! 6. Pack everything into the atomic types
     357             :       !-----------------------------------------------------------------------------
     358             :       CALL topology_coordinate_pack(particle_set, atomic_kind_set, &
     359             :                                     molecule_kind_set, molecule_set, topology, subsys_section=subsys_section, &
     360         538 :                                     force_env_section=force_env_section, ignore_outside_box=ignore_outside_box)
     361             : 
     362             :       !-----------------------------------------------------------------------------
     363             :       !-----------------------------------------------------------------------------
     364             :       ! 7. Cleanup the topology structure type
     365             :       !-----------------------------------------------------------------------------
     366         538 :       CALL deallocate_topology(topology)
     367             : 
     368             :       !-----------------------------------------------------------------------------
     369             :       !-----------------------------------------------------------------------------
     370             :       ! 8. Allocate new subsys
     371             :       !-----------------------------------------------------------------------------
     372        4842 :       ALLOCATE (small_subsys)
     373         538 :       CALL para_env%retain()
     374         538 :       small_subsys%para_env => para_env
     375         538 :       CALL particle_list_create(particles, els_ptr=particle_set)
     376         538 :       CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
     377         538 :       CALL molecule_list_create(mols, els_ptr=molecule_set)
     378         538 :       CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
     379             :       CALL cp_subsys_set(small_subsys, particles=particles, atomic_kinds=atomic_kinds, &
     380         538 :                          molecules=mols, molecule_kinds=mol_kinds)
     381         538 :       CALL particle_list_release(particles)
     382         538 :       CALL atomic_kind_list_release(atomic_kinds)
     383         538 :       CALL molecule_list_release(mols)
     384         538 :       CALL molecule_kind_list_release(mol_kinds)
     385             : 
     386      123202 :       ALLOCATE (small_subsys%virial)
     387         538 :       CALL atprop_create(small_subsys%atprop)
     388         538 :       CALL cp_result_create(small_subsys%results)
     389         538 :    END SUBROUTINE create_small_subsys
     390             : 
     391             : END MODULE cp_subsys_methods