LCOV - code coverage report
Current view: top level - src - topology_xtl.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 77.7 % 157 122
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 1 1

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief  Handles XTL (Molecular Simulations, Inc (MSI)) files
      10              : !> \author Teodoro Laino [tlaino]
      11              : !> \date   05.2009
      12              : ! **************************************************************************************************
      13              : MODULE topology_xtl
      14              :    USE cell_methods,                    ONLY: cell_create,&
      15              :                                               set_cell_param,&
      16              :                                               write_cell
      17              :    USE cell_types,                      ONLY: cell_release,&
      18              :                                               cell_type,&
      19              :                                               pbc,&
      20              :                                               scaled_to_real
      21              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      22              :                                               cp_logger_type,&
      23              :                                               cp_to_string
      24              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      25              :                                               cp_print_key_unit_nr
      26              :    USE cp_parser_methods,               ONLY: parser_get_next_line,&
      27              :                                               parser_get_object,&
      28              :                                               parser_search_string
      29              :    USE cp_parser_types,                 ONLY: cp_parser_type,&
      30              :                                               parser_create,&
      31              :                                               parser_release
      32              :    USE cp_units,                        ONLY: cp_unit_to_cp2k
      33              :    USE input_section_types,             ONLY: section_get_rval,&
      34              :                                               section_vals_type
      35              :    USE kinds,                           ONLY: default_string_length,&
      36              :                                               dp
      37              :    USE memory_utilities,                ONLY: reallocate
      38              :    USE message_passing,                 ONLY: mp_para_env_type
      39              :    USE string_table,                    ONLY: id2str,&
      40              :                                               s2s,&
      41              :                                               str2id
      42              :    USE topology_types,                  ONLY: atom_info_type,&
      43              :                                               topology_parameters_type
      44              : #include "./base/base_uses.f90"
      45              : 
      46              :    IMPLICIT NONE
      47              : 
      48              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_xtl'
      49              : 
      50              :    PRIVATE
      51              :    PUBLIC :: read_coordinate_xtl
      52              : 
      53              : CONTAINS
      54              : 
      55              : ! **************************************************************************************************
      56              : !> \brief  Performs the real task of reading the proper information from the XTL
      57              : !>         file
      58              : !> \param topology ...
      59              : !> \param para_env ...
      60              : !> \param subsys_section ...
      61              : !> \date   05.2009
      62              : !> \par    Format Information implemented:
      63              : !>            TITLE
      64              : !>            DIMENSION
      65              : !>            CELL
      66              : !>            SYMMETRY
      67              : !>            SYM MAT
      68              : !>            ATOMS
      69              : !>            EOF
      70              : !>
      71              : !> \author Teodoro Laino [tlaino]
      72              : ! **************************************************************************************************
      73            6 :    SUBROUTINE read_coordinate_xtl(topology, para_env, subsys_section)
      74              :       TYPE(topology_parameters_type)                     :: topology
      75              :       TYPE(mp_para_env_type), POINTER                    :: para_env
      76              :       TYPE(section_vals_type), POINTER                   :: subsys_section
      77              : 
      78              :       CHARACTER(len=*), PARAMETER :: routineN = 'read_coordinate_xtl'
      79              :       INTEGER, PARAMETER                                 :: nblock = 1000
      80              :       REAL(KIND=dp), PARAMETER                           :: threshold = 1.0E-6_dp
      81              : 
      82              :       CHARACTER(LEN=default_string_length)               :: strtmp
      83              :       INTEGER                                            :: dimensions, handle, icol, ii, isym, iw, &
      84              :                                                             jj, natom, natom_orig, newsize
      85              :       INTEGER, DIMENSION(3)                              :: periodic
      86              :       LOGICAL                                            :: check, found, my_end
      87              :       REAL(KIND=dp)                                      :: pfactor, threshold2
      88              :       REAL(KIND=dp), DIMENSION(3)                        :: cell_angles, cell_lengths, r, r1, r2, s, &
      89              :                                                             transl_vec
      90              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: rot_mat
      91              :       TYPE(atom_info_type), POINTER                      :: atom_info
      92              :       TYPE(cell_type), POINTER                           :: cell
      93              :       TYPE(cp_logger_type), POINTER                      :: logger
      94              :       TYPE(cp_parser_type)                               :: parser
      95              : 
      96            2 :       NULLIFY (logger)
      97            4 :       logger => cp_get_default_logger()
      98              :       iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/XTL_INFO", &
      99            2 :                                 extension=".subsysLog")
     100            2 :       CALL timeset(routineN, handle)
     101              : 
     102            2 :       pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
     103              :       ! Element is assigned on the basis of the atm_name
     104            2 :       topology%aa_element = .TRUE.
     105              : 
     106            2 :       atom_info => topology%atom_info
     107            2 :       CALL reallocate(atom_info%id_molname, 1, nblock)
     108            2 :       CALL reallocate(atom_info%id_resname, 1, nblock)
     109            2 :       CALL reallocate(atom_info%resid, 1, nblock)
     110            2 :       CALL reallocate(atom_info%id_atmname, 1, nblock)
     111            2 :       CALL reallocate(atom_info%r, 1, 3, 1, nblock)
     112            2 :       CALL reallocate(atom_info%atm_mass, 1, nblock)
     113            2 :       CALL reallocate(atom_info%atm_charge, 1, nblock)
     114            2 :       CALL reallocate(atom_info%occup, 1, nblock)
     115            2 :       CALL reallocate(atom_info%beta, 1, nblock)
     116            2 :       CALL reallocate(atom_info%id_element, 1, nblock)
     117              : 
     118            2 :       IF (iw > 0) WRITE (iw, *) "    Reading in XTL file ", TRIM(topology%coord_file_name)
     119            2 :       CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
     120              : 
     121              :       ! Check for TITLE
     122              :       CALL parser_search_string(parser, "TITLE", ignore_case=.FALSE., found=found, &
     123            2 :                                 begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
     124            2 :       IF (found) THEN
     125            2 :          IF (iw > 0) WRITE (iw, '(/,A)') " XTL_INFO| TITLE :: "//TRIM(parser%input_line(parser%icol:))
     126              :       END IF
     127              : 
     128              :       ! Check for   _chemical_formula_sum
     129              :       CALL parser_search_string(parser, "DIMENSION", ignore_case=.FALSE., found=found, &
     130            2 :                                 begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
     131            2 :       IF (found) THEN
     132            2 :          IF (iw > 0) WRITE (iw, '(A)') " XTL_INFO| DIMENSION :: "//TRIM(parser%input_line(parser%icol:))
     133            2 :          CALL parser_get_object(parser, dimensions)
     134            2 :          IF (dimensions /= 3) THEN
     135            0 :             CPABORT("XTL file with working DIMENSION different from 3 cannot be parsed!")
     136              :          END IF
     137              :       ELSE
     138              :          ! Assuming by default we work in 3D-periodic systems
     139            0 :          dimensions = 3
     140              :       END IF
     141              : 
     142              :       ! Parsing cell infos
     143            8 :       periodic = 1
     144              :       ! Check for   _cell_length_a
     145              :       CALL parser_search_string(parser, "CELL", ignore_case=.FALSE., found=found, &
     146            2 :                                 begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
     147            2 :       IF (.NOT. found) &
     148            0 :          CPABORT("The field CELL was not found in XTL file! ")
     149            2 :       CALL parser_get_next_line(parser, 1)
     150              :       ! CELL LENGTH A
     151            2 :       CALL parser_get_object(parser, cell_lengths(1))
     152            2 :       cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
     153              :       ! CELL LENGTH B
     154            2 :       CALL parser_get_object(parser, cell_lengths(2))
     155            2 :       cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
     156              :       ! CELL LENGTH C
     157            2 :       CALL parser_get_object(parser, cell_lengths(3))
     158            2 :       cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
     159              : 
     160              :       ! CELL ANGLE ALPHA
     161            2 :       CALL parser_get_object(parser, cell_angles(1))
     162            2 :       cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
     163              :       ! CELL ANGLE BETA
     164            2 :       CALL parser_get_object(parser, cell_angles(2))
     165            2 :       cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
     166              :       ! CELL ANGLE GAMMA
     167            2 :       CALL parser_get_object(parser, cell_angles(3))
     168            2 :       cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")
     169              : 
     170              :       ! Create cell
     171            2 :       NULLIFY (cell)
     172            2 :       CALL cell_create(cell, tag="CELL_XTL")
     173              :       CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
     174            2 :                           do_init_cell=.TRUE.)
     175            2 :       CALL write_cell(cell, subsys_section)
     176              : 
     177              :       ! Parse atoms info and fractional coordinates
     178              :       ! Check for   _atom_site_label
     179              :       CALL parser_search_string(parser, "ATOMS", ignore_case=.FALSE., found=found, &
     180            2 :                                 begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
     181            2 :       IF (.NOT. found) &
     182            0 :          CPABORT("The field ATOMS was not found in XTL file! ")
     183            2 :       CALL parser_get_next_line(parser, 1)
     184              :       ! Paranoic syntax check.. if this fails one should improve the description of XTL files
     185            2 :       found = (INDEX(parser%input_line, "NAME       X          Y          Z") /= 0)
     186            2 :       IF (.NOT. found) &
     187            0 :          CPABORT("The field ATOMS in XTL file, is not followed by name and coordinates tags! ")
     188            2 :       CALL parser_get_next_line(parser, 1)
     189              :       ! Parse real info
     190            2 :       natom = 0
     191           62 :       DO WHILE (INDEX(parser%input_line, "EOF") == 0)
     192           60 :          natom = natom + 1
     193              :          ! Resize in case needed
     194           60 :          IF (natom > SIZE(atom_info%id_molname)) THEN
     195            0 :             newsize = INT(pfactor*natom)
     196            0 :             CALL reallocate(atom_info%id_molname, 1, newsize)
     197            0 :             CALL reallocate(atom_info%id_resname, 1, newsize)
     198            0 :             CALL reallocate(atom_info%resid, 1, newsize)
     199            0 :             CALL reallocate(atom_info%id_atmname, 1, newsize)
     200            0 :             CALL reallocate(atom_info%r, 1, 3, 1, newsize)
     201            0 :             CALL reallocate(atom_info%atm_mass, 1, newsize)
     202            0 :             CALL reallocate(atom_info%atm_charge, 1, newsize)
     203            0 :             CALL reallocate(atom_info%occup, 1, newsize)
     204            0 :             CALL reallocate(atom_info%beta, 1, newsize)
     205            0 :             CALL reallocate(atom_info%id_element, 1, newsize)
     206              :          END IF
     207              :          ! NAME
     208           60 :          CALL parser_get_object(parser, strtmp)
     209           60 :          atom_info%id_atmname(natom) = str2id(strtmp)
     210           60 :          atom_info%id_molname(natom) = str2id(s2s("MOL"//TRIM(ADJUSTL(cp_to_string(natom)))))
     211           60 :          atom_info%id_resname(natom) = atom_info%id_molname(natom)
     212           60 :          atom_info%resid(natom) = 1
     213           60 :          atom_info%id_element(natom) = atom_info%id_atmname(natom)
     214              :          ! X
     215           60 :          CALL parser_get_object(parser, atom_info%r(1, natom))
     216              :          ! Y
     217           60 :          CALL parser_get_object(parser, atom_info%r(2, natom))
     218              :          ! Z
     219           60 :          CALL parser_get_object(parser, atom_info%r(3, natom))
     220          240 :          s = atom_info%r(1:3, natom)
     221           60 :          CALL scaled_to_real(atom_info%r(1:3, natom), s, cell)
     222           60 :          CALL parser_get_next_line(parser, 1, at_end=my_end)
     223           62 :          IF (my_end) EXIT
     224              :       END DO
     225              :       !
     226            2 :       threshold2 = threshold*threshold
     227              :       ! Preliminary check: check if atoms provided are really unique.. this is a paranoic
     228              :       ! check since they should be REALLY unique.. anyway..
     229           62 :       DO ii = 1, natom
     230          240 :          r1 = atom_info%r(1:3, ii)
     231          932 :          DO jj = ii + 1, natom
     232         3480 :             r2 = atom_info%r(1:3, jj)
     233         3480 :             r = pbc(r1 - r2, cell)
     234              :             ! SQRT(DOT_PRODUCT(r, r)) >= threshold
     235         3480 :             check = (DOT_PRODUCT(r, r) >= threshold2)
     236          930 :             CPASSERT(check)
     237              :          END DO
     238              :       END DO
     239              :       ! Parse Symmetry Group and generation elements..
     240              :       ! Check for SYMMETRY
     241              :       CALL parser_search_string(parser, "SYMMETRY", ignore_case=.FALSE., found=found, &
     242            2 :                                 begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
     243            2 :       IF (found) THEN
     244            2 :          IF (iw > 0) WRITE (iw, '(A)') " XTL_INFO| Symmetry Infos :: "//TRIM(parser%input_line(parser%icol:))
     245              :       END IF
     246              : 
     247              :       ! Check for SYM MAT
     248              :       CALL parser_search_string(parser, "SYM MAT", ignore_case=.FALSE., found=found, &
     249            2 :                                 begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
     250            2 :       CPWARN_IF(.NOT. found, "The field SYM MAT was not found in XTL file! ")
     251            2 :       IF (iw > 0) WRITE (iw, '(A,I0)') " XTL_INFO| Number of atoms before applying symmetry operations :: ", natom
     252           32 :       IF (iw > 0) WRITE (iw, '(A10,1X,3F12.6)') (TRIM(id2str(atom_info%id_atmname(ii))), atom_info%r(1:3, ii), ii=1, natom)
     253            2 :       IF (found) THEN
     254              :          ! Apply symmetry elements and generate the whole set of atoms in the unit cell
     255            2 :          isym = 0
     256            2 :          natom_orig = natom
     257            4 :          DO WHILE (found)
     258            2 :             isym = isym + 1
     259            2 :             icol = INDEX(parser%input_line, "SYM MAT") + 8
     260            8 :             READ (parser%input_line(icol:), *) ((rot_mat(ii, jj), jj=1, 3), ii=1, 3), transl_vec(1:3)
     261           62 :             Loop_over_unique_atoms: DO ii = 1, natom_orig
     262              :                ! Rotate and apply translation
     263          960 :                r1 = MATMUL(rot_mat, atom_info%r(1:3, ii)) + transl_vec
     264              :                ! Verify if this atom is really unique..
     265           60 :                check = .TRUE.
     266          930 :                DO jj = 1, natom
     267         3720 :                   r2 = atom_info%r(1:3, jj)
     268         3720 :                   r = pbc(r1 - r2, cell)
     269              :                   ! SQRT(DOT_PRODUCT(r, r)) <= threshold
     270         3720 :                   IF (DOT_PRODUCT(r, r) <= threshold2) THEN
     271              :                      check = .FALSE.
     272              :                      EXIT
     273              :                   END IF
     274              :                END DO
     275              :                ! If the atom generated is unique let's add to the atom set..
     276           62 :                IF (check) THEN
     277            0 :                   natom = natom + 1
     278              :                   ! Resize in case needed
     279            0 :                   IF (natom > SIZE(atom_info%id_molname)) THEN
     280            0 :                      newsize = INT(pfactor*natom)
     281            0 :                      CALL reallocate(atom_info%id_molname, 1, newsize)
     282            0 :                      CALL reallocate(atom_info%id_resname, 1, newsize)
     283            0 :                      CALL reallocate(atom_info%resid, 1, newsize)
     284            0 :                      CALL reallocate(atom_info%id_atmname, 1, newsize)
     285            0 :                      CALL reallocate(atom_info%r, 1, 3, 1, newsize)
     286            0 :                      CALL reallocate(atom_info%atm_mass, 1, newsize)
     287            0 :                      CALL reallocate(atom_info%atm_charge, 1, newsize)
     288            0 :                      CALL reallocate(atom_info%occup, 1, newsize)
     289            0 :                      CALL reallocate(atom_info%beta, 1, newsize)
     290            0 :                      CALL reallocate(atom_info%id_element, 1, newsize)
     291              :                   END IF
     292            0 :                   atom_info%id_atmname(natom) = atom_info%id_atmname(ii)
     293            0 :                   atom_info%id_molname(natom) = atom_info%id_molname(ii)
     294            0 :                   atom_info%id_resname(natom) = atom_info%id_resname(ii)
     295            0 :                   atom_info%resid(natom) = atom_info%resid(ii)
     296            0 :                   atom_info%id_element(natom) = atom_info%id_element(ii)
     297            0 :                   atom_info%r(1:3, natom) = r1
     298              :                END IF
     299              :             END DO Loop_over_unique_atoms
     300              :             CALL parser_search_string(parser, "SYM MAT", ignore_case=.FALSE., found=found, &
     301            4 :                                       begin_line=.FALSE., search_from_begin_of_file=.FALSE.)
     302              :          END DO
     303              :       END IF
     304            2 :       IF (iw > 0) WRITE (iw, '(A,I0)') " XTL_INFO| Number of symmetry operations :: ", isym
     305            2 :       IF (iw > 0) WRITE (iw, '(A,I0)') " XTL_INFO| Number of total atoms :: ", natom
     306           32 :       IF (iw > 0) WRITE (iw, '(A10,1X,3F12.6)') (TRIM(id2str(atom_info%id_atmname(ii))), atom_info%r(1:3, ii), ii=1, natom)
     307              : 
     308              :       ! Releasse local cell type and parser
     309            2 :       CALL cell_release(cell)
     310            2 :       CALL parser_release(parser)
     311              : 
     312              :       ! Reallocate all structures with the exact NATOM size
     313            2 :       CALL reallocate(atom_info%id_molname, 1, natom)
     314            2 :       CALL reallocate(atom_info%id_resname, 1, natom)
     315            2 :       CALL reallocate(atom_info%resid, 1, natom)
     316            2 :       CALL reallocate(atom_info%id_atmname, 1, natom)
     317            2 :       CALL reallocate(atom_info%r, 1, 3, 1, natom)
     318            2 :       CALL reallocate(atom_info%atm_mass, 1, natom)
     319            2 :       CALL reallocate(atom_info%atm_charge, 1, natom)
     320            2 :       CALL reallocate(atom_info%occup, 1, natom)
     321            2 :       CALL reallocate(atom_info%beta, 1, natom)
     322            2 :       CALL reallocate(atom_info%id_element, 1, natom)
     323              : 
     324            2 :       topology%natoms = natom
     325            2 :       topology%molname_generated = .TRUE.
     326              :       CALL cp_print_key_finished_output(iw, logger, subsys_section, &
     327            2 :                                         "PRINT%TOPOLOGY_INFO/XTL_INFO")
     328            2 :       CALL timestop(handle)
     329            6 :    END SUBROUTINE read_coordinate_xtl
     330              : 
     331              : END MODULE topology_xtl
        

Generated by: LCOV version 2.0-1