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