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