Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Handles the multiple unit cell option regarding atomic coordinates
      10             : !> \author Teodoro Laino [tlaino] - 05.2009
      11             : ! **************************************************************************************************
      12             : MODULE topology_multiple_unit_cell
      13             :    USE cell_types,                      ONLY: cell_type
      14             :    USE input_section_types,             ONLY: section_vals_get,&
      15             :                                               section_vals_get_subs_vals,&
      16             :                                               section_vals_remove_values,&
      17             :                                               section_vals_type,&
      18             :                                               section_vals_val_get,&
      19             :                                               section_vals_val_set
      20             :    USE kinds,                           ONLY: default_string_length,&
      21             :                                               dp
      22             :    USE memory_utilities,                ONLY: reallocate
      23             :    USE topology_types,                  ONLY: topology_parameters_type
      24             : #include "./base/base_uses.f90"
      25             : 
      26             :    IMPLICIT NONE
      27             : 
      28             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_multiple_unit_cell'
      29             : 
      30             :    PRIVATE
      31             : 
      32             : ! *** Public parameters ***
      33             :    PUBLIC :: topology_muc
      34             : 
      35             : CONTAINS
      36             : 
      37             : ! **************************************************************************************************
      38             : !> \brief Handles the multiple_unit_cell for the atomic coordinates..
      39             : !> \param topology ...
      40             : !> \param subsys_section ...
      41             : !> \author Teodoro Laino [tlaino] - 05.2009
      42             : ! **************************************************************************************************
      43        9512 :    SUBROUTINE topology_muc(topology, subsys_section)
      44             :       TYPE(topology_parameters_type), INTENT(INOUT)      :: topology
      45             :       TYPE(section_vals_type), POINTER                   :: subsys_section
      46             : 
      47             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'topology_muc'
      48             : 
      49             :       CHARACTER(LEN=default_string_length)               :: unit_str
      50             :       INTEGER                                            :: handle, i, ind, j, k, m, n, natoms, nrep
      51        9512 :       INTEGER, DIMENSION(:), POINTER                     :: iwork, multiple_unit_cell
      52             :       LOGICAL                                            :: check, explicit, scale
      53             :       REAL(KIND=dp), DIMENSION(3)                        :: trsl, trsl_i, trsl_j, trsl_k
      54             :       TYPE(cell_type), POINTER                           :: cell
      55             :       TYPE(section_vals_type), POINTER                   :: work_section
      56             : 
      57        9512 :       CALL timeset(routineN, handle)
      58        9512 :       NULLIFY (multiple_unit_cell, iwork, cell)
      59             :       CALL section_vals_val_get(subsys_section, "TOPOLOGY%MULTIPLE_UNIT_CELL", &
      60        9512 :                                 i_vals=multiple_unit_cell)
      61             :       ! Fail is one of the value is set to zero..
      62       38048 :       IF (ANY(multiple_unit_cell <= 0)) &
      63             :          CALL cp_abort(__LOCATION__, "SUBSYS%TOPOLOGY%MULTIPLE_UNIT_CELL accepts "// &
      64           0 :                        "only integer values larger than 0! A value of 0 or negative is meaningless!")
      65       37614 :       IF (ANY(multiple_unit_cell /= 1)) THEN
      66             :          ! Check that the setup between CELL and TOPOLOGY is the same..
      67             :          CALL section_vals_val_get(subsys_section, "CELL%MULTIPLE_UNIT_CELL", &
      68         148 :                                    i_vals=iwork)
      69         592 :          IF (ANY(iwork /= multiple_unit_cell)) &
      70             :             CALL cp_abort(__LOCATION__, "SUBSYS%TOPOLOGY%MULTIPLE_UNIT_CELL and "// &
      71             :                           "SUBSYS%CELL%MULTIPLE_UNIT_CELL have been "// &
      72           0 :                           "setup to two different values!! Correct this error!")
      73         148 :          cell => topology%cell_muc
      74         592 :          natoms = topology%natoms*PRODUCT(multiple_unit_cell)
      75             : 
      76             :          ! Check, if velocities are provided, that they are consistent in number with the atoms...
      77         148 :          work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
      78         148 :          CALL section_vals_get(work_section, explicit=explicit)
      79         148 :          IF (explicit) THEN
      80           0 :             CALL section_vals_val_get(work_section, '_DEFAULT_KEYWORD_', n_rep_val=nrep)
      81           0 :             check = nrep == natoms
      82           0 :             IF (.NOT. check) &
      83             :                CALL cp_abort(__LOCATION__, &
      84           0 :                              "Number of available entries in VELOCITY section is not compatible with the number of atoms!")
      85             :          END IF
      86             : 
      87         148 :          CALL reallocate(topology%atom_info%id_molname, 1, natoms)
      88         148 :          CALL reallocate(topology%atom_info%id_resname, 1, natoms)
      89         148 :          CALL reallocate(topology%atom_info%resid, 1, natoms)
      90         148 :          CALL reallocate(topology%atom_info%id_atmname, 1, natoms)
      91         148 :          CALL reallocate(topology%atom_info%r, 1, 3, 1, natoms)
      92         148 :          CALL reallocate(topology%atom_info%atm_mass, 1, natoms)
      93         148 :          CALL reallocate(topology%atom_info%atm_charge, 1, natoms)
      94         148 :          CALL reallocate(topology%atom_info%occup, 1, natoms)
      95         148 :          CALL reallocate(topology%atom_info%beta, 1, natoms)
      96         148 :          CALL reallocate(topology%atom_info%id_element, 1, natoms)
      97         148 :          ind = 0
      98         556 :          DO k = 1, multiple_unit_cell(3)
      99        1632 :             trsl_k = cell%hmat(:, 3)*REAL(k - 1, KIND=dp)
     100        1870 :             DO j = 1, multiple_unit_cell(2)
     101        5256 :                trsl_j = cell%hmat(:, 2)*REAL(j - 1, KIND=dp)
     102        6698 :                DO i = 1, multiple_unit_cell(1)
     103       19904 :                   trsl_i = cell%hmat(:, 1)*REAL(i - 1, KIND=dp)
     104       19904 :                   trsl = trsl_i + trsl_j + trsl_k
     105        4976 :                   ind = ind + 1
     106        4976 :                   IF (ind == 1) CYCLE
     107             : 
     108             :                   ! loop over atoms
     109        4828 :                   n = (ind - 1)*topology%natoms
     110       37284 :                   DO m = 1, topology%natoms
     111       31142 :                      topology%atom_info%id_atmname(n + m) = topology%atom_info%id_atmname(m)
     112       31142 :                      topology%atom_info%r(1, n + m) = topology%atom_info%r(1, m) + trsl(1)
     113       31142 :                      topology%atom_info%r(2, n + m) = topology%atom_info%r(2, m) + trsl(2)
     114       31142 :                      topology%atom_info%r(3, n + m) = topology%atom_info%r(3, m) + trsl(3)
     115       31142 :                      topology%atom_info%id_molname(n + m) = topology%atom_info%id_molname(m)
     116       31142 :                      topology%atom_info%id_resname(n + m) = topology%atom_info%id_resname(m)
     117       31142 :                      topology%atom_info%resid(n + m) = topology%atom_info%resid(m)
     118       31142 :                      topology%atom_info%id_element(n + m) = topology%atom_info%id_element(m)
     119       31142 :                      topology%atom_info%atm_mass(n + m) = topology%atom_info%atm_mass(m)
     120       36118 :                      topology%atom_info%atm_charge(n + m) = topology%atom_info%atm_charge(m)
     121             :                   END DO
     122             :                END DO
     123             :             END DO
     124             :          END DO
     125         148 :          topology%natoms = natoms
     126             : 
     127             :          ! Deallocate the coordinate section (will be rebuilt later with the whole atomic set)
     128         148 :          work_section => section_vals_get_subs_vals(subsys_section, "COORD")
     129         148 :          CALL section_vals_get(work_section, explicit=explicit)
     130         148 :          IF (explicit) THEN
     131         140 :             CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
     132         140 :             CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
     133             :          END IF
     134         148 :          CALL section_vals_remove_values(work_section)
     135         148 :          IF (explicit) THEN
     136         140 :             CALL section_vals_val_set(work_section, "UNIT", c_val=unit_str)
     137         140 :             CALL section_vals_val_set(work_section, "SCALED", l_val=scale)
     138             :          END IF
     139             :       END IF
     140        9512 :       CALL timestop(handle)
     141        9512 :    END SUBROUTINE topology_muc
     142             : 
     143             : END MODULE topology_multiple_unit_cell