atoms_input.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 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_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) &
                  cpabort("incorrectly formatted line in coord section'"//line_att//"'")
               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) &
                  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)//"> ")
               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) &
                     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)//"> ")
               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
         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) &
                  cpabort("incorrectly formatted line in coord section'"//line_att//"'")
               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) &
                        cpabort("incorrectly formatted line in coord section'"//line_att//"'")
                     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(r(:, ishell), shell_particle_set(sh_index)%r(:), cell)
               ELSE
                  shell_particle_set(sh_index)%r(:) = r(:, ishell)*unit_conv_shell
               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(rc(:, ishell), core_particle_set(sh_index)%r(:), cell)
                  ELSE
                     core_particle_set(sh_index)%r(:) = rc(:, ishell)*unit_conv_core
                  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