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 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 10274 : 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 10274 : 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 10274 : CALL timeset(routineN, handle)
58 :
59 10274 : NULLIFY (multiple_unit_cell, iwork, cell)
60 :
61 : ! Store original number of atoms for the molecule generation in any case
62 10274 : topology%natom_muc = topology%natoms
63 :
64 : CALL section_vals_val_get(subsys_section, "TOPOLOGY%MULTIPLE_UNIT_CELL", &
65 10274 : i_vals=multiple_unit_cell)
66 :
67 : ! Fail is one of the value is set to zero..
68 41096 : IF (ANY(multiple_unit_cell <= 0)) &
69 : CALL cp_abort(__LOCATION__, "SUBSYS%TOPOLOGY%MULTIPLE_UNIT_CELL accepts "// &
70 0 : "only integer values greater than zero.")
71 :
72 40662 : IF (ANY(multiple_unit_cell /= 1)) THEN
73 :
74 : ! Check that the setup between CELL and TOPOLOGY is the same
75 : CALL section_vals_val_get(subsys_section, "CELL%MULTIPLE_UNIT_CELL", &
76 148 : i_vals=iwork)
77 592 : IF (ANY(iwork /= multiple_unit_cell)) &
78 : CALL cp_abort(__LOCATION__, "The input parameters for "// &
79 : "SUBSYS%TOPOLOGY%MULTIPLE_UNIT_CELL and "// &
80 0 : "SUBSYS%CELL%MULTIPLE_UNIT_CELL have to agree.")
81 :
82 148 : cell => topology%cell_muc
83 592 : natoms = topology%natoms*PRODUCT(multiple_unit_cell)
84 :
85 : ! Check, if velocities are provided, that they are consistent in number with the atoms...
86 148 : work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
87 148 : CALL section_vals_get(work_section, explicit=explicit)
88 148 : IF (explicit) THEN
89 0 : CALL section_vals_val_get(work_section, '_DEFAULT_KEYWORD_', n_rep_val=nrep)
90 0 : check = nrep == natoms
91 0 : IF (.NOT. check) &
92 : CALL cp_abort(__LOCATION__, "The number of available entries in the "// &
93 0 : "VELOCITY section is not compatible with the number of atoms.")
94 : END IF
95 :
96 148 : CALL reallocate(topology%atom_info%id_molname, 1, natoms)
97 148 : CALL reallocate(topology%atom_info%id_resname, 1, natoms)
98 148 : CALL reallocate(topology%atom_info%resid, 1, natoms)
99 148 : CALL reallocate(topology%atom_info%id_atmname, 1, natoms)
100 148 : CALL reallocate(topology%atom_info%r, 1, 3, 1, natoms)
101 148 : CALL reallocate(topology%atom_info%atm_mass, 1, natoms)
102 148 : CALL reallocate(topology%atom_info%atm_charge, 1, natoms)
103 148 : CALL reallocate(topology%atom_info%occup, 1, natoms)
104 148 : CALL reallocate(topology%atom_info%beta, 1, natoms)
105 148 : CALL reallocate(topology%atom_info%id_element, 1, natoms)
106 :
107 148 : ind = 0
108 556 : DO k = 1, multiple_unit_cell(3)
109 1632 : trsl_k = cell%hmat(:, 3)*REAL(k - 1, KIND=dp)
110 1870 : DO j = 1, multiple_unit_cell(2)
111 5256 : trsl_j = cell%hmat(:, 2)*REAL(j - 1, KIND=dp)
112 6698 : DO i = 1, multiple_unit_cell(1)
113 19904 : trsl_i = cell%hmat(:, 1)*REAL(i - 1, KIND=dp)
114 19904 : trsl = trsl_i + trsl_j + trsl_k
115 4976 : ind = ind + 1
116 4976 : IF (ind == 1) CYCLE
117 : ! Loop over all atoms
118 4828 : n = (ind - 1)*topology%natoms
119 37284 : DO m = 1, topology%natoms
120 31142 : topology%atom_info%id_atmname(n + m) = topology%atom_info%id_atmname(m)
121 31142 : topology%atom_info%r(1, n + m) = topology%atom_info%r(1, m) + trsl(1)
122 31142 : topology%atom_info%r(2, n + m) = topology%atom_info%r(2, m) + trsl(2)
123 31142 : topology%atom_info%r(3, n + m) = topology%atom_info%r(3, m) + trsl(3)
124 31142 : topology%atom_info%id_molname(n + m) = topology%atom_info%id_molname(m)
125 31142 : topology%atom_info%id_resname(n + m) = topology%atom_info%id_resname(m)
126 31142 : topology%atom_info%resid(n + m) = topology%atom_info%resid(m)
127 31142 : topology%atom_info%id_element(n + m) = topology%atom_info%id_element(m)
128 31142 : topology%atom_info%atm_mass(n + m) = topology%atom_info%atm_mass(m)
129 36118 : topology%atom_info%atm_charge(n + m) = topology%atom_info%atm_charge(m)
130 : END DO
131 : END DO
132 : END DO
133 : END DO
134 : ! Store the new total number of atoms
135 148 : topology%natoms = natoms
136 :
137 : ! Deallocate the coordinate section (will be rebuilt later with the whole atomic set)
138 148 : work_section => section_vals_get_subs_vals(subsys_section, "COORD")
139 148 : CALL section_vals_get(work_section, explicit=explicit)
140 148 : IF (explicit) THEN
141 140 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
142 140 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
143 : END IF
144 148 : CALL section_vals_remove_values(work_section)
145 148 : IF (explicit) THEN
146 140 : CALL section_vals_val_set(work_section, "UNIT", c_val=unit_str)
147 140 : CALL section_vals_val_set(work_section, "SCALED", l_val=scale)
148 : END IF
149 : END IF
150 :
151 10274 : CALL timestop(handle)
152 :
153 10274 : END SUBROUTINE topology_muc
154 :
155 : END MODULE topology_multiple_unit_cell
|
