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 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       35372 :    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        8843 :       INTEGER, DIMENSION(:), POINTER                     :: seed_vals
     109             :       LOGICAL                                            :: atomic_energy, atomic_stress, &
     110             :                                                             my_use_motion_section, &
     111             :                                                             pv_availability, pv_diagonal, &
     112             :                                                             pv_numerical
     113             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     114        8843 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     115             :       TYPE(molecule_kind_list_type), POINTER             :: mol_kinds
     116        8843 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     117             :       TYPE(molecule_list_type), POINTER                  :: mols
     118        8843 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     119             :       TYPE(particle_list_type), POINTER                  :: particles
     120        8843 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     121             :       TYPE(section_vals_type), POINTER                   :: colvar_section, my_force_env_section, &
     122             :                                                             my_subsys_section
     123             : 
     124           0 :       CPASSERT(.NOT. ASSOCIATED(subsys))
     125       88430 :       ALLOCATE (subsys)
     126             : 
     127        8843 :       CALL para_env%retain()
     128        8843 :       subsys%para_env => para_env
     129             : 
     130        8843 :       my_use_motion_section = .FALSE.
     131        8843 :       IF (PRESENT(use_motion_section)) &
     132        8841 :          my_use_motion_section = use_motion_section
     133             : 
     134        8843 :       my_force_env_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
     135        8843 :       IF (PRESENT(force_env_section)) &
     136        8841 :          my_force_env_section => force_env_section
     137             : 
     138        8843 :       my_subsys_section => section_vals_get_subs_vals(my_force_env_section, "SUBSYS")
     139        8843 :       IF (PRESENT(subsys_section)) &
     140        8687 :          my_subsys_section => subsys_section
     141             : 
     142        8843 :       CALL section_vals_val_get(my_subsys_section, "SEED", i_vals=seed_vals)
     143        8843 :       IF (SIZE(seed_vals) == 1) THEN
     144       79551 :          subsys%seed(:, :) = REAL(seed_vals(1), KIND=dp)
     145           4 :       ELSE IF (SIZE(seed_vals) == 6) THEN
     146          60 :          subsys%seed(1:3, 1:2) = RESHAPE(REAL(seed_vals(:), KIND=dp), (/3, 2/))
     147             :       ELSE
     148           0 :          CPABORT("Supply exactly 1 or 6 arguments for SEED in &SUBSYS only!")
     149             :       END IF
     150             : 
     151        8843 :       colvar_section => section_vals_get_subs_vals(my_subsys_section, "COLVAR")
     152             : 
     153        8843 :       CALL cp_subsys_read_colvar(subsys, colvar_section)
     154             : 
     155             :       !   *** Read the particle coordinates and allocate the atomic kind, ***
     156             :       !   *** the molecule kind, and the molecule data structures         ***
     157             :       CALL topology_control(atomic_kind_set, particle_set, molecule_kind_set, molecule_set, &
     158             :                             subsys%colvar_p, subsys%gci, root_section, para_env, &
     159             :                             force_env_section=my_force_env_section, &
     160             :                             subsys_section=my_subsys_section, use_motion_section=my_use_motion_section, &
     161       15047 :                             qmmm=qmmm, qmmm_env=qmmm_env, exclusions=exclusions, elkind=elkind)
     162             : 
     163        8843 :       CALL particle_list_create(particles, els_ptr=particle_set)
     164        8843 :       CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
     165        8843 :       CALL molecule_list_create(mols, els_ptr=molecule_set)
     166        8843 :       CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
     167             : 
     168             :       CALL cp_subsys_set(subsys, particles=particles, atomic_kinds=atomic_kinds, &
     169        8843 :                          molecules=mols, molecule_kinds=mol_kinds)
     170             : 
     171        8843 :       CALL particle_list_release(particles)
     172        8843 :       CALL atomic_kind_list_release(atomic_kinds)
     173        8843 :       CALL molecule_list_release(mols)
     174        8843 :       CALL molecule_kind_list_release(mol_kinds)
     175             : 
     176             :       ! Should we compute the virial?
     177        8843 :       CALL section_vals_val_get(my_force_env_section, "STRESS_TENSOR", i_val=stress_tensor)
     178        8077 :       SELECT CASE (stress_tensor)
     179             :       CASE (do_stress_none)
     180        8077 :          pv_availability = .FALSE.
     181        8077 :          pv_numerical = .FALSE.
     182        8077 :          pv_diagonal = .FALSE.
     183             :       CASE (do_stress_analytical)
     184         754 :          pv_availability = .TRUE.
     185         754 :          pv_numerical = .FALSE.
     186         754 :          pv_diagonal = .FALSE.
     187             :       CASE (do_stress_numerical)
     188           2 :          pv_availability = .TRUE.
     189           2 :          pv_numerical = .TRUE.
     190           2 :          pv_diagonal = .FALSE.
     191             :       CASE (do_stress_diagonal_anal)
     192           0 :          pv_availability = .TRUE.
     193           0 :          pv_numerical = .FALSE.
     194           0 :          pv_diagonal = .TRUE.
     195             :       CASE (do_stress_diagonal_numer)
     196          10 :          pv_availability = .TRUE.
     197          10 :          pv_numerical = .TRUE.
     198        8843 :          pv_diagonal = .TRUE.
     199             :       END SELECT
     200             : 
     201     2193064 :       ALLOCATE (subsys%virial)
     202             :       CALL virial_set(virial=subsys%virial, &
     203             :                       pv_availability=pv_availability, &
     204             :                       pv_numer=pv_numerical, &
     205        8843 :                       pv_diagonal=pv_diagonal)
     206             : 
     207             :       ! Should we compute atomic properties?
     208        8843 :       CALL atprop_create(subsys%atprop)
     209        8843 :       CALL section_vals_val_get(my_force_env_section, "PROPERTIES%ATOMIC%ENERGY", l_val=atomic_energy)
     210        8843 :       subsys%atprop%energy = atomic_energy
     211        8843 :       CALL section_vals_val_get(my_force_env_section, "PROPERTIES%ATOMIC%PRESSURE", l_val=atomic_stress)
     212        8843 :       IF (atomic_stress) THEN
     213          22 :          CPASSERT(pv_availability)
     214          22 :          CPASSERT(.NOT. pv_numerical)
     215             :       END IF
     216        8843 :       subsys%atprop%stress = atomic_stress
     217             : 
     218        8843 :       CALL cp_result_create(subsys%results)
     219        8843 :    END SUBROUTINE cp_subsys_create
     220             : 
     221             : ! **************************************************************************************************
     222             : !> \brief reads the colvar section of the colvar
     223             : !> \param subsys ...
     224             : !> \param colvar_section ...
     225             : !> \par History
     226             : !>      2006.01 Joost VandeVondele
     227             : ! **************************************************************************************************
     228        8843 :    SUBROUTINE cp_subsys_read_colvar(subsys, colvar_section)
     229             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     230             :       TYPE(section_vals_type), POINTER                   :: colvar_section
     231             : 
     232             :       INTEGER                                            :: ig, ncol
     233             : 
     234        8843 :       CALL section_vals_get(colvar_section, n_repetition=ncol)
     235       18420 :       ALLOCATE (subsys%colvar_p(ncol))
     236        9305 :       DO ig = 1, ncol
     237         462 :          NULLIFY (subsys%colvar_p(ig)%colvar)
     238        9305 :          CALL colvar_read(subsys%colvar_p(ig)%colvar, ig, colvar_section, subsys%para_env)
     239             :       END DO
     240        8843 :    END SUBROUTINE cp_subsys_read_colvar
     241             : 
     242             : ! **************************************************************************************************
     243             : !> \brief updates the molecule information of the given subsys
     244             : !> \param small_subsys the subsys to create
     245             : !> \param big_subsys the superset of small_subsys
     246             : !> \param small_cell ...
     247             : !> \param small_para_env the parallel environment for the new (small)
     248             : !>        subsys
     249             : !> \param sub_atom_index indexes of the atoms that should be in small_subsys
     250             : !> \param sub_atom_kind_name ...
     251             : !> \param para_env ...
     252             : !> \param force_env_section ...
     253             : !> \param subsys_section ...
     254             : !> \param ignore_outside_box ...
     255             : !> \par History
     256             : !>      05.2004 created [fawzi]
     257             : !> \author Fawzi Mohamed, Teodoro Laino
     258             : !> \note
     259             : !>      not really ready to be used with different para_envs for the small
     260             : !>      and big part
     261             : ! **************************************************************************************************
     262         528 :    SUBROUTINE create_small_subsys(small_subsys, big_subsys, small_cell, &
     263         528 :                                   small_para_env, sub_atom_index, sub_atom_kind_name, &
     264             :                                   para_env, force_env_section, subsys_section, ignore_outside_box)
     265             : 
     266             :       TYPE(cp_subsys_type), POINTER                      :: small_subsys, big_subsys
     267             :       TYPE(cell_type), POINTER                           :: small_cell
     268             :       TYPE(mp_para_env_type), POINTER                    :: small_para_env
     269             :       INTEGER, DIMENSION(:), INTENT(in)                  :: sub_atom_index
     270             :       CHARACTER(len=default_string_length), &
     271             :          DIMENSION(:), INTENT(in)                        :: sub_atom_kind_name
     272             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     273             :       TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
     274             :       LOGICAL, INTENT(in), OPTIONAL                      :: ignore_outside_box
     275             : 
     276             :       CHARACTER(len=default_string_length)               :: my_element, strtmp1
     277             :       INTEGER                                            :: iat, id_, nat
     278             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     279         528 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     280             :       TYPE(molecule_kind_list_type), POINTER             :: mol_kinds
     281         528 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     282             :       TYPE(molecule_list_type), POINTER                  :: mols
     283         528 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     284             :       TYPE(particle_list_type), POINTER                  :: particles
     285         528 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     286             :       TYPE(topology_parameters_type)                     :: topology
     287             : 
     288         528 :       NULLIFY (mol_kinds, mols, particles, atomic_kinds, atomic_kind_set, particle_set, &
     289         528 :                molecule_kind_set, molecule_set, particles, atomic_kinds)
     290             : 
     291           0 :       CPASSERT(.NOT. ASSOCIATED(small_subsys))
     292         528 :       CPASSERT(ASSOCIATED(big_subsys))
     293         528 :       IF (big_subsys%para_env /= small_para_env) &
     294           0 :          CPABORT("big_subsys%para_env==small_para_env")
     295             : 
     296             :       !-----------------------------------------------------------------------------
     297             :       !-----------------------------------------------------------------------------
     298             :       ! 1. Initialize the topology structure type
     299             :       !-----------------------------------------------------------------------------
     300         528 :       CALL init_topology(topology)
     301             : 
     302             :       !-----------------------------------------------------------------------------
     303             :       !-----------------------------------------------------------------------------
     304             :       ! 2. Get the cell info
     305             :       !-----------------------------------------------------------------------------
     306         528 :       topology%cell => small_cell
     307         528 :       CALL cell_retain(small_cell)
     308             : 
     309             :       !-----------------------------------------------------------------------------
     310             :       !-----------------------------------------------------------------------------
     311             :       ! 3. Initialize atom coords from the bigger system
     312             :       !-----------------------------------------------------------------------------
     313         528 :       nat = SIZE(sub_atom_index)
     314         528 :       topology%natoms = nat
     315         528 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%r))
     316         528 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_atmname))
     317         528 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_molname))
     318         528 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_resname))
     319         528 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_mass))
     320         528 :       CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_charge))
     321             :       ALLOCATE (topology%atom_info%r(3, nat), topology%atom_info%id_atmname(nat), &
     322             :                 topology%atom_info%id_molname(nat), topology%atom_info%id_resname(nat), &
     323             :                 topology%atom_info%id_element(nat), topology%atom_info%atm_mass(nat), &
     324        7920 :                 topology%atom_info%atm_charge(nat))
     325             : 
     326         528 :       CALL cp_subsys_get(big_subsys, particles=particles)
     327        4192 :       DO iat = 1, nat
     328       14656 :          topology%atom_info%r(:, iat) = particles%els(sub_atom_index(iat))%r
     329        3664 :          topology%atom_info%id_atmname(iat) = str2id(s2s(sub_atom_kind_name(iat)))
     330        3664 :          topology%atom_info%id_molname(iat) = topology%atom_info%id_atmname(iat)
     331        3664 :          topology%atom_info%id_resname(iat) = topology%atom_info%id_atmname(iat)
     332             :          !
     333             :          ! Defining element
     334             :          !
     335        3664 :          id_ = INDEX(id2str(topology%atom_info%id_atmname(iat)), "_") - 1
     336        3664 :          IF (id_ == -1) id_ = LEN_TRIM(id2str(topology%atom_info%id_atmname(iat)))
     337        3664 :          strtmp1 = id2str(topology%atom_info%id_atmname(iat))
     338        3664 :          strtmp1 = strtmp1(1:id_)
     339             :          CALL check_subsys_element(strtmp1, strtmp1, my_element, &
     340        3664 :                                    subsys_section, use_mm_map_first=.FALSE.)
     341        3664 :          topology%atom_info%id_element(iat) = str2id(s2s(my_element))
     342        3664 :          topology%atom_info%atm_mass(iat) = 0._dp
     343        4192 :          topology%atom_info%atm_charge(iat) = 0._dp
     344             :       END DO
     345         528 :       topology%conn_type = do_conn_off
     346             : 
     347             :       !-----------------------------------------------------------------------------
     348             :       !-----------------------------------------------------------------------------
     349             :       ! 4. Read in or generate the molecular connectivity
     350             :       !-----------------------------------------------------------------------------
     351             :       CALL connectivity_control(topology, para_env, subsys_section=subsys_section, &
     352         528 :                                 force_env_section=force_env_section)
     353             : 
     354             :       !-----------------------------------------------------------------------------
     355             :       !-----------------------------------------------------------------------------
     356             :       ! 5. Pack everything into the molecular types
     357             :       !-----------------------------------------------------------------------------
     358             :       CALL topology_connectivity_pack(molecule_kind_set, molecule_set, &
     359         528 :                                       topology, subsys_section=subsys_section)
     360             : 
     361             :       !-----------------------------------------------------------------------------
     362             :       !-----------------------------------------------------------------------------
     363             :       ! 6. Pack everything into the atomic types
     364             :       !-----------------------------------------------------------------------------
     365             :       CALL topology_coordinate_pack(particle_set, atomic_kind_set, &
     366             :                                     molecule_kind_set, molecule_set, topology, subsys_section=subsys_section, &
     367         528 :                                     force_env_section=force_env_section, ignore_outside_box=ignore_outside_box)
     368             : 
     369             :       !-----------------------------------------------------------------------------
     370             :       !-----------------------------------------------------------------------------
     371             :       ! 7. Cleanup the topology structure type
     372             :       !-----------------------------------------------------------------------------
     373         528 :       CALL deallocate_topology(topology)
     374             : 
     375             :       !-----------------------------------------------------------------------------
     376             :       !-----------------------------------------------------------------------------
     377             :       ! 8. Allocate new subsys
     378             :       !-----------------------------------------------------------------------------
     379        4752 :       ALLOCATE (small_subsys)
     380         528 :       CALL para_env%retain()
     381         528 :       small_subsys%para_env => para_env
     382         528 :       CALL particle_list_create(particles, els_ptr=particle_set)
     383         528 :       CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
     384         528 :       CALL molecule_list_create(mols, els_ptr=molecule_set)
     385         528 :       CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
     386             :       CALL cp_subsys_set(small_subsys, particles=particles, atomic_kinds=atomic_kinds, &
     387         528 :                          molecules=mols, molecule_kinds=mol_kinds)
     388         528 :       CALL particle_list_release(particles)
     389         528 :       CALL atomic_kind_list_release(atomic_kinds)
     390         528 :       CALL molecule_list_release(mols)
     391         528 :       CALL molecule_kind_list_release(mol_kinds)
     392             : 
     393      120912 :       ALLOCATE (small_subsys%virial)
     394         528 :       CALL atprop_create(small_subsys%atprop)
     395         528 :       CALL cp_result_create(small_subsys%results)
     396         528 :    END SUBROUTINE create_small_subsys
     397             : 
     398             : END MODULE cp_subsys_methods