1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> cjm, Feb-20-2001 : added all the extended variables to 9!> system_type 10!> gt 23-09-2002 : major changes. Pointer part is allocated/deallocated 11!> and initialized here. Atomic coordinates can now be 12!> read also from &COORD section in the input file. 13!> If &COORD is not found, .dat file is read. 14!> If & coord is found and .NOT. 'INIT', parsing of the .dat 15!> is performed to get the proper coords/vel/eta variables 16!> CJM 31-7-03 : Major rewrite. No more atype 17! ************************************************************************************************** 18MODULE atoms_input 19 USE atomic_kind_types, ONLY: atomic_kind_type,& 20 get_atomic_kind 21 USE cell_types, ONLY: cell_type,& 22 pbc,& 23 scaled_to_real 24 USE cp_linked_list_input, ONLY: cp_sll_val_next,& 25 cp_sll_val_type 26 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,& 27 cp_to_string 28 USE cp_parser_methods, ONLY: read_float_object 29 USE cp_units, ONLY: cp_unit_to_cp2k 30 USE input_section_types, ONLY: section_vals_get,& 31 section_vals_get_subs_vals,& 32 section_vals_list_get,& 33 section_vals_remove_values,& 34 section_vals_type,& 35 section_vals_val_get 36 USE input_val_types, ONLY: val_get,& 37 val_type 38 USE kinds, ONLY: default_string_length,& 39 dp 40 USE memory_utilities, ONLY: reallocate 41 USE particle_types, ONLY: particle_type 42 USE shell_potential_types, ONLY: shell_kind_type 43 USE string_table, ONLY: id2str,& 44 s2s,& 45 str2id 46 USE string_utilities, ONLY: uppercase 47 USE topology_types, ONLY: atom_info_type,& 48 topology_parameters_type 49#include "./base/base_uses.f90" 50 51 IMPLICIT NONE 52 53 PRIVATE 54 PUBLIC :: read_atoms_input, read_shell_coord_input 55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atoms_input' 56 57CONTAINS 58 59! ************************************************************************************************** 60!> \brief ... 61!> \param topology ... 62!> \param overwrite ... 63!> \param subsys_section ... 64!> \param save_mem ... 65!> \author CJM 66! ************************************************************************************************** 67 SUBROUTINE read_atoms_input(topology, overwrite, subsys_section, save_mem) 68 69 TYPE(topology_parameters_type) :: topology 70 LOGICAL, INTENT(IN), OPTIONAL :: overwrite 71 TYPE(section_vals_type), POINTER :: subsys_section 72 LOGICAL, INTENT(IN), OPTIONAL :: save_mem 73 74 CHARACTER(len=*), PARAMETER :: routineN = 'read_atoms_input', & 75 routineP = moduleN//':'//routineN 76 77 CHARACTER(len=2*default_string_length) :: line_att 78 CHARACTER(len=default_string_length) :: error_message, my_default_index, strtmp, & 79 unit_str 80 INTEGER :: default_id, end_c, handle, iatom, j, & 81 natom, output_unit, start_c, wrd 82 LOGICAL :: explicit, is_ok, my_overwrite, & 83 my_save_mem, scaled_coordinates 84 REAL(KIND=dp) :: r0(3), unit_conv 85 TYPE(atom_info_type), POINTER :: atom_info 86 TYPE(cell_type), POINTER :: cell 87 TYPE(cp_sll_val_type), POINTER :: list 88 TYPE(section_vals_type), POINTER :: coord_section 89 TYPE(val_type), POINTER :: val 90 91 my_overwrite = .FALSE. 92 my_save_mem = .FALSE. 93 error_message = "" 94 output_unit = cp_logger_get_default_io_unit() 95 IF (PRESENT(overwrite)) my_overwrite = overwrite 96 IF (PRESENT(save_mem)) my_save_mem = save_mem 97 NULLIFY (coord_section) 98 coord_section => section_vals_get_subs_vals(subsys_section, "COORD") 99 CALL section_vals_get(coord_section, explicit=explicit) 100 IF (.NOT. explicit) RETURN 101 102 CALL timeset(routineN, handle) 103 !----------------------------------------------------------------------------- 104 !----------------------------------------------------------------------------- 105 ! 1. get cell and topology%atom_info 106 !----------------------------------------------------------------------------- 107 atom_info => topology%atom_info 108 cell => topology%cell_muc 109 CALL section_vals_val_get(coord_section, "UNIT", c_val=unit_str) 110 CALL section_vals_val_get(coord_section, "SCALED", l_val=scaled_coordinates) 111 unit_conv = cp_unit_to_cp2k(1.0_dp, TRIM(unit_str)) 112 113 !----------------------------------------------------------------------------- 114 !----------------------------------------------------------------------------- 115 ! 2. Read in the coordinates from &COORD section in the input file 116 !----------------------------------------------------------------------------- 117 CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", & 118 n_rep_val=natom) 119 topology%natoms = natom 120 IF (my_overwrite) THEN 121 CPASSERT(SIZE(atom_info%r, 2) == natom) 122 CALL cp_warn(__LOCATION__, & 123 "Overwriting coordinates. Active coordinates read from &COORD section."// & 124 " Active coordinates READ from &COORD section ") 125 CALL section_vals_list_get(coord_section, "_DEFAULT_KEYWORD_", list=list) 126 DO iatom = 1, natom 127 is_ok = cp_sll_val_next(list, val) 128 CALL val_get(val, c_val=line_att) 129 ! Read name and atomic coordinates 130 start_c = 1 131 DO wrd = 1, 4 132 DO j = start_c, LEN(line_att) 133 IF (line_att(j:j) /= ' ') THEN 134 start_c = j 135 EXIT 136 END IF 137 END DO 138 end_c = LEN(line_att) + 1 139 DO j = start_c, LEN(line_att) 140 IF (line_att(j:j) == ' ') THEN 141 end_c = j 142 EXIT 143 END IF 144 END DO 145 IF (LEN_TRIM(line_att(start_c:end_c - 1)) == 0) & 146 CPABORT("incorrectly formatted line in coord section'"//line_att//"'") 147 IF (wrd == 1) THEN 148 atom_info%id_atmname(iatom) = str2id(s2s(line_att(start_c:end_c - 1))) 149 ELSE 150 READ (line_att(start_c:end_c - 1), *) atom_info%r(wrd - 1, iatom) 151 END IF 152 start_c = end_c 153 END DO 154 END DO 155 ELSE 156 ! Element is assigned on the basis of the atm_name 157 topology%aa_element = .TRUE. 158 159 CALL reallocate(atom_info%id_molname, 1, natom) 160 CALL reallocate(atom_info%id_resname, 1, natom) 161 CALL reallocate(atom_info%resid, 1, natom) 162 CALL reallocate(atom_info%id_atmname, 1, natom) 163 CALL reallocate(atom_info%id_element, 1, natom) 164 CALL reallocate(atom_info%r, 1, 3, 1, natom) 165 CALL reallocate(atom_info%atm_mass, 1, natom) 166 CALL reallocate(atom_info%atm_charge, 1, natom) 167 168 CALL section_vals_list_get(coord_section, "_DEFAULT_KEYWORD_", list=list) 169 DO iatom = 1, natom 170 ! we use only the first default_string_length characters of each line 171 is_ok = cp_sll_val_next(list, val) 172 CALL val_get(val, c_val=line_att) 173 default_id = str2id(s2s("")) 174 atom_info%id_molname(iatom) = default_id 175 atom_info%id_resname(iatom) = default_id 176 atom_info%resid(iatom) = 1 177 atom_info%id_atmname(iatom) = default_id 178 atom_info%id_element(iatom) = default_id 179 topology%molname_generated = .TRUE. 180 ! Read name and atomic coordinates 181 start_c = 1 182 DO wrd = 1, 6 183 DO j = start_c, LEN(line_att) 184 IF (line_att(j:j) /= ' ') THEN 185 start_c = j 186 EXIT 187 END IF 188 END DO 189 end_c = LEN(line_att) + 1 190 DO j = start_c, LEN(line_att) 191 IF (line_att(j:j) == ' ') THEN 192 end_c = j 193 EXIT 194 END IF 195 END DO 196 IF (LEN_TRIM(line_att(start_c:end_c - 1)) == 0) & 197 CALL cp_abort(__LOCATION__, & 198 "Incorrectly formatted input line for atom "// & 199 TRIM(ADJUSTL(cp_to_string(iatom)))// & 200 " found in COORD section. Input line: <"// & 201 TRIM(line_att)//"> ") 202 SELECT CASE (wrd) 203 CASE (1) 204 atom_info%id_atmname(iatom) = str2id(s2s(line_att(start_c:end_c - 1))) 205 CASE (2:4) 206 CALL read_float_object(line_att(start_c:end_c - 1), & 207 atom_info%r(wrd - 1, iatom), error_message) 208 IF (LEN_TRIM(error_message) /= 0) & 209 CALL cp_abort(__LOCATION__, & 210 "Incorrectly formatted input line for atom "// & 211 TRIM(ADJUSTL(cp_to_string(iatom)))// & 212 " found in COORD section. "//TRIM(error_message)// & 213 " Input line: <"//TRIM(line_att)//"> ") 214 CASE (5) 215 READ (line_att(start_c:end_c - 1), *) strtmp 216 atom_info%id_molname(iatom) = str2id(strtmp) 217 atom_info%id_resname(iatom) = atom_info%id_molname(iatom) 218 topology%molname_generated = .FALSE. 219 CASE (6) 220 READ (line_att(start_c:end_c - 1), *) strtmp 221 atom_info%id_resname(iatom) = str2id(strtmp) 222 END SELECT 223 start_c = end_c 224 IF (start_c > LEN_TRIM(line_att)) EXIT 225 END DO 226 IF (topology%molname_generated) THEN 227 ! Use defaults, if no molname was specified 228 WRITE (my_default_index, '(I0)') iatom 229 atom_info%id_molname(iatom) = str2id(s2s(TRIM(id2str(atom_info%id_atmname(iatom)))//TRIM(my_default_index))) 230 atom_info%id_resname(iatom) = atom_info%id_molname(iatom) 231 END IF 232 atom_info%id_element(iatom) = atom_info%id_atmname(iatom) 233 atom_info%atm_mass(iatom) = 0.0_dp 234 atom_info%atm_charge(iatom) = -HUGE(0.0_dp) 235 END DO 236 END IF 237 !----------------------------------------------------------------------------- 238 !----------------------------------------------------------------------------- 239 ! 3. Convert coordinates into internal cp2k coordinates 240 !----------------------------------------------------------------------------- 241 DO iatom = 1, natom 242 IF (scaled_coordinates) THEN 243 r0 = atom_info%r(:, iatom) 244 CALL scaled_to_real(atom_info%r(:, iatom), r0, cell) 245 ELSE 246 atom_info%r(:, iatom) = atom_info%r(:, iatom)*unit_conv 247 END IF 248 END DO 249 IF (my_save_mem) CALL section_vals_remove_values(coord_section) 250 251 CALL timestop(handle) 252 END SUBROUTINE read_atoms_input 253 254! ************************************************************************************************** 255!> \brief ... 256!> \param particle_set ... 257!> \param shell_particle_set ... 258!> \param cell ... 259!> \param subsys_section ... 260!> \param core_particle_set ... 261!> \param save_mem ... 262!> \author MI 263! ************************************************************************************************** 264 SUBROUTINE read_shell_coord_input(particle_set, shell_particle_set, cell, & 265 subsys_section, core_particle_set, save_mem) 266 267 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, shell_particle_set 268 TYPE(cell_type), POINTER :: cell 269 TYPE(section_vals_type), POINTER :: subsys_section 270 TYPE(particle_type), DIMENSION(:), OPTIONAL, & 271 POINTER :: core_particle_set 272 LOGICAL, INTENT(IN), OPTIONAL :: save_mem 273 274 CHARACTER(len=*), PARAMETER :: routineN = 'read_shell_coord_input', & 275 routineP = moduleN//':'//routineN 276 277 CHARACTER(len=2*default_string_length) :: line_att 278 CHARACTER(len=default_string_length) :: name_kind, unit_str 279 CHARACTER(len=default_string_length), & 280 ALLOCATABLE, DIMENSION(:) :: at_name, at_name_c 281 INTEGER :: end_c, handle, ishell, j, nshell, & 282 output_unit, sh_index, start_c, wrd 283 INTEGER, ALLOCATABLE, DIMENSION(:) :: at_index, at_index_c 284 LOGICAL :: core_scaled_coordinates, explicit, & 285 is_ok, is_shell, my_save_mem, & 286 shell_scaled_coordinates 287 REAL(KIND=dp) :: dab, mass_com, rab(3), unit_conv_core, & 288 unit_conv_shell 289 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r, rc 290 TYPE(atomic_kind_type), POINTER :: atomic_kind 291 TYPE(cp_sll_val_type), POINTER :: list 292 TYPE(section_vals_type), POINTER :: core_coord_section, shell_coord_section 293 TYPE(shell_kind_type), POINTER :: shell 294 TYPE(val_type), POINTER :: val 295 296 my_save_mem = .FALSE. 297 NULLIFY (atomic_kind, list, shell_coord_section, shell, val) 298 output_unit = cp_logger_get_default_io_unit() 299 300 IF (PRESENT(save_mem)) my_save_mem = save_mem 301 NULLIFY (shell_coord_section, core_coord_section) 302 shell_coord_section => section_vals_get_subs_vals(subsys_section, "SHELL_COORD") 303 CALL section_vals_get(shell_coord_section, explicit=explicit) 304 IF (.NOT. explicit) RETURN 305 306 CALL timeset(routineN, handle) 307 CPASSERT(ASSOCIATED(particle_set)) 308 !----------------------------------------------------------------------------- 309 !----------------------------------------------------------------------------- 310 ! 2. Read in the coordinates from &SHELL_COORD section in the input file 311 !----------------------------------------------------------------------------- 312 CALL section_vals_val_get(shell_coord_section, "UNIT", c_val=unit_str) 313 CALL section_vals_val_get(shell_coord_section, "SCALED", l_val=shell_scaled_coordinates) 314 unit_conv_shell = cp_unit_to_cp2k(1.0_dp, TRIM(unit_str)) 315 CALL section_vals_val_get(shell_coord_section, "_DEFAULT_KEYWORD_", & 316 n_rep_val=nshell) 317 318 IF (ASSOCIATED(shell_particle_set)) THEN 319 CPASSERT((SIZE(shell_particle_set, 1) == nshell)) 320 ALLOCATE (r(3, nshell), at_name(nshell), at_index(nshell)) 321 CALL cp_warn(__LOCATION__, & 322 "Overwriting shell coordinates. "// & 323 "Active coordinates READ from &SHELL_COORD section. ") 324 CALL section_vals_list_get(shell_coord_section, "_DEFAULT_KEYWORD_", list=list) 325 DO ishell = 1, nshell 326 ! we use only the first default_string_length characters of each line 327 is_ok = cp_sll_val_next(list, val) 328 CALL val_get(val, c_val=line_att) 329 start_c = 1 330 DO wrd = 1, 5 331 DO j = start_c, LEN(line_att) 332 IF (line_att(j:j) /= ' ') THEN 333 start_c = j 334 EXIT 335 END IF 336 END DO 337 end_c = LEN(line_att) + 1 338 DO j = start_c, LEN(line_att) 339 IF (line_att(j:j) == ' ') THEN 340 end_c = j 341 EXIT 342 END IF 343 END DO 344 IF (wrd /= 5 .AND. end_c >= LEN(line_att) + 1) & 345 CPABORT("incorrectly formatted line in coord section'"//line_att//"'") 346 IF (wrd == 1) THEN 347 at_name(ishell) = line_att(start_c:end_c - 1) 348 CALL uppercase(at_name(ishell)) 349 ELSE IF (wrd == 5) THEN 350 READ (line_att(start_c:end_c - 1), *) at_index(ishell) 351 ELSE 352 READ (line_att(start_c:end_c - 1), *) r(wrd - 1, ishell) 353 END IF 354 start_c = end_c 355 END DO 356 END DO 357 358 IF (PRESENT(core_particle_set)) THEN 359 CPASSERT(ASSOCIATED(core_particle_set)) 360 core_coord_section => section_vals_get_subs_vals(subsys_section, "CORE_COORD") 361 CALL section_vals_get(core_coord_section, explicit=explicit) 362 IF (explicit) THEN 363 CALL section_vals_val_get(core_coord_section, "UNIT", c_val=unit_str) 364 CALL section_vals_val_get(core_coord_section, "SCALED", l_val=core_scaled_coordinates) 365 unit_conv_core = cp_unit_to_cp2k(1.0_dp, TRIM(unit_str)) 366 CALL section_vals_val_get(core_coord_section, "_DEFAULT_KEYWORD_", & 367 n_rep_val=nshell) 368 369 CPASSERT((SIZE(core_particle_set, 1) == nshell)) 370 ALLOCATE (rc(3, nshell), at_name_c(nshell), at_index_c(nshell)) 371 CALL cp_warn(__LOCATION__, & 372 "Overwriting cores coordinates. "// & 373 " Active coordinates READ from &CORE_COORD section. ") 374 CALL section_vals_list_get(core_coord_section, "_DEFAULT_KEYWORD_", list=list) 375 DO ishell = 1, nshell 376 ! we use only the first default_string_length characters of each line 377 is_ok = cp_sll_val_next(list, val) 378 CALL val_get(val, c_val=line_att) 379 start_c = 1 380 DO wrd = 1, 5 381 DO j = start_c, LEN(line_att) 382 IF (line_att(j:j) /= ' ') THEN 383 start_c = j 384 EXIT 385 END IF 386 END DO 387 end_c = LEN(line_att) + 1 388 DO j = start_c, LEN(line_att) 389 IF (line_att(j:j) == ' ') THEN 390 end_c = j 391 EXIT 392 END IF 393 END DO 394 IF (wrd /= 5 .AND. end_c >= LEN(line_att) + 1) & 395 CPABORT("incorrectly formatted line in coord section'"//line_att//"'") 396 IF (wrd == 1) THEN 397 at_name_c(ishell) = line_att(start_c:end_c - 1) 398 CALL uppercase(at_name_c(ishell)) 399 ELSE IF (wrd == 5) THEN 400 READ (line_att(start_c:end_c - 1), *) at_index_c(ishell) 401 ELSE 402 READ (line_att(start_c:end_c - 1), *) rc(wrd - 1, ishell) 403 END IF 404 start_c = end_c 405 END DO 406 END DO 407 IF (my_save_mem) CALL section_vals_remove_values(core_coord_section) 408 END IF ! explicit 409 END IF ! core_particle_set 410 411 !----------------------------------------------------------------------------- 412 ! 3. Check corrispondence and convert coordinates into internal cp2k coordinates 413 !----------------------------------------------------------------------------- 414 DO ishell = 1, nshell 415 atomic_kind => particle_set(at_index(ishell))%atomic_kind 416 CALL get_atomic_kind(atomic_kind=atomic_kind, & 417 name=name_kind, shell_active=is_shell, mass=mass_com, shell=shell) 418 CALL uppercase(name_kind) 419 IF ((TRIM(at_name(ishell)) == TRIM(name_kind)) .AND. is_shell) THEN 420 sh_index = particle_set(at_index(ishell))%shell_index 421 IF (shell_scaled_coordinates) THEN 422 CALL scaled_to_real(r(:, ishell), shell_particle_set(sh_index)%r(:), cell) 423 ELSE 424 shell_particle_set(sh_index)%r(:) = r(:, ishell)*unit_conv_shell 425 END IF 426 shell_particle_set(sh_index)%atom_index = at_index(ishell) 427 428 IF (PRESENT(core_particle_set) .AND. .NOT. explicit) THEN 429 core_particle_set(sh_index)%r(1) = (mass_com*particle_set(at_index(ishell))%r(1) - & 430 shell%mass_shell*shell_particle_set(sh_index)%r(1))/shell%mass_core 431 core_particle_set(sh_index)%r(2) = (mass_com*particle_set(at_index(ishell))%r(2) - & 432 shell%mass_shell*shell_particle_set(sh_index)%r(2))/shell%mass_core 433 core_particle_set(sh_index)%r(3) = (mass_com*particle_set(at_index(ishell))%r(3) - & 434 shell%mass_shell*shell_particle_set(sh_index)%r(3))/shell%mass_core 435 core_particle_set(sh_index)%atom_index = at_index(ishell) 436 rab = pbc(shell_particle_set(sh_index)%r, core_particle_set(sh_index)%r, cell) 437 ELSE IF (explicit) THEN 438 IF (core_scaled_coordinates) THEN 439 CALL scaled_to_real(rc(:, ishell), core_particle_set(sh_index)%r(:), cell) 440 ELSE 441 core_particle_set(sh_index)%r(:) = rc(:, ishell)*unit_conv_core 442 END IF 443 core_particle_set(sh_index)%atom_index = at_index_c(ishell) 444 rab = pbc(shell_particle_set(sh_index)%r, core_particle_set(sh_index)%r, cell) 445 CPASSERT(TRIM(at_name(ishell)) == TRIM(at_name_c(ishell))) 446 CPASSERT(at_index(ishell) == at_index_c(ishell)) 447 ELSE 448 rab = pbc(shell_particle_set(sh_index)%r, particle_set(at_index(ishell))%r, cell) 449 END IF 450 451 dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)) 452 IF (shell%max_dist > 0.0_dp .AND. shell%max_dist < dab) THEN 453 IF (output_unit > 0) THEN 454 WRITE (output_unit, *) "WARNING : shell and core for atom ", at_index(ishell), " seem to be too distant. " 455 END IF 456 END IF 457 458 ELSE 459 CPABORT("shell coordinate assigned to the wrong atom. check the shell indexes in the input") 460 END IF 461 END DO 462 DEALLOCATE (r, at_index, at_name) 463 DEALLOCATE (rc, at_index_c, at_name_c) 464 465 END IF 466 467 IF (my_save_mem) CALL section_vals_remove_values(shell_coord_section) 468 469 CALL timestop(handle) 470 471 END SUBROUTINE read_shell_coord_input 472 473END MODULE atoms_input 474