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 32841 : 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 10947 : 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 10947 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
115 : TYPE(molecule_kind_list_type), POINTER :: mol_kinds
116 10947 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
117 : TYPE(molecule_list_type), POINTER :: mols
118 10947 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
119 : TYPE(particle_list_type), POINTER :: particles
120 10947 : 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 109470 : ALLOCATE (subsys)
126 :
127 10947 : CALL para_env%retain()
128 10947 : subsys%para_env => para_env
129 :
130 10947 : my_use_motion_section = .FALSE.
131 10947 : IF (PRESENT(use_motion_section)) THEN
132 10939 : my_use_motion_section = use_motion_section
133 : END IF
134 :
135 10947 : my_force_env_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
136 10947 : IF (PRESENT(force_env_section)) THEN
137 10939 : my_force_env_section => force_env_section
138 : END IF
139 :
140 10947 : my_subsys_section => section_vals_get_subs_vals(my_force_env_section, "SUBSYS")
141 10947 : IF (PRESENT(subsys_section)) THEN
142 10779 : my_subsys_section => subsys_section
143 : END IF
144 :
145 10947 : CALL section_vals_val_get(my_subsys_section, "SEED", i_vals=seed_vals)
146 10947 : IF (SIZE(seed_vals) == 1) THEN
147 98487 : subsys%seed(:, :) = REAL(seed_vals(1), KIND=dp)
148 4 : ELSE IF (SIZE(seed_vals) == 6) THEN
149 60 : subsys%seed(1:3, 1:2) = RESHAPE(REAL(seed_vals(:), KIND=dp), [3, 2])
150 : ELSE
151 0 : CPABORT("Supply exactly 1 or 6 arguments for SEED in &SUBSYS only!")
152 : END IF
153 :
154 10947 : colvar_section => section_vals_get_subs_vals(my_subsys_section, "COLVAR")
155 :
156 10947 : CALL cp_subsys_read_colvar(subsys, colvar_section)
157 :
158 : ! *** Read the particle coordinates and allocate the atomic kind, ***
159 : ! *** the molecule kind, and the molecule data structures ***
160 : CALL topology_control(atomic_kind_set, particle_set, molecule_kind_set, molecule_set, &
161 : subsys%colvar_p, subsys%gci, root_section, para_env, &
162 : force_env_section=my_force_env_section, &
163 : subsys_section=my_subsys_section, use_motion_section=my_use_motion_section, &
164 : qmmm=qmmm, qmmm_env=qmmm_env, exclusions=exclusions, elkind=elkind, &
165 10947 : subsys=subsys)
166 :
167 10947 : CALL particle_list_create(particles, els_ptr=particle_set)
168 10947 : CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
169 10947 : CALL molecule_list_create(mols, els_ptr=molecule_set)
170 10947 : CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
171 :
172 : CALL cp_subsys_set(subsys, particles=particles, atomic_kinds=atomic_kinds, &
173 10947 : molecules=mols, molecule_kinds=mol_kinds)
174 :
175 10947 : CALL particle_list_release(particles)
176 10947 : CALL atomic_kind_list_release(atomic_kinds)
177 10947 : CALL molecule_list_release(mols)
178 10947 : CALL molecule_kind_list_release(mol_kinds)
179 :
180 : ! Should we compute the virial?
181 10947 : CALL section_vals_val_get(my_force_env_section, "STRESS_TENSOR", i_val=stress_tensor)
182 9881 : SELECT CASE (stress_tensor)
183 : CASE (do_stress_none)
184 9881 : pv_availability = .FALSE.
185 9881 : pv_numerical = .FALSE.
186 9881 : pv_diagonal = .FALSE.
187 : CASE (do_stress_analytical)
188 1012 : pv_availability = .TRUE.
189 1012 : pv_numerical = .FALSE.
190 1012 : pv_diagonal = .FALSE.
191 : CASE (do_stress_numerical)
192 2 : pv_availability = .TRUE.
193 2 : pv_numerical = .TRUE.
194 2 : pv_diagonal = .FALSE.
195 : CASE (do_stress_diagonal_anal)
196 42 : pv_availability = .TRUE.
197 42 : pv_numerical = .FALSE.
198 42 : pv_diagonal = .TRUE.
199 : CASE (do_stress_diagonal_numer)
200 10 : pv_availability = .TRUE.
201 10 : pv_numerical = .TRUE.
202 10947 : pv_diagonal = .TRUE.
203 : END SELECT
204 :
205 2714856 : ALLOCATE (subsys%virial)
206 : CALL virial_set(virial=subsys%virial, &
207 : pv_availability=pv_availability, &
208 : pv_numer=pv_numerical, &
209 10947 : pv_diagonal=pv_diagonal)
210 :
211 : ! Should we compute atomic properties?
212 10947 : CALL atprop_create(subsys%atprop)
213 10947 : CALL section_vals_val_get(my_force_env_section, "PROPERTIES%ATOMIC%ENERGY", l_val=atomic_energy)
214 10947 : subsys%atprop%energy = atomic_energy
215 :
216 10947 : CALL cp_result_create(subsys%results)
217 10947 : END SUBROUTINE cp_subsys_create
218 :
219 : ! **************************************************************************************************
220 : !> \brief reads the colvar section of the colvar
221 : !> \param subsys ...
222 : !> \param colvar_section ...
223 : !> \par History
224 : !> 2006.01 Joost VandeVondele
225 : ! **************************************************************************************************
226 10947 : SUBROUTINE cp_subsys_read_colvar(subsys, colvar_section)
227 : TYPE(cp_subsys_type), POINTER :: subsys
228 : TYPE(section_vals_type), POINTER :: colvar_section
229 :
230 : INTEGER :: ig, ncol
231 :
232 10947 : CALL section_vals_get(colvar_section, n_repetition=ncol)
233 22636 : ALLOCATE (subsys%colvar_p(ncol))
234 11413 : DO ig = 1, ncol
235 466 : NULLIFY (subsys%colvar_p(ig)%colvar)
236 11413 : CALL colvar_read(subsys%colvar_p(ig)%colvar, ig, colvar_section, subsys%para_env, subsys%cell)
237 : END DO
238 10947 : END SUBROUTINE cp_subsys_read_colvar
239 :
240 : ! **************************************************************************************************
241 : !> \brief updates the molecule information of the given subsys
242 : !> \param small_subsys the subsys to create
243 : !> \param big_subsys the superset of small_subsys
244 : !> \param small_cell the cell of small_subsys
245 : !> \param small_para_env the parallel environment for the new (small)
246 : !> subsys
247 : !> \param sub_atom_index indexes of the atoms that should be in small_subsys
248 : !> \param sub_atom_kind_name ...
249 : !> \param para_env ...
250 : !> \param force_env_section ...
251 : !> \param subsys_section ...
252 : !> \param ignore_outside_box ...
253 : !> \par History
254 : !> 05.2004 created [fawzi]
255 : !> \author Fawzi Mohamed, Teodoro Laino
256 : !> \note
257 : !> not really ready to be used with different para_envs for the small
258 : !> and big part
259 : !> qmmm_env_create() should be the only instance where this subroutine is called
260 : !> with small_cell distinct from big_subsys%cell; otherwise it could be possible
261 : !> to drop the small_cell as input parameter entirely.
262 : ! **************************************************************************************************
263 550 : SUBROUTINE create_small_subsys(small_subsys, big_subsys, small_cell, &
264 550 : small_para_env, sub_atom_index, sub_atom_kind_name, &
265 : para_env, force_env_section, subsys_section, ignore_outside_box)
266 :
267 : TYPE(cp_subsys_type), POINTER :: small_subsys, big_subsys
268 : TYPE(cell_type), POINTER :: small_cell
269 : TYPE(mp_para_env_type), POINTER :: small_para_env
270 : INTEGER, DIMENSION(:), INTENT(in) :: sub_atom_index
271 : CHARACTER(len=default_string_length), &
272 : DIMENSION(:), INTENT(in) :: sub_atom_kind_name
273 : TYPE(mp_para_env_type), POINTER :: para_env
274 : TYPE(section_vals_type), POINTER :: force_env_section, subsys_section
275 : LOGICAL, INTENT(in), OPTIONAL :: ignore_outside_box
276 :
277 : CHARACTER(len=default_string_length) :: my_element, strtmp1
278 : INTEGER :: iat, id_, nat
279 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
280 550 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
281 : TYPE(molecule_kind_list_type), POINTER :: mol_kinds
282 550 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
283 : TYPE(molecule_list_type), POINTER :: mols
284 550 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
285 : TYPE(particle_list_type), POINTER :: particles
286 550 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
287 : TYPE(topology_parameters_type) :: topology
288 :
289 550 : NULLIFY (mol_kinds, mols, particles, atomic_kinds, atomic_kind_set, particle_set, &
290 550 : molecule_kind_set, molecule_set, particles, atomic_kinds)
291 :
292 0 : CPASSERT(.NOT. ASSOCIATED(small_subsys))
293 550 : CPASSERT(ASSOCIATED(big_subsys))
294 550 : IF (big_subsys%para_env /= small_para_env) THEN
295 0 : CPABORT("big_subsys%para_env==small_para_env")
296 : END IF
297 :
298 : !-----------------------------------------------------------------------------
299 : !-----------------------------------------------------------------------------
300 : ! 1. Initialize the topology structure type
301 : !-----------------------------------------------------------------------------
302 550 : CALL init_topology(topology)
303 :
304 : !-----------------------------------------------------------------------------
305 : !-----------------------------------------------------------------------------
306 : ! 2. Get the cell info
307 : !-----------------------------------------------------------------------------
308 550 : topology%cell => small_cell
309 550 : CALL cell_retain(small_cell)
310 550 : CALL write_cell(small_cell, subsys_section, tag="CELL")
311 550 : CALL write_cell(small_cell, subsys_section, tag="CELL_REF")
312 :
313 : !-----------------------------------------------------------------------------
314 : !-----------------------------------------------------------------------------
315 : ! 3. Initialize atom coords from the bigger system
316 : !-----------------------------------------------------------------------------
317 550 : nat = SIZE(sub_atom_index)
318 550 : topology%natoms = nat
319 550 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%r))
320 550 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_atmname))
321 550 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_molname))
322 550 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_resname))
323 550 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_mass))
324 550 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_charge))
325 : ALLOCATE (topology%atom_info%r(3, nat), topology%atom_info%id_atmname(nat), &
326 : topology%atom_info%id_molname(nat), topology%atom_info%id_resname(nat), &
327 : topology%atom_info%id_element(nat), topology%atom_info%atm_mass(nat), &
328 6050 : topology%atom_info%atm_charge(nat))
329 :
330 550 : CALL cp_subsys_get(big_subsys, particles=particles)
331 4290 : DO iat = 1, nat
332 14960 : topology%atom_info%r(:, iat) = particles%els(sub_atom_index(iat))%r
333 3740 : topology%atom_info%id_atmname(iat) = str2id(s2s(sub_atom_kind_name(iat)))
334 3740 : topology%atom_info%id_molname(iat) = topology%atom_info%id_atmname(iat)
335 3740 : topology%atom_info%id_resname(iat) = topology%atom_info%id_atmname(iat)
336 : !
337 : ! Defining element
338 : !
339 3740 : id_ = INDEX(id2str(topology%atom_info%id_atmname(iat)), "_") - 1
340 3740 : IF (id_ == -1) id_ = LEN_TRIM(id2str(topology%atom_info%id_atmname(iat)))
341 3740 : strtmp1 = id2str(topology%atom_info%id_atmname(iat))
342 3740 : strtmp1 = strtmp1(1:id_)
343 : CALL check_subsys_element(strtmp1, strtmp1, my_element, &
344 3740 : subsys_section, use_mm_map_first=.FALSE.)
345 3740 : topology%atom_info%id_element(iat) = str2id(s2s(my_element))
346 3740 : topology%atom_info%atm_mass(iat) = 0._dp
347 4290 : topology%atom_info%atm_charge(iat) = 0._dp
348 : END DO
349 550 : topology%conn_type = do_conn_off
350 :
351 : !-----------------------------------------------------------------------------
352 : !-----------------------------------------------------------------------------
353 : ! 4. Read in or generate the molecular connectivity
354 : !-----------------------------------------------------------------------------
355 : CALL connectivity_control(topology, para_env, subsys_section=subsys_section, &
356 550 : force_env_section=force_env_section)
357 :
358 : !-----------------------------------------------------------------------------
359 : !-----------------------------------------------------------------------------
360 : ! 5. Pack everything into the molecular types
361 : !-----------------------------------------------------------------------------
362 : CALL topology_connectivity_pack(molecule_kind_set, molecule_set, &
363 550 : topology, subsys_section=subsys_section)
364 :
365 : !-----------------------------------------------------------------------------
366 : !-----------------------------------------------------------------------------
367 : ! 6. Pack everything into the atomic types
368 : !-----------------------------------------------------------------------------
369 : CALL topology_coordinate_pack(particle_set, atomic_kind_set, &
370 : molecule_kind_set, molecule_set, topology, subsys_section=subsys_section, &
371 550 : force_env_section=force_env_section, ignore_outside_box=ignore_outside_box)
372 :
373 : !-----------------------------------------------------------------------------
374 : !-----------------------------------------------------------------------------
375 : ! 7. Cleanup the topology structure type
376 : !-----------------------------------------------------------------------------
377 550 : CALL deallocate_topology(topology)
378 :
379 : !-----------------------------------------------------------------------------
380 : !-----------------------------------------------------------------------------
381 : ! 8. Allocate new subsys
382 : !-----------------------------------------------------------------------------
383 4950 : ALLOCATE (small_subsys)
384 550 : CALL para_env%retain()
385 550 : small_subsys%para_env => para_env
386 550 : CALL particle_list_create(particles, els_ptr=particle_set)
387 550 : CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
388 550 : CALL molecule_list_create(mols, els_ptr=molecule_set)
389 550 : CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
390 : CALL cp_subsys_set(small_subsys, particles=particles, atomic_kinds=atomic_kinds, &
391 : molecules=mols, molecule_kinds=mol_kinds, cell=small_cell, &
392 550 : cell_ref=small_cell, use_ref_cell=.FALSE.)
393 550 : CALL particle_list_release(particles)
394 550 : CALL atomic_kind_list_release(atomic_kinds)
395 550 : CALL molecule_list_release(mols)
396 550 : CALL molecule_kind_list_release(mol_kinds)
397 :
398 125950 : ALLOCATE (small_subsys%virial)
399 550 : CALL atprop_create(small_subsys%atprop)
400 550 : CALL cp_result_create(small_subsys%results)
401 550 : END SUBROUTINE create_small_subsys
402 :
403 : END MODULE cp_subsys_methods
|
