Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 35516 : 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 8879 : INTEGER, DIMENSION(:), POINTER :: seed_vals
109 : LOGICAL :: atomic_energy, atomic_stress, &
110 : my_use_motion_section, &
111 : pv_availability, pv_diagonal, &
112 : pv_numerical
113 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
114 8879 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
115 : TYPE(molecule_kind_list_type), POINTER :: mol_kinds
116 8879 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
117 : TYPE(molecule_list_type), POINTER :: mols
118 8879 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
119 : TYPE(particle_list_type), POINTER :: particles
120 8879 : 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 88790 : ALLOCATE (subsys)
126 :
127 8879 : CALL para_env%retain()
128 8879 : subsys%para_env => para_env
129 :
130 8879 : my_use_motion_section = .FALSE.
131 8879 : IF (PRESENT(use_motion_section)) &
132 8877 : my_use_motion_section = use_motion_section
133 :
134 8879 : my_force_env_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
135 8879 : IF (PRESENT(force_env_section)) &
136 8877 : my_force_env_section => force_env_section
137 :
138 8879 : my_subsys_section => section_vals_get_subs_vals(my_force_env_section, "SUBSYS")
139 8879 : IF (PRESENT(subsys_section)) &
140 8723 : my_subsys_section => subsys_section
141 :
142 8879 : CALL section_vals_val_get(my_subsys_section, "SEED", i_vals=seed_vals)
143 8879 : IF (SIZE(seed_vals) == 1) THEN
144 79875 : 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 8879 : colvar_section => section_vals_get_subs_vals(my_subsys_section, "COLVAR")
152 :
153 8879 : 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 8879 : qmmm=qmmm, qmmm_env=qmmm_env, exclusions=exclusions, elkind=elkind)
162 :
163 8879 : CALL particle_list_create(particles, els_ptr=particle_set)
164 8879 : CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
165 8879 : CALL molecule_list_create(mols, els_ptr=molecule_set)
166 8879 : CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
167 :
168 : CALL cp_subsys_set(subsys, particles=particles, atomic_kinds=atomic_kinds, &
169 8879 : molecules=mols, molecule_kinds=mol_kinds)
170 :
171 8879 : CALL particle_list_release(particles)
172 8879 : CALL atomic_kind_list_release(atomic_kinds)
173 8879 : CALL molecule_list_release(mols)
174 8879 : CALL molecule_kind_list_release(mol_kinds)
175 :
176 : ! Should we compute the virial?
177 8879 : CALL section_vals_val_get(my_force_env_section, "STRESS_TENSOR", i_val=stress_tensor)
178 8109 : SELECT CASE (stress_tensor)
179 : CASE (do_stress_none)
180 8109 : pv_availability = .FALSE.
181 8109 : pv_numerical = .FALSE.
182 8109 : pv_diagonal = .FALSE.
183 : CASE (do_stress_analytical)
184 758 : pv_availability = .TRUE.
185 758 : pv_numerical = .FALSE.
186 758 : pv_diagonal = .FALSE.
187 : CASE (do_stress_numerical)
188 2 : pv_availability = .TRUE.
189 2 : pv_numerical = .TRUE.
190 2 : pv_diagonal = .FALSE.
191 : CASE (do_stress_diagonal_anal)
192 0 : pv_availability = .TRUE.
193 0 : pv_numerical = .FALSE.
194 0 : pv_diagonal = .TRUE.
195 : CASE (do_stress_diagonal_numer)
196 10 : pv_availability = .TRUE.
197 10 : pv_numerical = .TRUE.
198 8879 : pv_diagonal = .TRUE.
199 : END SELECT
200 :
201 2201992 : ALLOCATE (subsys%virial)
202 : CALL virial_set(virial=subsys%virial, &
203 : pv_availability=pv_availability, &
204 : pv_numer=pv_numerical, &
205 8879 : pv_diagonal=pv_diagonal)
206 :
207 : ! Should we compute atomic properties?
208 8879 : CALL atprop_create(subsys%atprop)
209 8879 : CALL section_vals_val_get(my_force_env_section, "PROPERTIES%ATOMIC%ENERGY", l_val=atomic_energy)
210 8879 : subsys%atprop%energy = atomic_energy
211 8879 : CALL section_vals_val_get(my_force_env_section, "PROPERTIES%ATOMIC%PRESSURE", l_val=atomic_stress)
212 8879 : IF (atomic_stress) THEN
213 22 : CPASSERT(pv_availability)
214 22 : CPASSERT(.NOT. pv_numerical)
215 : END IF
216 8879 : subsys%atprop%stress = atomic_stress
217 :
218 8879 : CALL cp_result_create(subsys%results)
219 8879 : END SUBROUTINE cp_subsys_create
220 :
221 : ! **************************************************************************************************
222 : !> \brief reads the colvar section of the colvar
223 : !> \param subsys ...
224 : !> \param colvar_section ...
225 : !> \par History
226 : !> 2006.01 Joost VandeVondele
227 : ! **************************************************************************************************
228 8879 : SUBROUTINE cp_subsys_read_colvar(subsys, colvar_section)
229 : TYPE(cp_subsys_type), POINTER :: subsys
230 : TYPE(section_vals_type), POINTER :: colvar_section
231 :
232 : INTEGER :: ig, ncol
233 :
234 8879 : CALL section_vals_get(colvar_section, n_repetition=ncol)
235 18492 : ALLOCATE (subsys%colvar_p(ncol))
236 9341 : DO ig = 1, ncol
237 462 : NULLIFY (subsys%colvar_p(ig)%colvar)
238 9341 : CALL colvar_read(subsys%colvar_p(ig)%colvar, ig, colvar_section, subsys%para_env)
239 : END DO
240 8879 : END SUBROUTINE cp_subsys_read_colvar
241 :
242 : ! **************************************************************************************************
243 : !> \brief updates the molecule information of the given subsys
244 : !> \param small_subsys the subsys to create
245 : !> \param big_subsys the superset of small_subsys
246 : !> \param small_cell ...
247 : !> \param small_para_env the parallel environment for the new (small)
248 : !> subsys
249 : !> \param sub_atom_index indexes of the atoms that should be in small_subsys
250 : !> \param sub_atom_kind_name ...
251 : !> \param para_env ...
252 : !> \param force_env_section ...
253 : !> \param subsys_section ...
254 : !> \param ignore_outside_box ...
255 : !> \par History
256 : !> 05.2004 created [fawzi]
257 : !> \author Fawzi Mohamed, Teodoro Laino
258 : !> \note
259 : !> not really ready to be used with different para_envs for the small
260 : !> and big part
261 : ! **************************************************************************************************
262 528 : SUBROUTINE create_small_subsys(small_subsys, big_subsys, small_cell, &
263 528 : small_para_env, sub_atom_index, sub_atom_kind_name, &
264 : para_env, force_env_section, subsys_section, ignore_outside_box)
265 :
266 : TYPE(cp_subsys_type), POINTER :: small_subsys, big_subsys
267 : TYPE(cell_type), POINTER :: small_cell
268 : TYPE(mp_para_env_type), POINTER :: small_para_env
269 : INTEGER, DIMENSION(:), INTENT(in) :: sub_atom_index
270 : CHARACTER(len=default_string_length), &
271 : DIMENSION(:), INTENT(in) :: sub_atom_kind_name
272 : TYPE(mp_para_env_type), POINTER :: para_env
273 : TYPE(section_vals_type), POINTER :: force_env_section, subsys_section
274 : LOGICAL, INTENT(in), OPTIONAL :: ignore_outside_box
275 :
276 : CHARACTER(len=default_string_length) :: my_element, strtmp1
277 : INTEGER :: iat, id_, nat
278 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
279 528 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
280 : TYPE(molecule_kind_list_type), POINTER :: mol_kinds
281 528 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
282 : TYPE(molecule_list_type), POINTER :: mols
283 528 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
284 : TYPE(particle_list_type), POINTER :: particles
285 528 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
286 : TYPE(topology_parameters_type) :: topology
287 :
288 528 : NULLIFY (mol_kinds, mols, particles, atomic_kinds, atomic_kind_set, particle_set, &
289 528 : molecule_kind_set, molecule_set, particles, atomic_kinds)
290 :
291 0 : CPASSERT(.NOT. ASSOCIATED(small_subsys))
292 528 : CPASSERT(ASSOCIATED(big_subsys))
293 528 : IF (big_subsys%para_env /= small_para_env) &
294 0 : CPABORT("big_subsys%para_env==small_para_env")
295 :
296 : !-----------------------------------------------------------------------------
297 : !-----------------------------------------------------------------------------
298 : ! 1. Initialize the topology structure type
299 : !-----------------------------------------------------------------------------
300 528 : CALL init_topology(topology)
301 :
302 : !-----------------------------------------------------------------------------
303 : !-----------------------------------------------------------------------------
304 : ! 2. Get the cell info
305 : !-----------------------------------------------------------------------------
306 528 : topology%cell => small_cell
307 528 : CALL cell_retain(small_cell)
308 :
309 : !-----------------------------------------------------------------------------
310 : !-----------------------------------------------------------------------------
311 : ! 3. Initialize atom coords from the bigger system
312 : !-----------------------------------------------------------------------------
313 528 : nat = SIZE(sub_atom_index)
314 528 : topology%natoms = nat
315 528 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%r))
316 528 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_atmname))
317 528 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_molname))
318 528 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_resname))
319 528 : CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_mass))
320 528 : 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 5808 : topology%atom_info%atm_charge(nat))
325 :
326 528 : CALL cp_subsys_get(big_subsys, particles=particles)
327 4192 : DO iat = 1, nat
328 14656 : topology%atom_info%r(:, iat) = particles%els(sub_atom_index(iat))%r
329 3664 : topology%atom_info%id_atmname(iat) = str2id(s2s(sub_atom_kind_name(iat)))
330 3664 : topology%atom_info%id_molname(iat) = topology%atom_info%id_atmname(iat)
331 3664 : topology%atom_info%id_resname(iat) = topology%atom_info%id_atmname(iat)
332 : !
333 : ! Defining element
334 : !
335 3664 : id_ = INDEX(id2str(topology%atom_info%id_atmname(iat)), "_") - 1
336 3664 : IF (id_ == -1) id_ = LEN_TRIM(id2str(topology%atom_info%id_atmname(iat)))
337 3664 : strtmp1 = id2str(topology%atom_info%id_atmname(iat))
338 3664 : strtmp1 = strtmp1(1:id_)
339 : CALL check_subsys_element(strtmp1, strtmp1, my_element, &
340 3664 : subsys_section, use_mm_map_first=.FALSE.)
341 3664 : topology%atom_info%id_element(iat) = str2id(s2s(my_element))
342 3664 : topology%atom_info%atm_mass(iat) = 0._dp
343 4192 : topology%atom_info%atm_charge(iat) = 0._dp
344 : END DO
345 528 : 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 528 : 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 528 : 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 528 : force_env_section=force_env_section, ignore_outside_box=ignore_outside_box)
368 :
369 : !-----------------------------------------------------------------------------
370 : !-----------------------------------------------------------------------------
371 : ! 7. Cleanup the topology structure type
372 : !-----------------------------------------------------------------------------
373 528 : CALL deallocate_topology(topology)
374 :
375 : !-----------------------------------------------------------------------------
376 : !-----------------------------------------------------------------------------
377 : ! 8. Allocate new subsys
378 : !-----------------------------------------------------------------------------
379 4752 : ALLOCATE (small_subsys)
380 528 : CALL para_env%retain()
381 528 : small_subsys%para_env => para_env
382 528 : CALL particle_list_create(particles, els_ptr=particle_set)
383 528 : CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
384 528 : CALL molecule_list_create(mols, els_ptr=molecule_set)
385 528 : 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 528 : molecules=mols, molecule_kinds=mol_kinds)
388 528 : CALL particle_list_release(particles)
389 528 : CALL atomic_kind_list_release(atomic_kinds)
390 528 : CALL molecule_list_release(mols)
391 528 : CALL molecule_kind_list_release(mol_kinds)
392 :
393 120912 : ALLOCATE (small_subsys%virial)
394 528 : CALL atprop_create(small_subsys%atprop)
395 528 : CALL cp_result_create(small_subsys%results)
396 528 : END SUBROUTINE create_small_subsys
397 :
398 : END MODULE cp_subsys_methods
|