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