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 Initialize a small environment for a particular calculation
10 : !> \par History
11 : !> 5.2004 created [fawzi]
12 : !> 9.2007 cleaned [tlaino] - University of Zurich
13 : !> \author Teodoro Laino
14 : ! **************************************************************************************************
15 : MODULE cp_subsys_methods
16 : USE atomic_kind_list_types, ONLY: atomic_kind_list_create,&
17 : atomic_kind_list_release,&
18 : atomic_kind_list_type
19 : USE atomic_kind_types, ONLY: atomic_kind_type
20 : USE atprop_types, ONLY: atprop_create
21 : USE cell_methods, ONLY: write_cell
22 : USE cell_types, ONLY: cell_retain,&
23 : cell_type
24 : USE colvar_methods, ONLY: colvar_read
25 : USE cp_result_types, ONLY: cp_result_create
26 : USE cp_subsys_types, ONLY: cp_subsys_get,&
27 : cp_subsys_set,&
28 : cp_subsys_type
29 : USE exclusion_types, ONLY: exclusion_type
30 : USE input_constants, ONLY: do_conn_off,&
31 : do_stress_analytical,&
32 : do_stress_diagonal_anal,&
33 : do_stress_diagonal_numer,&
34 : do_stress_none,&
35 : do_stress_numerical
36 : USE input_section_types, ONLY: section_vals_get,&
37 : section_vals_get_subs_vals,&
38 : section_vals_type,&
39 : section_vals_val_get
40 : USE kinds, ONLY: default_string_length,&
41 : dp
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE molecule_kind_list_types, ONLY: molecule_kind_list_create,&
44 : molecule_kind_list_release,&
45 : molecule_kind_list_type
46 : USE molecule_kind_types, ONLY: molecule_kind_type
47 : USE molecule_list_types, ONLY: molecule_list_create,&
48 : molecule_list_release,&
49 : molecule_list_type
50 : USE molecule_types, ONLY: molecule_type
51 : USE particle_list_types, ONLY: particle_list_create,&
52 : particle_list_release,&
53 : particle_list_type
54 : USE particle_types, ONLY: particle_type
55 : USE qmmm_types_low, ONLY: qmmm_env_mm_type
56 : USE string_table, ONLY: id2str,&
57 : s2s,&
58 : str2id
59 : USE topology, ONLY: connectivity_control,&
60 : topology_control
61 : USE topology_connectivity_util, ONLY: topology_connectivity_pack
62 : USE topology_coordinate_util, ONLY: topology_coordinate_pack
63 : USE topology_types, ONLY: deallocate_topology,&
64 : init_topology,&
65 : topology_parameters_type
66 : USE topology_util, ONLY: check_subsys_element
67 : USE virial_types, ONLY: virial_set
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 : PRIVATE
72 :
73 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
74 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_subsys_methods'
75 :
76 : PUBLIC :: create_small_subsys, cp_subsys_create
77 :
78 : CONTAINS
79 :
80 : ! **************************************************************************************************
81 : !> \brief Creates allocates and fills subsys from given input.
82 : !> \param subsys ...
83 : !> \param para_env ...
84 : !> \param root_section ...
85 : !> \param force_env_section ...
86 : !> \param subsys_section ...
87 : !> \param use_motion_section ...
88 : !> \param qmmm ...
89 : !> \param qmmm_env ...
90 : !> \param exclusions ...
91 : !> \param elkind ...
92 : !> \author Ole Schuett
93 : ! **************************************************************************************************
94 30087 : SUBROUTINE cp_subsys_create(subsys, para_env, &
95 : root_section, force_env_section, subsys_section, &
96 : use_motion_section, qmmm, qmmm_env, exclusions, elkind)
97 : TYPE(cp_subsys_type), POINTER :: subsys
98 : TYPE(mp_para_env_type), POINTER :: para_env
99 : TYPE(section_vals_type), POINTER :: root_section
100 : TYPE(section_vals_type), OPTIONAL, POINTER :: force_env_section, subsys_section
101 : LOGICAL, INTENT(IN), OPTIONAL :: use_motion_section
102 : LOGICAL, OPTIONAL :: qmmm
103 : TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env
104 : TYPE(exclusion_type), DIMENSION(:), OPTIONAL, &
105 : POINTER :: exclusions
106 : LOGICAL, INTENT(IN), OPTIONAL :: elkind
107 :
108 : INTEGER :: stress_tensor
109 10029 : INTEGER, DIMENSION(:), POINTER :: seed_vals
110 : LOGICAL :: atomic_energy, my_use_motion_section, &
111 : pv_availability, pv_diagonal, &
112 : pv_numerical
113 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
114 10029 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
115 : TYPE(molecule_kind_list_type), POINTER :: mol_kinds
116 10029 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
117 : TYPE(molecule_list_type), POINTER :: mols
118 10029 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
119 : TYPE(particle_list_type), POINTER :: particles
120 10029 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
121 : TYPE(section_vals_type), POINTER :: colvar_section, my_force_env_section, &
122 : my_subsys_section
123 :
124 0 : CPASSERT(.NOT. ASSOCIATED(subsys))
125 100290 : ALLOCATE (subsys)
126 :
127 10029 : CALL para_env%retain()
128 10029 : subsys%para_env => para_env
129 :
130 10029 : my_use_motion_section = .FALSE.
131 10029 : IF (PRESENT(use_motion_section)) &
132 10027 : my_use_motion_section = use_motion_section
133 :
134 10029 : my_force_env_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
135 10029 : IF (PRESENT(force_env_section)) &
136 10027 : my_force_env_section => force_env_section
137 :
138 10029 : my_subsys_section => section_vals_get_subs_vals(my_force_env_section, "SUBSYS")
139 10029 : IF (PRESENT(subsys_section)) &
140 9873 : my_subsys_section => subsys_section
141 :
142 10029 : CALL section_vals_val_get(my_subsys_section, "SEED", i_vals=seed_vals)
143 10029 : IF (SIZE(seed_vals) == 1) THEN
144 90225 : subsys%seed(:, :) = REAL(seed_vals(1), KIND=dp)
145 4 : ELSE IF (SIZE(seed_vals) == 6) THEN
146 60 : subsys%seed(1:3, 1:2) = RESHAPE(REAL(seed_vals(:), KIND=dp), [3, 2])
147 : ELSE
148 0 : CPABORT("Supply exactly 1 or 6 arguments for SEED in &SUBSYS only!")
149 : END IF
150 :
151 10029 : colvar_section => section_vals_get_subs_vals(my_subsys_section, "COLVAR")
152 :
153 10029 : CALL cp_subsys_read_colvar(subsys, colvar_section)
154 :
155 : ! *** Read the particle coordinates and allocate the atomic kind, ***
156 : ! *** the molecule kind, and the molecule data structures ***
157 : CALL topology_control(atomic_kind_set, particle_set, molecule_kind_set, molecule_set, &
158 : subsys%colvar_p, subsys%gci, root_section, para_env, &
159 : force_env_section=my_force_env_section, &
160 : subsys_section=my_subsys_section, use_motion_section=my_use_motion_section, &
161 : qmmm=qmmm, qmmm_env=qmmm_env, exclusions=exclusions, elkind=elkind, &
162 10029 : subsys=subsys)
163 :
164 10029 : CALL particle_list_create(particles, els_ptr=particle_set)
165 10029 : CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
166 10029 : CALL molecule_list_create(mols, els_ptr=molecule_set)
167 10029 : CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
168 :
169 : CALL cp_subsys_set(subsys, particles=particles, atomic_kinds=atomic_kinds, &
170 10029 : molecules=mols, molecule_kinds=mol_kinds)
171 :
172 10029 : CALL particle_list_release(particles)
173 10029 : CALL atomic_kind_list_release(atomic_kinds)
174 10029 : CALL molecule_list_release(mols)
175 10029 : CALL molecule_kind_list_release(mol_kinds)
176 :
177 : ! Should we compute the virial?
178 10029 : CALL section_vals_val_get(my_force_env_section, "STRESS_TENSOR", i_val=stress_tensor)
179 9197 : SELECT CASE (stress_tensor)
180 : CASE (do_stress_none)
181 9197 : pv_availability = .FALSE.
182 9197 : pv_numerical = .FALSE.
183 9197 : pv_diagonal = .FALSE.
184 : CASE (do_stress_analytical)
185 820 : pv_availability = .TRUE.
186 820 : pv_numerical = .FALSE.
187 820 : pv_diagonal = .FALSE.
188 : CASE (do_stress_numerical)
189 2 : pv_availability = .TRUE.
190 2 : pv_numerical = .TRUE.
191 2 : pv_diagonal = .FALSE.
192 : CASE (do_stress_diagonal_anal)
193 0 : pv_availability = .TRUE.
194 0 : pv_numerical = .FALSE.
195 0 : pv_diagonal = .TRUE.
196 : CASE (do_stress_diagonal_numer)
197 10 : pv_availability = .TRUE.
198 10 : pv_numerical = .TRUE.
199 10029 : pv_diagonal = .TRUE.
200 : END SELECT
201 :
202 2487192 : ALLOCATE (subsys%virial)
203 : CALL virial_set(virial=subsys%virial, &
204 : pv_availability=pv_availability, &
205 : pv_numer=pv_numerical, &
206 10029 : pv_diagonal=pv_diagonal)
207 :
208 : ! Should we compute atomic properties?
209 10029 : CALL atprop_create(subsys%atprop)
210 10029 : CALL section_vals_val_get(my_force_env_section, "PROPERTIES%ATOMIC%ENERGY", l_val=atomic_energy)
211 10029 : subsys%atprop%energy = atomic_energy
212 :
213 10029 : CALL cp_result_create(subsys%results)
214 10029 : END SUBROUTINE cp_subsys_create
215 :
216 : ! **************************************************************************************************
217 : !> \brief reads the colvar section of the colvar
218 : !> \param subsys ...
219 : !> \param colvar_section ...
220 : !> \par History
221 : !> 2006.01 Joost VandeVondele
222 : ! **************************************************************************************************
223 10029 : SUBROUTINE cp_subsys_read_colvar(subsys, colvar_section)
224 : TYPE(cp_subsys_type), POINTER :: subsys
225 : TYPE(section_vals_type), POINTER :: colvar_section
226 :
227 : INTEGER :: ig, ncol
228 :
229 10029 : CALL section_vals_get(colvar_section, n_repetition=ncol)
230 20800 : ALLOCATE (subsys%colvar_p(ncol))
231 10495 : DO ig = 1, ncol
232 466 : NULLIFY (subsys%colvar_p(ig)%colvar)
233 10495 : CALL colvar_read(subsys%colvar_p(ig)%colvar, ig, colvar_section, subsys%para_env)
234 : END DO
235 10029 : END SUBROUTINE cp_subsys_read_colvar
236 :
237 : ! **************************************************************************************************
238 : !> \brief updates the molecule information of the given subsys
239 : !> \param small_subsys the subsys to create
240 : !> \param big_subsys the superset of small_subsys
241 : !> \param small_cell the cell of small_subsys
242 : !> \param small_para_env the parallel environment for the new (small)
243 : !> subsys
244 : !> \param sub_atom_index indexes of the atoms that should be in small_subsys
245 : !> \param sub_atom_kind_name ...
246 : !> \param para_env ...
247 : !> \param force_env_section ...
248 : !> \param subsys_section ...
249 : !> \param ignore_outside_box ...
250 : !> \par History
251 : !> 05.2004 created [fawzi]
252 : !> \author Fawzi Mohamed, Teodoro Laino
253 : !> \note
254 : !> not really ready to be used with different para_envs for the small
255 : !> and big part
256 : !> qmmm_env_create() should be the only instance where this subroutine is called
257 : !> with small_cell distinct from big_subsys%cell; otherwise it could be possible
258 : !> to drop the small_cell as input parameter entirely.
259 : ! **************************************************************************************************
260 538 : SUBROUTINE create_small_subsys(small_subsys, big_subsys, small_cell, &
261 538 : small_para_env, sub_atom_index, sub_atom_kind_name, &
262 : para_env, force_env_section, subsys_section, ignore_outside_box)
263 :
264 : TYPE(cp_subsys_type), POINTER :: small_subsys, big_subsys
265 : TYPE(cell_type), POINTER :: small_cell
266 : TYPE(mp_para_env_type), POINTER :: small_para_env
267 : INTEGER, DIMENSION(:), INTENT(in) :: sub_atom_index
268 : CHARACTER(len=default_string_length), &
269 : DIMENSION(:), INTENT(in) :: sub_atom_kind_name
270 : TYPE(mp_para_env_type), POINTER :: para_env
271 : TYPE(section_vals_type), POINTER :: force_env_section, subsys_section
272 : LOGICAL, INTENT(in), OPTIONAL :: ignore_outside_box
273 :
274 : CHARACTER(len=default_string_length) :: my_element, strtmp1
275 : INTEGER :: iat, id_, nat
276 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
277 538 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
278 : TYPE(molecule_kind_list_type), POINTER :: mol_kinds
279 538 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
280 : TYPE(molecule_list_type), POINTER :: mols
281 538 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
282 : TYPE(particle_list_type), POINTER :: particles
283 538 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
284 : TYPE(topology_parameters_type) :: topology
285 :
286 538 : NULLIFY (mol_kinds, mols, particles, atomic_kinds, atomic_kind_set, particle_set, &
287 538 : molecule_kind_set, molecule_set, particles, atomic_kinds)
288 :
289 0 : CPASSERT(.NOT. ASSOCIATED(small_subsys))
290 538 : CPASSERT(ASSOCIATED(big_subsys))
291 538 : IF (big_subsys%para_env /= small_para_env) &
292 0 : CPABORT("big_subsys%para_env==small_para_env")
293 :
294 : !-----------------------------------------------------------------------------
295 : !-----------------------------------------------------------------------------
296 : ! 1. Initialize the topology structure type
297 : !-----------------------------------------------------------------------------
298 538 : CALL init_topology(topology)
299 :
300 : !-----------------------------------------------------------------------------
301 : !-----------------------------------------------------------------------------
302 : ! 2. Get the cell info
303 : !-----------------------------------------------------------------------------
304 538 : topology%cell => small_cell
305 538 : CALL cell_retain(small_cell)
306 538 : CALL write_cell(small_cell, subsys_section, tag="CELL")
307 538 : CALL write_cell(small_cell, subsys_section, tag="CELL_REF")
308 :
309 : !-----------------------------------------------------------------------------
310 : !-----------------------------------------------------------------------------
311 : ! 3. Initialize atom coords from the bigger system
312 : !-----------------------------------------------------------------------------
313 538 : nat = SIZE(sub_atom_index)
314 538 : topology%natoms = nat
315 538 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%r))
316 538 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_atmname))
317 538 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_molname))
318 538 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_resname))
319 538 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_mass))
320 538 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_charge))
321 : ALLOCATE (topology%atom_info%r(3, nat), topology%atom_info%id_atmname(nat), &
322 : topology%atom_info%id_molname(nat), topology%atom_info%id_resname(nat), &
323 : topology%atom_info%id_element(nat), topology%atom_info%atm_mass(nat), &
324 5918 : topology%atom_info%atm_charge(nat))
325 :
326 538 : CALL cp_subsys_get(big_subsys, particles=particles)
327 4250 : DO iat = 1, nat
328 14848 : topology%atom_info%r(:, iat) = particles%els(sub_atom_index(iat))%r
329 3712 : topology%atom_info%id_atmname(iat) = str2id(s2s(sub_atom_kind_name(iat)))
330 3712 : topology%atom_info%id_molname(iat) = topology%atom_info%id_atmname(iat)
331 3712 : topology%atom_info%id_resname(iat) = topology%atom_info%id_atmname(iat)
332 : !
333 : ! Defining element
334 : !
335 3712 : id_ = INDEX(id2str(topology%atom_info%id_atmname(iat)), "_") - 1
336 3712 : IF (id_ == -1) id_ = LEN_TRIM(id2str(topology%atom_info%id_atmname(iat)))
337 3712 : strtmp1 = id2str(topology%atom_info%id_atmname(iat))
338 3712 : strtmp1 = strtmp1(1:id_)
339 : CALL check_subsys_element(strtmp1, strtmp1, my_element, &
340 3712 : subsys_section, use_mm_map_first=.FALSE.)
341 3712 : topology%atom_info%id_element(iat) = str2id(s2s(my_element))
342 3712 : topology%atom_info%atm_mass(iat) = 0._dp
343 4250 : topology%atom_info%atm_charge(iat) = 0._dp
344 : END DO
345 538 : topology%conn_type = do_conn_off
346 :
347 : !-----------------------------------------------------------------------------
348 : !-----------------------------------------------------------------------------
349 : ! 4. Read in or generate the molecular connectivity
350 : !-----------------------------------------------------------------------------
351 : CALL connectivity_control(topology, para_env, subsys_section=subsys_section, &
352 538 : force_env_section=force_env_section)
353 :
354 : !-----------------------------------------------------------------------------
355 : !-----------------------------------------------------------------------------
356 : ! 5. Pack everything into the molecular types
357 : !-----------------------------------------------------------------------------
358 : CALL topology_connectivity_pack(molecule_kind_set, molecule_set, &
359 538 : topology, subsys_section=subsys_section)
360 :
361 : !-----------------------------------------------------------------------------
362 : !-----------------------------------------------------------------------------
363 : ! 6. Pack everything into the atomic types
364 : !-----------------------------------------------------------------------------
365 : CALL topology_coordinate_pack(particle_set, atomic_kind_set, &
366 : molecule_kind_set, molecule_set, topology, subsys_section=subsys_section, &
367 538 : force_env_section=force_env_section, ignore_outside_box=ignore_outside_box)
368 :
369 : !-----------------------------------------------------------------------------
370 : !-----------------------------------------------------------------------------
371 : ! 7. Cleanup the topology structure type
372 : !-----------------------------------------------------------------------------
373 538 : CALL deallocate_topology(topology)
374 :
375 : !-----------------------------------------------------------------------------
376 : !-----------------------------------------------------------------------------
377 : ! 8. Allocate new subsys
378 : !-----------------------------------------------------------------------------
379 4842 : ALLOCATE (small_subsys)
380 538 : CALL para_env%retain()
381 538 : small_subsys%para_env => para_env
382 538 : CALL particle_list_create(particles, els_ptr=particle_set)
383 538 : CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
384 538 : CALL molecule_list_create(mols, els_ptr=molecule_set)
385 538 : CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
386 : CALL cp_subsys_set(small_subsys, particles=particles, atomic_kinds=atomic_kinds, &
387 : molecules=mols, molecule_kinds=mol_kinds, cell=small_cell, &
388 538 : cell_ref=small_cell, use_ref_cell=.FALSE.)
389 538 : CALL particle_list_release(particles)
390 538 : CALL atomic_kind_list_release(atomic_kinds)
391 538 : CALL molecule_list_release(mols)
392 538 : CALL molecule_kind_list_release(mol_kinds)
393 :
394 123202 : ALLOCATE (small_subsys%virial)
395 538 : CALL atprop_create(small_subsys%atprop)
396 538 : CALL cp_result_create(small_subsys%results)
397 538 : END SUBROUTINE create_small_subsys
398 :
399 : END MODULE cp_subsys_methods
|
