atoms_input.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \par History
!>    cjm, Feb-20-2001 : added all the extended variables to
!>    system_type
!>    gt 23-09-2002 : major changes. Pointer part is allocated/deallocated
!>                    and initialized here. Atomic coordinates can now be
!>                    read also from &COORD section in the input file.
!>                    If &COORD is not found, .dat file is read.
!>                    If & coord is found and .NOT. 'INIT', parsing of the .dat
!>                    is performed to get the proper coords/vel/eta variables
!>     CJM 31-7-03  : Major rewrite.  No more atype
! **************************************************************************************************
MODULE atoms_input
   USE atomic_kind_types,               ONLY: atomic_kind_type,&
                                              get_atomic_kind
   USE cell_types,                      ONLY: cell_transform_input_cartesian,&
                                              cell_type,&
                                              pbc,&
                                              scaled_to_real
   USE cp_linked_list_input,            ONLY: cp_sll_val_next,&
                                              cp_sll_val_type
   USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit,&
                                              cp_to_string
   USE cp_parser_methods,               ONLY: read_float_object
   USE cp_units,                        ONLY: cp_unit_to_cp2k
   USE input_section_types,             ONLY: section_vals_get,&
                                              section_vals_get_subs_vals,&
                                              section_vals_list_get,&
                                              section_vals_remove_values,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE input_val_types,                 ONLY: val_get,&
                                              val_type
   USE kinds,                           ONLY: default_string_length,&
                                              dp
   USE memory_utilities,                ONLY: reallocate
   USE particle_types,                  ONLY: particle_type
   USE shell_potential_types,           ONLY: shell_kind_type
   USE string_table,                    ONLY: id2str,&
                                              s2s,&
                                              str2id
   USE string_utilities,                ONLY: uppercase
   USE topology_types,                  ONLY: atom_info_type,&
                                              topology_parameters_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
   PUBLIC :: read_atoms_input, read_shell_coord_input
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atoms_input'
 
CONTAINS
 
! **************************************************************************************************
!> \brief ...
!> \param topology ...
!> \param overwrite ...
!> \param subsys_section ...
!> \param save_mem ...
!> \author CJM
! **************************************************************************************************
   SUBROUTINE read_atoms_input(topology, overwrite, subsys_section, save_mem)
 
      TYPE(topology_parameters_type)                     :: topology
      LOGICAL, INTENT(IN), OPTIONAL                      :: overwrite
      TYPE(section_vals_type), POINTER                   :: subsys_section
      LOGICAL, INTENT(IN), OPTIONAL                      :: save_mem
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'read_atoms_input'
 
      CHARACTER(len=2*default_string_length)             :: line_att
      CHARACTER(len=default_string_length)               :: error_message, my_default_index, strtmp, &
                                                            unit_str
      INTEGER                                            :: default_id, end_c, handle, iatom, j, &
                                                            natom, output_unit, start_c, wrd
      LOGICAL                                            :: explicit, is_ok, my_overwrite, &
                                                            my_save_mem, scaled_coordinates
      REAL(kind=dp)                                      :: r0(3), unit_conv
      TYPE(atom_info_type), POINTER                      :: atom_info
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_sll_val_type), POINTER                     :: list
      TYPE(section_vals_type), POINTER                   :: coord_section
      TYPE(val_type), POINTER                            :: val
 
      my_overwrite = .false.
      my_save_mem = .false.
      error_message = ""
      output_unit = cp_logger_get_default_io_unit()
      IF (PRESENT(overwrite)) my_overwrite = overwrite
      IF (PRESENT(save_mem)) my_save_mem = save_mem
      NULLIFY (coord_section)
      coord_section => section_vals_get_subs_vals(subsys_section, "COORD")
      CALL section_vals_get(coord_section, explicit=explicit)
      IF (.NOT. explicit) RETURN
 
      CALL timeset(routinen, handle)
      !-----------------------------------------------------------------------------
      !-----------------------------------------------------------------------------
      ! 1. get cell and topology%atom_info
      !-----------------------------------------------------------------------------
      atom_info => topology%atom_info
      cell => topology%cell_muc
      CALL section_vals_val_get(coord_section, "UNIT", c_val=unit_str)
      CALL section_vals_val_get(coord_section, "SCALED", l_val=scaled_coordinates)
      unit_conv = cp_unit_to_cp2k(1.0_dp, trim(unit_str))
 
      !-----------------------------------------------------------------------------
      !-----------------------------------------------------------------------------
      ! 2. Read in the coordinates from &COORD section in the input file
      !-----------------------------------------------------------------------------
      CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", &
                                n_rep_val=natom)
      topology%natoms = natom
      IF (my_overwrite) THEN
         cpassert(SIZE(atom_info%r, 2) == natom)
         CALL cp_warn(__location__, &
                      "Overwriting coordinates. Active coordinates read from &COORD section."// &
                      " Active coordinates READ from &COORD section ")
         CALL section_vals_list_get(coord_section, "_DEFAULT_KEYWORD_", list=list)
         DO iatom = 1, natom
            is_ok = cp_sll_val_next(list, val)
            CALL val_get(val, c_val=line_att)
            ! Read name and atomic coordinates
            start_c = 1
            DO wrd = 1, 4
               DO j = start_c, len(line_att)
                  IF (line_att(j:j) /= ' ') THEN
                     start_c = j
                     EXIT
                  END IF
               END DO
               end_c = len(line_att) + 1
               DO j = start_c, len(line_att)
                  IF (line_att(j:j) == ' ') THEN
                     end_c = j
                     EXIT
                  END IF
               END DO
               IF (len_trim(line_att(start_c:end_c - 1)) == 0) THEN
                  cpabort("incorrectly formatted line in coord section'"//line_att//"'")
               END IF
               IF (wrd == 1) THEN
                  atom_info%id_atmname(iatom) = str2id(s2s(line_att(start_c:end_c - 1)))
               ELSE
                  READ (line_att(start_c:end_c - 1), *) atom_info%r(wrd - 1, iatom)
               END IF
               start_c = end_c
            END DO
         END DO
      ELSE
         ! Element is assigned on the basis of the atm_name
         topology%aa_element = .true.
 
         CALL reallocate(atom_info%id_molname, 1, natom)
         CALL reallocate(atom_info%id_resname, 1, natom)
         CALL reallocate(atom_info%resid, 1, natom)
         CALL reallocate(atom_info%id_atmname, 1, natom)
         CALL reallocate(atom_info%id_element, 1, natom)
         CALL reallocate(atom_info%r, 1, 3, 1, natom)
         CALL reallocate(atom_info%atm_mass, 1, natom)
         CALL reallocate(atom_info%atm_charge, 1, natom)
 
         CALL section_vals_list_get(coord_section, "_DEFAULT_KEYWORD_", list=list)
         DO iatom = 1, natom
            ! we use only the first default_string_length characters of each line
            is_ok = cp_sll_val_next(list, val)
            CALL val_get(val, c_val=line_att)
            default_id = str2id(s2s(""))
            atom_info%id_molname(iatom) = default_id
            atom_info%id_resname(iatom) = default_id
            atom_info%resid(iatom) = 1
            atom_info%id_atmname(iatom) = default_id
            atom_info%id_element(iatom) = default_id
            topology%molname_generated = .true.
            ! Read name and atomic coordinates
            start_c = 1
            DO wrd = 1, 6
               DO j = start_c, len(line_att)
                  IF (line_att(j:j) /= ' ') THEN
                     start_c = j
                     EXIT
                  END IF
               END DO
               end_c = len(line_att) + 1
               DO j = start_c, len(line_att)
                  IF (line_att(j:j) == ' ') THEN
                     end_c = j
                     EXIT
                  END IF
               END DO
               IF (len_trim(line_att(start_c:end_c - 1)) == 0) THEN
                  CALL cp_abort(__location__, &
                                "Incorrectly formatted input line for atom "// &
                                trim(adjustl(cp_to_string(iatom)))// &
                                " found in COORD section. Input line: <"// &
                                trim(line_att)//"> ")
               END IF
               SELECT CASE (wrd)
               CASE (1)
                  atom_info%id_atmname(iatom) = str2id(s2s(line_att(start_c:end_c - 1)))
               CASE (2:4)
                  CALL read_float_object(line_att(start_c:end_c - 1), &
                                         atom_info%r(wrd - 1, iatom), error_message)
                  IF (len_trim(error_message) /= 0) THEN
                     CALL cp_abort(__location__, &
                                   "Incorrectly formatted input line for atom "// &
                                   trim(adjustl(cp_to_string(iatom)))// &
                                   " found in COORD section. "//trim(error_message)// &
                                   " Input line: <"//trim(line_att)//"> ")
                  END IF
               CASE (5)
                  READ (line_att(start_c:end_c - 1), *) strtmp
                  atom_info%id_molname(iatom) = str2id(strtmp)
                  atom_info%id_resname(iatom) = atom_info%id_molname(iatom)
                  topology%molname_generated = .false.
               CASE (6)
                  READ (line_att(start_c:end_c - 1), *) strtmp
                  atom_info%id_resname(iatom) = str2id(strtmp)
               END SELECT
               start_c = end_c
               IF (start_c > len_trim(line_att)) EXIT
            END DO
            IF (topology%molname_generated) THEN
               ! Use defaults, if no molname was specified
               WRITE (my_default_index, '(I0)') iatom
               atom_info%id_molname(iatom) = str2id(s2s(trim(id2str(atom_info%id_atmname(iatom)))//trim(my_default_index)))
               atom_info%id_resname(iatom) = atom_info%id_molname(iatom)
            END IF
            atom_info%id_element(iatom) = atom_info%id_atmname(iatom)
            atom_info%atm_mass(iatom) = 0.0_dp
            atom_info%atm_charge(iatom) = -huge(0.0_dp)
         END DO
      END IF
      !-----------------------------------------------------------------------------
      !-----------------------------------------------------------------------------
      ! 3. Convert coordinates into internal cp2k coordinates
      !-----------------------------------------------------------------------------
      DO iatom = 1, natom
         IF (scaled_coordinates) THEN
            r0 = atom_info%r(:, iatom)
            CALL scaled_to_real(atom_info%r(:, iatom), r0, cell)
         ELSE
            atom_info%r(:, iatom) = atom_info%r(:, iatom)*unit_conv
            CALL cell_transform_input_cartesian(cell, atom_info%r(:, iatom))
         END IF
      END DO
      IF (my_save_mem) CALL section_vals_remove_values(coord_section)
 
      CALL timestop(handle)
   END SUBROUTINE read_atoms_input
 
! **************************************************************************************************
!> \brief ...
!> \param particle_set ...
!> \param shell_particle_set ...
!> \param cell ...
!> \param subsys_section ...
!> \param core_particle_set ...
!> \param save_mem ...
!> \author MI
! **************************************************************************************************
   SUBROUTINE read_shell_coord_input(particle_set, shell_particle_set, cell, &
                                     subsys_section, core_particle_set, save_mem)
 
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set, shell_particle_set
      TYPE(cell_type), POINTER                           :: cell
      TYPE(section_vals_type), POINTER                   :: subsys_section
      TYPE(particle_type), DIMENSION(:), OPTIONAL, &
         POINTER                                         :: core_particle_set
      LOGICAL, INTENT(IN), OPTIONAL                      :: save_mem
 
      CHARACTER(len=*), PARAMETER :: routinen = 'read_shell_coord_input'
 
      CHARACTER(len=2*default_string_length)             :: line_att
      CHARACTER(len=default_string_length)               :: name_kind, unit_str
      CHARACTER(len=default_string_length), &
         ALLOCATABLE, DIMENSION(:)                       :: at_name, at_name_c
      INTEGER                                            :: end_c, handle, ishell, j, nshell, &
                                                            output_unit, sh_index, start_c, wrd
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: at_index, at_index_c
      LOGICAL                                            :: core_scaled_coordinates, explicit, &
                                                            is_ok, is_shell, my_save_mem, &
                                                            shell_scaled_coordinates
      REAL(kind=dp)                                      :: dab, mass_com, rab(3), unit_conv_core, &
                                                            unit_conv_shell
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: r, rc
      TYPE(atomic_kind_type), POINTER                    :: atomic_kind
      TYPE(cp_sll_val_type), POINTER                     :: list
      TYPE(section_vals_type), POINTER                   :: core_coord_section, shell_coord_section
      TYPE(shell_kind_type), POINTER                     :: shell
      TYPE(val_type), POINTER                            :: val
 
      my_save_mem = .false.
      NULLIFY (atomic_kind, list, shell_coord_section, shell, val)
      output_unit = cp_logger_get_default_io_unit()
 
      IF (PRESENT(save_mem)) my_save_mem = save_mem
      NULLIFY (shell_coord_section, core_coord_section)
      shell_coord_section => section_vals_get_subs_vals(subsys_section, "SHELL_COORD")
      CALL section_vals_get(shell_coord_section, explicit=explicit)
      IF (.NOT. explicit) RETURN
 
      CALL timeset(routinen, handle)
      cpassert(ASSOCIATED(particle_set))
      !-----------------------------------------------------------------------------
      !-----------------------------------------------------------------------------
      ! 2. Read in the coordinates from &SHELL_COORD section in the input file
      !-----------------------------------------------------------------------------
      CALL section_vals_val_get(shell_coord_section, "UNIT", c_val=unit_str)
      CALL section_vals_val_get(shell_coord_section, "SCALED", l_val=shell_scaled_coordinates)
      unit_conv_shell = cp_unit_to_cp2k(1.0_dp, trim(unit_str))
      CALL section_vals_val_get(shell_coord_section, "_DEFAULT_KEYWORD_", &
                                n_rep_val=nshell)
 
      IF (ASSOCIATED(shell_particle_set)) THEN
         cpassert((SIZE(shell_particle_set, 1) == nshell))
         ALLOCATE (r(3, nshell), at_name(nshell), at_index(nshell))
         CALL cp_warn(__location__, &
                      "Overwriting shell coordinates. "// &
                      "Active coordinates READ from &SHELL_COORD section. ")
         CALL section_vals_list_get(shell_coord_section, "_DEFAULT_KEYWORD_", list=list)
         DO ishell = 1, nshell
            ! we use only the first default_string_length characters of each line
            is_ok = cp_sll_val_next(list, val)
            CALL val_get(val, c_val=line_att)
            start_c = 1
            DO wrd = 1, 5
               DO j = start_c, len(line_att)
                  IF (line_att(j:j) /= ' ') THEN
                     start_c = j
                     EXIT
                  END IF
               END DO
               end_c = len(line_att) + 1
               DO j = start_c, len(line_att)
                  IF (line_att(j:j) == ' ') THEN
                     end_c = j
                     EXIT
                  END IF
               END DO
               IF (wrd /= 5 .AND. end_c >= len(line_att) + 1) THEN
                  cpabort("incorrectly formatted line in coord section'"//line_att//"'")
               END IF
               IF (wrd == 1) THEN
                  at_name(ishell) = line_att(start_c:end_c - 1)
                  CALL uppercase(at_name(ishell))
               ELSE IF (wrd == 5) THEN
                  READ (line_att(start_c:end_c - 1), *) at_index(ishell)
               ELSE
                  READ (line_att(start_c:end_c - 1), *) r(wrd - 1, ishell)
               END IF
               start_c = end_c
            END DO
         END DO
 
         IF (PRESENT(core_particle_set)) THEN
            cpassert(ASSOCIATED(core_particle_set))
            core_coord_section => section_vals_get_subs_vals(subsys_section, "CORE_COORD")
            CALL section_vals_get(core_coord_section, explicit=explicit)
            IF (explicit) THEN
               CALL section_vals_val_get(core_coord_section, "UNIT", c_val=unit_str)
               CALL section_vals_val_get(core_coord_section, "SCALED", l_val=core_scaled_coordinates)
               unit_conv_core = cp_unit_to_cp2k(1.0_dp, trim(unit_str))
               CALL section_vals_val_get(core_coord_section, "_DEFAULT_KEYWORD_", &
                                         n_rep_val=nshell)
 
               cpassert((SIZE(core_particle_set, 1) == nshell))
               ALLOCATE (rc(3, nshell), at_name_c(nshell), at_index_c(nshell))
               CALL cp_warn(__location__, &
                            "Overwriting cores coordinates. "// &
                            "Active coordinates READ from &CORE_COORD section. ")
               CALL section_vals_list_get(core_coord_section, "_DEFAULT_KEYWORD_", list=list)
               DO ishell = 1, nshell
                  ! we use only the first default_string_length characters of each line
                  is_ok = cp_sll_val_next(list, val)
                  CALL val_get(val, c_val=line_att)
                  start_c = 1
                  DO wrd = 1, 5
                     DO j = start_c, len(line_att)
                        IF (line_att(j:j) /= ' ') THEN
                           start_c = j
                           EXIT
                        END IF
                     END DO
                     end_c = len(line_att) + 1
                     DO j = start_c, len(line_att)
                        IF (line_att(j:j) == ' ') THEN
                           end_c = j
                           EXIT
                        END IF
                     END DO
                     IF (wrd /= 5 .AND. end_c >= len(line_att) + 1) THEN
                        cpabort("incorrectly formatted line in coord section'"//line_att//"'")
                     END IF
                     IF (wrd == 1) THEN
                        at_name_c(ishell) = line_att(start_c:end_c - 1)
                        CALL uppercase(at_name_c(ishell))
                     ELSE IF (wrd == 5) THEN
                        READ (line_att(start_c:end_c - 1), *) at_index_c(ishell)
                     ELSE
                        READ (line_att(start_c:end_c - 1), *) rc(wrd - 1, ishell)
                     END IF
                     start_c = end_c
                  END DO
               END DO
               IF (my_save_mem) CALL section_vals_remove_values(core_coord_section)
            END IF ! explicit
         END IF ! core_particle_set
 
         !-----------------------------------------------------------------------------
         ! 3. Check corrispondence and convert coordinates into internal cp2k coordinates
         !-----------------------------------------------------------------------------
         DO ishell = 1, nshell
            atomic_kind => particle_set(at_index(ishell))%atomic_kind
            CALL get_atomic_kind(atomic_kind=atomic_kind, &
                                 name=name_kind, shell_active=is_shell, mass=mass_com, shell=shell)
            CALL uppercase(name_kind)
            IF ((trim(at_name(ishell)) == trim(name_kind)) .AND. is_shell) THEN
               sh_index = particle_set(at_index(ishell))%shell_index
               IF (shell_scaled_coordinates) THEN
                  CALL scaled_to_real(shell_particle_set(sh_index)%r(:), r(:, ishell), cell)
               ELSE
                  shell_particle_set(sh_index)%r(:) = r(:, ishell)*unit_conv_shell
                  CALL cell_transform_input_cartesian(cell, shell_particle_set(sh_index)%r(:))
               END IF
               shell_particle_set(sh_index)%atom_index = at_index(ishell)
 
               IF (PRESENT(core_particle_set) .AND. .NOT. explicit) THEN
                  core_particle_set(sh_index)%r(1) = (mass_com*particle_set(at_index(ishell))%r(1) - &
                                                      shell%mass_shell*shell_particle_set(sh_index)%r(1))/shell%mass_core
                  core_particle_set(sh_index)%r(2) = (mass_com*particle_set(at_index(ishell))%r(2) - &
                                                      shell%mass_shell*shell_particle_set(sh_index)%r(2))/shell%mass_core
                  core_particle_set(sh_index)%r(3) = (mass_com*particle_set(at_index(ishell))%r(3) - &
                                                      shell%mass_shell*shell_particle_set(sh_index)%r(3))/shell%mass_core
                  core_particle_set(sh_index)%atom_index = at_index(ishell)
                  rab = pbc(shell_particle_set(sh_index)%r, core_particle_set(sh_index)%r, cell)
               ELSE IF (explicit) THEN
                  IF (core_scaled_coordinates) THEN
                     CALL scaled_to_real(core_particle_set(sh_index)%r(:), rc(:, ishell), cell)
                  ELSE
                     core_particle_set(sh_index)%r(:) = rc(:, ishell)*unit_conv_core
                     CALL cell_transform_input_cartesian(cell, core_particle_set(sh_index)%r(:))
                  END IF
                  core_particle_set(sh_index)%atom_index = at_index_c(ishell)
                  rab = pbc(shell_particle_set(sh_index)%r, core_particle_set(sh_index)%r, cell)
                  cpassert(trim(at_name(ishell)) == trim(at_name_c(ishell)))
                  cpassert(at_index(ishell) == at_index_c(ishell))
               ELSE
                  rab = pbc(shell_particle_set(sh_index)%r, particle_set(at_index(ishell))%r, cell)
               END IF
 
               dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
               IF (shell%max_dist > 0.0_dp .AND. shell%max_dist < dab) THEN
                  IF (output_unit > 0) THEN
                     WRITE (output_unit, *) "WARNING : shell and core for atom ", at_index(ishell), " seem to be too distant."
                  END IF
               END IF
 
            ELSE
               cpabort("shell coordinate assigned to the wrong atom. check the shell indexes in the input")
            END IF
         END DO
         DEALLOCATE (r, at_index, at_name)
         DEALLOCATE (rc, at_index_c, at_name_c)
 
      END IF
 
      IF (my_save_mem) CALL section_vals_remove_values(shell_coord_section)
 
      CALL timestop(handle)
 
   END SUBROUTINE read_shell_coord_input
 
END MODULE atoms_input