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 : MODULE input_restart_force_eval
9 :
10 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
11 : USE atomic_kind_types, ONLY: get_atomic_kind
12 : USE cell_types, ONLY: cell_type,&
13 : real_to_scaled
14 : USE cp_control_types, ONLY: dft_control_type
15 : USE cp_files, ONLY: close_file,&
16 : open_file
17 : USE cp_linked_list_input, ONLY: cp_sll_val_create,&
18 : cp_sll_val_get_length,&
19 : cp_sll_val_type
20 : USE cp_subsys_types, ONLY: cp_subsys_get,&
21 : cp_subsys_type
22 : USE cp_units, ONLY: cp_unit_from_cp2k
23 : USE distribution_1d_types, ONLY: distribution_1d_type
24 : USE eip_environment_types, ONLY: eip_env_get
25 : USE force_env_types, ONLY: force_env_get,&
26 : force_env_type,&
27 : multiple_fe_list,&
28 : use_eip_force,&
29 : use_nnp_force,&
30 : use_qmmm,&
31 : use_qs_force
32 : USE input_constants, ONLY: do_coord_cp2k,&
33 : ehrenfest,&
34 : mol_dyn_run,&
35 : mon_car_run,&
36 : pint_run,&
37 : use_rt_restart
38 : USE input_cp2k_restarts_util, ONLY: section_velocity_val_set
39 : USE input_restart_rng, ONLY: section_rng_val_set
40 : USE input_section_types, ONLY: &
41 : section_get_keyword_index, section_type, section_vals_add_values, section_vals_get, &
42 : section_vals_get_subs_vals, section_vals_get_subs_vals3, section_vals_remove_values, &
43 : section_vals_type, section_vals_val_get, section_vals_val_set, section_vals_val_unset
44 : USE input_val_types, ONLY: val_create,&
45 : val_release,&
46 : val_type
47 : USE kinds, ONLY: default_string_length,&
48 : dp
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE molecule_kind_types, ONLY: get_molecule_kind
51 : USE molecule_list_types, ONLY: molecule_list_type
52 : USE molecule_types, ONLY: get_molecule,&
53 : molecule_type
54 : USE multipole_types, ONLY: multipole_type
55 : USE nnp_environment_types, ONLY: nnp_env_get
56 : USE parallel_rng_types, ONLY: rng_record_length
57 : USE particle_list_types, ONLY: particle_list_type
58 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
59 : USE qs_environment_types, ONLY: get_qs_env
60 : USE string_utilities, ONLY: string_to_ascii
61 : USE virial_types, ONLY: virial_type
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 :
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'input_restart_force_eval'
69 :
70 : PUBLIC :: update_force_eval, update_subsys
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief Updates the force_eval section of the input file
76 : !> \param force_env ...
77 : !> \param root_section ...
78 : !> \param write_binary_restart_file ...
79 : !> \param respa ...
80 : !> \par History
81 : !> 01.2006 created [teo]
82 : !> \author Teodoro Laino
83 : ! **************************************************************************************************
84 14240 : SUBROUTINE update_force_eval(force_env, root_section, &
85 : write_binary_restart_file, respa)
86 :
87 : TYPE(force_env_type), POINTER :: force_env
88 : TYPE(section_vals_type), POINTER :: root_section
89 : LOGICAL, INTENT(IN) :: write_binary_restart_file
90 : LOGICAL, OPTIONAL :: respa
91 :
92 : INTEGER :: i, i_subsys, iforce_eval, myid, &
93 : nforce_eval
94 14240 : INTEGER, DIMENSION(:), POINTER :: i_force_eval
95 : LOGICAL :: is_present, multiple_subsys, &
96 : skip_vel_section
97 14240 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
98 : TYPE(cell_type), POINTER :: cell
99 : TYPE(cp_subsys_type), POINTER :: subsys
100 : TYPE(dft_control_type), POINTER :: dft_control
101 : TYPE(section_vals_type), POINTER :: cell_section, dft_section, &
102 : force_env_sections, qmmm_section, &
103 : rng_section, subsys_section, &
104 : tmp_section
105 : TYPE(virial_type), POINTER :: virial
106 :
107 14240 : NULLIFY (rng_section, subsys_section, cell_section, virial, subsys, cell, dft_control, work, tmp_section)
108 : ! If it's not a dynamical run don't update the velocity section
109 14240 : CALL section_vals_val_get(root_section, "GLOBAL%RUN_TYPE", i_val=myid)
110 : skip_vel_section = ((myid /= mol_dyn_run) .AND. &
111 : (myid /= mon_car_run) .AND. &
112 : (myid /= pint_run) .AND. &
113 14240 : (myid /= ehrenfest))
114 :
115 : ! Go on updatig the force_env_section
116 14240 : force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
117 14240 : CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, nforce_eval)
118 : ! The update of the input MUST be realized only on the main force_eval
119 : ! All the others will be left not updated because there is no real need to update them...
120 14240 : iforce_eval = 1
121 : subsys_section => section_vals_get_subs_vals3(force_env_sections, "SUBSYS", &
122 14240 : i_rep_section=i_force_eval(iforce_eval))
123 : CALL update_subsys(subsys_section, force_env, skip_vel_section, &
124 14240 : write_binary_restart_file)
125 :
126 : ! For RESPA we need to update all subsystems
127 14240 : IF (PRESENT(respa)) THEN
128 14218 : IF (respa) THEN
129 18 : DO i_subsys = 1, 2
130 12 : iforce_eval = i_subsys
131 : subsys_section => section_vals_get_subs_vals3(force_env_sections, "SUBSYS", &
132 12 : i_rep_section=i_force_eval(iforce_eval))
133 : CALL update_subsys(subsys_section, force_env, skip_vel_section, &
134 18 : write_binary_restart_file)
135 : END DO
136 : END IF
137 : END IF
138 :
139 14240 : rng_section => section_vals_get_subs_vals(subsys_section, "RNG_INIT")
140 14240 : CALL update_rng_particle(rng_section, force_env)
141 :
142 : qmmm_section => section_vals_get_subs_vals3(force_env_sections, "QMMM", &
143 14240 : i_rep_section=i_force_eval(iforce_eval))
144 14240 : CALL update_qmmm(qmmm_section, force_env)
145 :
146 : ! Non-mixed CDFT calculation: update constraint Lagrangian and counter
147 14240 : IF (nforce_eval == 1) THEN
148 14176 : IF (ASSOCIATED(force_env%qs_env)) THEN
149 2386 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
150 2386 : IF (dft_control%qs_control%cdft) THEN
151 : dft_section => section_vals_get_subs_vals3(force_env_sections, "DFT", &
152 16 : i_rep_section=i_force_eval(iforce_eval))
153 : CALL section_vals_val_set(dft_section, "QS%CDFT%COUNTER", &
154 16 : i_val=dft_control%qs_control%cdft_control%ienergy)
155 48 : ALLOCATE (work(SIZE(dft_control%qs_control%cdft_control%strength)))
156 80 : work = dft_control%qs_control%cdft_control%strength
157 : CALL section_vals_val_set(dft_section, "QS%CDFT%STRENGTH", &
158 16 : r_vals_ptr=work)
159 : END IF
160 : END IF
161 : END IF
162 : ! And now update only the cells of all other force_eval
163 : ! This is to make consistent for cell variable runs..
164 14240 : IF (nforce_eval > 1) THEN
165 64 : CALL force_env_get(force_env, subsys=subsys, cell=cell)
166 64 : CALL cp_subsys_get(subsys, virial=virial)
167 64 : CALL section_vals_val_get(root_section, "MULTIPLE_FORCE_EVALS%MULTIPLE_SUBSYS", l_val=multiple_subsys)
168 64 : IF (virial%pv_availability .AND. multiple_subsys) THEN
169 26 : DO iforce_eval = 2, nforce_eval
170 : subsys_section => section_vals_get_subs_vals3(force_env_sections, "SUBSYS", &
171 20 : i_rep_section=i_force_eval(iforce_eval))
172 20 : cell_section => section_vals_get_subs_vals(subsys_section, "CELL")
173 26 : CALL update_cell_section(cell, cell_section)
174 : END DO
175 : END IF
176 : ! With mixed CDFT, update value of constraint Lagrangian
177 64 : IF (ASSOCIATED(force_env%mixed_env)) THEN
178 58 : IF (force_env%mixed_env%do_mixed_cdft) THEN
179 36 : DO iforce_eval = 2, nforce_eval
180 : dft_section => section_vals_get_subs_vals3(force_env_sections, "DFT", &
181 24 : i_rep_section=i_force_eval(iforce_eval))
182 72 : ALLOCATE (work(SIZE(force_env%mixed_env%strength, 2)))
183 96 : work = force_env%mixed_env%strength(iforce_eval - 1, :)
184 : CALL section_vals_val_set(dft_section, "QS%CDFT%STRENGTH", &
185 24 : r_vals_ptr=work)
186 : CALL section_vals_val_set(dft_section, "QS%CDFT%COUNTER", &
187 36 : i_val=force_env%mixed_env%cdft_control%sim_step)
188 : END DO
189 : END IF
190 : END IF
191 : END IF
192 :
193 14240 : IF (ASSOCIATED(force_env%qs_env)) THEN
194 2390 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
195 : dft_section => section_vals_get_subs_vals(force_env_sections, "DFT", &
196 2390 : i_rep_section=i_force_eval(iforce_eval))
197 2390 : tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION")
198 2390 : CALL section_vals_get(tmp_section, explicit=is_present)
199 2390 : IF (is_present) THEN
200 : CALL section_vals_val_set(root_section, "FORCE_EVAL%DFT%REAL_TIME_PROPAGATION%INITIAL_WFN", &
201 158 : i_val=use_rt_restart)
202 158 : IF (ASSOCIATED(dft_control%efield_fields)) THEN
203 10 : tmp_section => section_vals_get_subs_vals(dft_section, "EFIELD")
204 10 : ALLOCATE (work(3))
205 20 : DO i = 1, SIZE(dft_control%efield_fields)
206 40 : work = dft_control%efield_fields(1)%efield%vec_pot_initial
207 : CALL section_vals_val_set(tmp_section, "VEC_POT_INITIAL", i_rep_section=i, &
208 20 : r_vals_ptr=work)
209 : END DO
210 : END IF
211 : END IF
212 : END IF
213 14240 : DEALLOCATE (i_force_eval)
214 :
215 14240 : END SUBROUTINE update_force_eval
216 :
217 : ! **************************************************************************************************
218 : !> \brief Updates the qmmm section if needed
219 : !> \param qmmm_section ...
220 : !> \param force_env ...
221 : !> \par History
222 : !> 08.2007 created [teo]
223 : !> \author Teodoro Laino
224 : ! **************************************************************************************************
225 14240 : SUBROUTINE update_qmmm(qmmm_section, force_env)
226 :
227 : TYPE(section_vals_type), POINTER :: qmmm_section
228 : TYPE(force_env_type), POINTER :: force_env
229 :
230 : LOGICAL :: explicit
231 14240 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
232 :
233 14494 : SELECT CASE (force_env%in_use)
234 : CASE (use_qmmm)
235 254 : CALL section_vals_get(qmmm_section, explicit=explicit)
236 254 : CPASSERT(explicit)
237 :
238 254 : ALLOCATE (work(3))
239 2032 : work = force_env%qmmm_env%qm%transl_v
240 508 : CALL section_vals_val_set(qmmm_section, "INITIAL_TRANSLATION_VECTOR", r_vals_ptr=work)
241 : END SELECT
242 :
243 14240 : END SUBROUTINE update_qmmm
244 :
245 : ! **************************************************************************************************
246 : !> \brief Updates the rng section of the input file
247 : !> Write current status of the parallel random number generator (RNG)
248 : !> \param rng_section ...
249 : !> \param force_env ...
250 : !> \par History
251 : !> 01.2006 created [teo]
252 : !> \author Teodoro Laino
253 : ! **************************************************************************************************
254 14240 : SUBROUTINE update_rng_particle(rng_section, force_env)
255 :
256 : TYPE(section_vals_type), POINTER :: rng_section
257 : TYPE(force_env_type), POINTER :: force_env
258 :
259 : CHARACTER(LEN=rng_record_length) :: rng_record
260 : INTEGER :: iparticle, iparticle_kind, &
261 : iparticle_local, nparticle, &
262 : nparticle_kind, nparticle_local
263 14240 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ascii
264 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
265 : TYPE(cp_subsys_type), POINTER :: subsys
266 : TYPE(distribution_1d_type), POINTER :: local_particles
267 : TYPE(mp_para_env_type), POINTER :: para_env
268 : TYPE(particle_list_type), POINTER :: particles
269 :
270 14240 : CALL force_env_get(force_env, subsys=subsys, para_env=para_env)
271 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
272 14240 : particles=particles)
273 :
274 14240 : IF (ASSOCIATED(local_particles%local_particle_set)) THEN
275 56 : nparticle_kind = atomic_kinds%n_els
276 56 : nparticle = particles%n_els
277 :
278 168 : ALLOCATE (ascii(rng_record_length, nparticle))
279 3915604 : ascii = 0
280 :
281 9078 : DO iparticle = 1, nparticle
282 159062 : DO iparticle_kind = 1, nparticle_kind
283 149984 : nparticle_local = local_particles%n_el(iparticle_kind)
284 13872895 : DO iparticle_local = 1, nparticle_local
285 13863873 : IF (iparticle == local_particles%list(iparticle_kind)%array(iparticle_local)) THEN
286 : CALL local_particles%local_particle_set(iparticle_kind)% &
287 4511 : rng(iparticle_local)%stream%dump(rng_record=rng_record)
288 4511 : CALL string_to_ascii(rng_record, ascii(:, iparticle))
289 : END IF
290 : END DO
291 : END DO
292 : END DO
293 :
294 56 : CALL para_env%sum(ascii)
295 :
296 56 : CALL section_rng_val_set(rng_section=rng_section, nsize=nparticle, ascii=ascii)
297 :
298 56 : DEALLOCATE (ascii)
299 : END IF
300 :
301 14240 : END SUBROUTINE update_rng_particle
302 :
303 : ! **************************************************************************************************
304 : !> \brief Updates the subsys section of the input file
305 : !> \param subsys_section ...
306 : !> \param force_env ...
307 : !> \param skip_vel_section ...
308 : !> \param write_binary_restart_file ...
309 : !> \par History
310 : !> 01.2006 created [teo]
311 : !> \author Teodoro Laino
312 : ! **************************************************************************************************
313 14278 : SUBROUTINE update_subsys(subsys_section, force_env, skip_vel_section, &
314 : write_binary_restart_file)
315 :
316 : TYPE(section_vals_type), POINTER :: subsys_section
317 : TYPE(force_env_type), POINTER :: force_env
318 : LOGICAL, INTENT(IN) :: skip_vel_section, &
319 : write_binary_restart_file
320 :
321 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_subsys'
322 :
323 : CHARACTER(LEN=default_string_length) :: coord_file_name, unit_str
324 : INTEGER :: coord_file_format, handle, output_unit
325 14278 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
326 : LOGICAL :: scale, use_ref_cell
327 : REAL(KIND=dp) :: conv_factor
328 : TYPE(cell_type), POINTER :: cell
329 : TYPE(cp_subsys_type), POINTER :: subsys
330 : TYPE(molecule_list_type), POINTER :: molecules
331 : TYPE(mp_para_env_type), POINTER :: para_env
332 : TYPE(multipole_type), POINTER :: multipoles
333 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
334 : shell_particles
335 : TYPE(section_vals_type), POINTER :: work_section
336 :
337 14278 : NULLIFY (work_section, core_particles, particles, shell_particles, &
338 14278 : subsys, cell, para_env, multipoles)
339 14278 : CALL timeset(routineN, handle)
340 14278 : CALL force_env_get(force_env, subsys=subsys, cell=cell, para_env=para_env)
341 :
342 : CALL cp_subsys_get(subsys, particles=particles, molecules=molecules, &
343 : shell_particles=shell_particles, core_particles=core_particles, &
344 14278 : multipoles=multipoles)
345 :
346 : ! Remove the multiple_unit_cell from the input structure.. at this point we have
347 : ! already all the information we need..
348 14278 : ALLOCATE (multiple_unit_cell(3))
349 57112 : multiple_unit_cell = 1
350 : CALL section_vals_val_set(subsys_section, "TOPOLOGY%MULTIPLE_UNIT_CELL", &
351 14278 : i_vals_ptr=multiple_unit_cell)
352 :
353 : ! Coordinates and Velocities
354 14278 : work_section => section_vals_get_subs_vals(subsys_section, "COORD")
355 14278 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
356 14278 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
357 14278 : conv_factor = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
358 14278 : IF (.NOT. write_binary_restart_file) THEN
359 : CALL section_vals_val_get(subsys_section, "TOPOLOGY%COORD_FILE_FORMAT", &
360 14142 : i_val=coord_file_format)
361 14142 : IF (coord_file_format == do_coord_cp2k) THEN
362 : CALL section_vals_val_get(subsys_section, "TOPOLOGY%COORD_FILE_NAME", &
363 0 : c_val=coord_file_name)
364 0 : output_unit = 0
365 0 : IF (para_env%is_source()) THEN
366 : CALL open_file(file_name=TRIM(ADJUSTL(coord_file_name)), &
367 : file_action="READWRITE", &
368 : file_form="FORMATTED", &
369 : file_position="REWIND", &
370 : file_status="UNKNOWN", &
371 0 : unit_number=output_unit)
372 : CALL dump_coordinates_cp2k(particles, molecules, cell, conv_factor, &
373 : output_unit=output_unit, &
374 : core_or_shell=.FALSE., &
375 0 : scaled_coordinates=scale)
376 0 : CALL close_file(unit_number=output_unit)
377 : END IF
378 : ELSE
379 : ! Remove COORD_FILE_NAME and COORD_FILE_FORMAT keywords from restart
380 14142 : CALL section_vals_val_unset(subsys_section, "TOPOLOGY%COORD_FILE_NAME")
381 14142 : CALL section_vals_val_unset(subsys_section, "TOPOLOGY%COORD_FILE_FORMAT")
382 : CALL section_coord_val_set(work_section, particles, molecules, conv_factor, scale, &
383 14142 : cell)
384 : END IF
385 : END IF
386 14278 : CALL section_vals_val_set(subsys_section, "TOPOLOGY%NUMBER_OF_ATOMS", i_val=particles%n_els)
387 14278 : work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
388 14278 : IF (.NOT. skip_vel_section) THEN
389 5574 : IF (.NOT. write_binary_restart_file) THEN
390 5438 : CALL section_velocity_val_set(work_section, particles, conv_factor=1.0_dp)
391 : END IF
392 : ELSE
393 8704 : CALL section_vals_remove_values(work_section)
394 : END IF
395 :
396 : ! Write restart input for core-shell model
397 14278 : IF (.NOT. write_binary_restart_file) THEN
398 14142 : IF (ASSOCIATED(shell_particles)) THEN
399 3302 : work_section => section_vals_get_subs_vals(subsys_section, "SHELL_COORD")
400 3302 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
401 3302 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
402 3302 : conv_factor = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
403 : CALL section_coord_val_set(work_section, shell_particles, molecules, &
404 3302 : conv_factor, scale, cell, shell=.TRUE.)
405 3302 : IF (.NOT. skip_vel_section) THEN
406 362 : work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
407 362 : CALL section_velocity_val_set(work_section, shell_particles, conv_factor=1.0_dp)
408 : END IF
409 : END IF
410 14142 : IF (ASSOCIATED(core_particles)) THEN
411 3302 : work_section => section_vals_get_subs_vals(subsys_section, "CORE_COORD")
412 3302 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
413 3302 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
414 3302 : conv_factor = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
415 : CALL section_coord_val_set(work_section, core_particles, molecules, &
416 3302 : conv_factor, scale, cell, shell=.TRUE.)
417 3302 : IF (.NOT. skip_vel_section) THEN
418 362 : work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
419 362 : CALL section_velocity_val_set(work_section, core_particles, conv_factor=1.0_dp)
420 : END IF
421 : END IF
422 : END IF
423 :
424 : ! Updating cell info
425 14278 : CALL force_env_get(force_env, cell=cell)
426 14278 : work_section => section_vals_get_subs_vals(subsys_section, "CELL")
427 14278 : CALL update_cell_section(cell, cell_section=work_section)
428 :
429 : ! Updating cell_ref info
430 14278 : use_ref_cell = .FALSE.
431 16676 : SELECT CASE (force_env%in_use)
432 : CASE (use_qs_force)
433 2398 : CALL get_qs_env(force_env%qs_env, cell_ref=cell, use_ref_cell=use_ref_cell)
434 : CASE (use_eip_force)
435 16 : CALL eip_env_get(force_env%eip_env, cell_ref=cell, use_ref_cell=use_ref_cell)
436 : CASE (use_nnp_force)
437 14278 : CALL nnp_env_get(force_env%nnp_env, cell_ref=cell, use_ref_cell=use_ref_cell)
438 : END SELECT
439 14278 : IF (use_ref_cell) THEN
440 34 : work_section => section_vals_get_subs_vals(subsys_section, "CELL%CELL_REF")
441 34 : CALL update_cell_section(cell, cell_section=work_section)
442 : END IF
443 :
444 : ! Updating multipoles
445 14278 : IF (ASSOCIATED(multipoles)) THEN
446 24 : work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES")
447 16 : DO
448 40 : IF (SIZE(work_section%values, 2) == 1) EXIT
449 16 : CALL section_vals_add_values(work_section)
450 : END DO
451 24 : IF (ASSOCIATED(multipoles%dipoles)) THEN
452 24 : work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES%DIPOLES")
453 24 : CALL update_dipoles_section(multipoles%dipoles, work_section)
454 : END IF
455 24 : IF (ASSOCIATED(multipoles%quadrupoles)) THEN
456 0 : work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES%QUADRUPOLES")
457 0 : CALL update_quadrupoles_section(multipoles%quadrupoles, work_section)
458 : END IF
459 : END IF
460 :
461 14278 : CALL timestop(handle)
462 :
463 28556 : END SUBROUTINE update_subsys
464 :
465 : ! **************************************************************************************************
466 : !> \brief Routine to update a cell section
467 : !> \param cell ...
468 : !> \param cell_section ...
469 : !> \author Ole Schuett
470 : ! **************************************************************************************************
471 14332 : SUBROUTINE update_cell_section(cell, cell_section)
472 :
473 : TYPE(cell_type), POINTER :: cell
474 : TYPE(section_vals_type), POINTER :: cell_section
475 :
476 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_cell_section'
477 :
478 : INTEGER :: handle
479 14332 : INTEGER, DIMENSION(:), POINTER :: iwork
480 14332 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
481 :
482 14332 : NULLIFY (work, iwork)
483 14332 : CALL timeset(routineN, handle)
484 :
485 : ! CELL VECTORS - A
486 14332 : ALLOCATE (work(3))
487 114656 : work(1:3) = cell%hmat(1:3, 1)
488 14332 : CALL section_vals_val_set(cell_section, "A", r_vals_ptr=work)
489 :
490 : ! CELL VECTORS - B
491 14332 : ALLOCATE (work(3))
492 114656 : work(1:3) = cell%hmat(1:3, 2)
493 14332 : CALL section_vals_val_set(cell_section, "B", r_vals_ptr=work)
494 :
495 : ! CELL VECTORS - C
496 14332 : ALLOCATE (work(3))
497 114656 : work(1:3) = cell%hmat(1:3, 3)
498 14332 : CALL section_vals_val_set(cell_section, "C", r_vals_ptr=work)
499 :
500 : ! MULTIPLE_UNIT_CELL
501 14332 : ALLOCATE (iwork(3))
502 57328 : iwork = 1
503 14332 : CALL section_vals_val_set(cell_section, "MULTIPLE_UNIT_CELL", i_vals_ptr=iwork)
504 :
505 : ! Unset unused or misleading information
506 14332 : CALL section_vals_val_unset(cell_section, "ABC")
507 14332 : CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA")
508 :
509 14332 : CALL timestop(handle)
510 :
511 14332 : END SUBROUTINE update_cell_section
512 :
513 : ! **************************************************************************************************
514 : !> \brief routine to dump coordinates.. fast implementation
515 : !> \param coord_section ...
516 : !> \param particles ...
517 : !> \param molecules ...
518 : !> \param conv_factor ...
519 : !> \param scale ...
520 : !> \param cell ...
521 : !> \param shell ...
522 : !> \par History
523 : !> 02.2006 created [teo]
524 : !> \author Teodoro Laino
525 : ! **************************************************************************************************
526 20746 : SUBROUTINE section_coord_val_set(coord_section, particles, molecules, conv_factor, &
527 : scale, cell, shell)
528 :
529 : TYPE(section_vals_type), POINTER :: coord_section
530 : TYPE(particle_list_type), POINTER :: particles
531 : TYPE(molecule_list_type), POINTER :: molecules
532 : REAL(KIND=dp) :: conv_factor
533 : LOGICAL, INTENT(IN) :: scale
534 : TYPE(cell_type), POINTER :: cell
535 : LOGICAL, INTENT(IN), OPTIONAL :: shell
536 :
537 : CHARACTER(LEN=*), PARAMETER :: routineN = 'section_coord_val_set'
538 :
539 : CHARACTER(LEN=2*default_string_length) :: line
540 : CHARACTER(LEN=default_string_length) :: my_tag, name
541 : INTEGER :: handle, ik, imol, irk, last_atom, Nlist
542 : LOGICAL :: ldum, molname_generated, my_shell
543 : REAL(KIND=dp), DIMENSION(3) :: r, s
544 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
545 : TYPE(molecule_type), POINTER :: molecule_now
546 : TYPE(section_type), POINTER :: section
547 : TYPE(val_type), POINTER :: my_val, old_val
548 :
549 20746 : CALL timeset(routineN, handle)
550 :
551 20746 : NULLIFY (my_val, old_val, section, vals)
552 20746 : my_shell = .FALSE.
553 20746 : IF (PRESENT(shell)) my_shell = shell
554 20746 : CPASSERT(ASSOCIATED(coord_section))
555 20746 : CPASSERT(coord_section%ref_count > 0)
556 20746 : section => coord_section%section
557 20746 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
558 20746 : IF (ik == -2) &
559 : CALL cp_abort(__LOCATION__, &
560 : "section "//TRIM(section%name)//" does not contain keyword "// &
561 0 : "_DEFAULT_KEYWORD_")
562 :
563 898 : DO
564 21644 : IF (SIZE(coord_section%values, 2) == 1) EXIT
565 898 : CALL section_vals_add_values(coord_section)
566 : END DO
567 20746 : vals => coord_section%values(ik, 1)%list
568 20746 : Nlist = 0
569 20746 : IF (ASSOCIATED(vals)) THEN
570 19702 : Nlist = cp_sll_val_get_length(vals)
571 : END IF
572 20746 : imol = 0
573 20746 : last_atom = 0
574 1593376 : DO irk = 1, particles%n_els
575 1572630 : CALL get_atomic_kind(particles%els(irk)%atomic_kind, name=name)
576 1593376 : IF (my_shell) THEN
577 2766288 : s = particles%els(irk)%r
578 691572 : IF (scale) THEN
579 0 : r = s
580 0 : CALL real_to_scaled(s, r, cell)
581 : ELSE
582 2766288 : s = s*conv_factor
583 : END IF
584 : WRITE (UNIT=line, FMT="(T7,A,3ES25.16E3,1X,I0)") &
585 691572 : TRIM(ADJUSTL(name)), s(1:3), particles%els(irk)%atom_index
586 691572 : CALL val_create(my_val, lc_val=line)
587 691572 : IF (Nlist /= 0) THEN
588 651672 : IF (irk == 1) THEN
589 6348 : new_pos => vals
590 : ELSE
591 645324 : new_pos => new_pos%rest
592 : END IF
593 651672 : old_val => new_pos%first_el
594 651672 : CALL val_release(old_val)
595 651672 : new_pos%first_el => my_val
596 : ELSE
597 39900 : IF (irk == 1) THEN
598 256 : NULLIFY (new_pos)
599 256 : CALL cp_sll_val_create(new_pos, first_el=my_val)
600 256 : vals => new_pos
601 : ELSE
602 39644 : NULLIFY (new_pos%rest)
603 39644 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
604 39644 : new_pos => new_pos%rest
605 : END IF
606 : END IF
607 691572 : NULLIFY (my_val)
608 : ELSE
609 881058 : IF (last_atom < irk) THEN
610 382688 : imol = imol + 1
611 382688 : molecule_now => molecules%els(imol)
612 382688 : CALL get_molecule(molecule_now, last_atom=last_atom)
613 : CALL get_molecule_kind(molecule_now%molecule_kind, molname_generated=molname_generated, &
614 382688 : name=my_tag)
615 382688 : IF (molname_generated) my_tag = ""
616 : END IF
617 881058 : ldum = qmmm_ff_precond_only_qm(my_tag)
618 881058 : ldum = qmmm_ff_precond_only_qm(name)
619 3524232 : s = particles%els(irk)%r
620 881058 : IF (scale) THEN
621 403870 : r = s
622 403870 : CALL real_to_scaled(s, r, cell)
623 : ELSE
624 1908752 : s = s*conv_factor
625 : END IF
626 881058 : IF (LEN_TRIM(my_tag) > 0) THEN
627 : WRITE (UNIT=line, FMT="(T7,A,3ES25.16E3,1X,A,1X,I0)") &
628 450594 : TRIM(ADJUSTL(name)), s(1:3), TRIM(my_tag), imol
629 : ELSE
630 : WRITE (UNIT=line, FMT="(T7,A,3ES25.16E3)") &
631 430464 : TRIM(ADJUSTL(name)), s(1:3)
632 : END IF
633 881058 : CALL val_create(my_val, lc_val=line)
634 :
635 881058 : IF (Nlist /= 0) THEN
636 679766 : IF (irk == 1) THEN
637 13354 : new_pos => vals
638 : ELSE
639 666412 : new_pos => new_pos%rest
640 : END IF
641 679766 : old_val => new_pos%first_el
642 679766 : CALL val_release(old_val)
643 679766 : new_pos%first_el => my_val
644 : ELSE
645 201292 : IF (irk == 1) THEN
646 788 : NULLIFY (new_pos)
647 788 : CALL cp_sll_val_create(new_pos, first_el=my_val)
648 788 : vals => new_pos
649 : ELSE
650 200504 : NULLIFY (new_pos%rest)
651 200504 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
652 200504 : new_pos => new_pos%rest
653 : END IF
654 : END IF
655 881058 : NULLIFY (my_val)
656 : END IF
657 : END DO
658 :
659 20746 : coord_section%values(ik, 1)%list => vals
660 :
661 20746 : CALL timestop(handle)
662 :
663 20746 : END SUBROUTINE section_coord_val_set
664 :
665 : ! **************************************************************************************************
666 : !> \brief routine to dump dipoles.. fast implementation
667 : !> \param dipoles ...
668 : !> \param dipoles_section ...
669 : !> \par History
670 : !> 12.2007 created [teo]
671 : !> \author Teodoro Laino
672 : ! **************************************************************************************************
673 24 : SUBROUTINE update_dipoles_section(dipoles, dipoles_section)
674 :
675 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dipoles
676 : TYPE(section_vals_type), POINTER :: dipoles_section
677 :
678 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_dipoles_section'
679 :
680 : INTEGER :: handle, ik, irk, Nlist, nloop
681 24 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
682 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
683 : TYPE(section_type), POINTER :: section
684 : TYPE(val_type), POINTER :: my_val, old_val
685 :
686 24 : CALL timeset(routineN, handle)
687 24 : NULLIFY (my_val, old_val, section, vals)
688 24 : CPASSERT(ASSOCIATED(dipoles_section))
689 24 : CPASSERT(dipoles_section%ref_count > 0)
690 24 : section => dipoles_section%section
691 24 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
692 24 : IF (ik == -2) &
693 : CALL cp_abort(__LOCATION__, &
694 : "section "//TRIM(section%name)//" does not contain keyword "// &
695 0 : "_DEFAULT_KEYWORD_")
696 :
697 : ! At least one of the two arguments must be present..
698 24 : nloop = SIZE(dipoles, 2)
699 16 : DO
700 40 : IF (SIZE(dipoles_section%values, 2) == 1) EXIT
701 16 : CALL section_vals_add_values(dipoles_section)
702 : END DO
703 24 : vals => dipoles_section%values(ik, 1)%list
704 24 : Nlist = 0
705 24 : IF (ASSOCIATED(vals)) THEN
706 8 : Nlist = cp_sll_val_get_length(vals)
707 : END IF
708 4656 : DO irk = 1, nloop
709 4632 : ALLOCATE (work(3))
710 : ! Always stored in A.U.
711 37056 : work = dipoles(1:3, irk)
712 4632 : CALL val_create(my_val, r_vals_ptr=work)
713 :
714 4632 : IF (Nlist /= 0) THEN
715 2154 : IF (irk == 1) THEN
716 8 : new_pos => vals
717 : ELSE
718 2146 : new_pos => new_pos%rest
719 : END IF
720 2154 : old_val => new_pos%first_el
721 2154 : CALL val_release(old_val)
722 2154 : new_pos%first_el => my_val
723 : ELSE
724 2478 : IF (irk == 1) THEN
725 16 : NULLIFY (new_pos)
726 16 : CALL cp_sll_val_create(new_pos, first_el=my_val)
727 16 : vals => new_pos
728 : ELSE
729 2462 : NULLIFY (new_pos%rest)
730 2462 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
731 2462 : new_pos => new_pos%rest
732 : END IF
733 : END IF
734 4656 : NULLIFY (my_val)
735 : END DO
736 24 : dipoles_section%values(ik, 1)%list => vals
737 :
738 24 : CALL timestop(handle)
739 :
740 24 : END SUBROUTINE update_dipoles_section
741 :
742 : ! **************************************************************************************************
743 : !> \brief routine to dump quadrupoles.. fast implementation
744 : !> \param quadrupoles ...
745 : !> \param quadrupoles_section ...
746 : !> \par History
747 : !> 12.2007 created [teo]
748 : !> \author Teodoro Laino
749 : ! **************************************************************************************************
750 0 : SUBROUTINE update_quadrupoles_section(quadrupoles, quadrupoles_section)
751 :
752 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: quadrupoles
753 : TYPE(section_vals_type), POINTER :: quadrupoles_section
754 :
755 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_quadrupoles_section'
756 :
757 : INTEGER :: handle, i, ik, ind, irk, j, Nlist, nloop
758 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
759 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
760 : TYPE(section_type), POINTER :: section
761 : TYPE(val_type), POINTER :: my_val, old_val
762 :
763 0 : CALL timeset(routineN, handle)
764 0 : NULLIFY (my_val, old_val, section, vals)
765 0 : CPASSERT(ASSOCIATED(quadrupoles_section))
766 0 : CPASSERT(quadrupoles_section%ref_count > 0)
767 0 : section => quadrupoles_section%section
768 0 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
769 0 : IF (ik == -2) &
770 : CALL cp_abort(__LOCATION__, &
771 : "section "//TRIM(section%name)//" does not contain keyword "// &
772 0 : "_DEFAULT_KEYWORD_")
773 :
774 : ! At least one of the two arguments must be present..
775 0 : nloop = SIZE(quadrupoles, 2)
776 0 : DO
777 0 : IF (SIZE(quadrupoles_section%values, 2) == 1) EXIT
778 0 : CALL section_vals_add_values(quadrupoles_section)
779 : END DO
780 0 : vals => quadrupoles_section%values(ik, 1)%list
781 0 : Nlist = 0
782 0 : IF (ASSOCIATED(vals)) THEN
783 0 : Nlist = cp_sll_val_get_length(vals)
784 : END IF
785 0 : DO irk = 1, nloop
786 0 : ALLOCATE (work(6))
787 : ! Always stored in A.U.
788 0 : ind = 0
789 0 : DO i = 1, 3
790 0 : DO j = i, 3
791 0 : ind = ind + 1
792 0 : work(ind) = quadrupoles(j, i, irk)
793 : END DO
794 : END DO
795 0 : CALL val_create(my_val, r_vals_ptr=work)
796 :
797 0 : IF (Nlist /= 0) THEN
798 0 : IF (irk == 1) THEN
799 0 : new_pos => vals
800 : ELSE
801 0 : new_pos => new_pos%rest
802 : END IF
803 0 : old_val => new_pos%first_el
804 0 : CALL val_release(old_val)
805 0 : new_pos%first_el => my_val
806 : ELSE
807 0 : IF (irk == 1) THEN
808 0 : NULLIFY (new_pos)
809 0 : CALL cp_sll_val_create(new_pos, first_el=my_val)
810 0 : vals => new_pos
811 : ELSE
812 0 : NULLIFY (new_pos%rest)
813 0 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
814 0 : new_pos => new_pos%rest
815 : END IF
816 : END IF
817 0 : NULLIFY (my_val)
818 : END DO
819 :
820 0 : quadrupoles_section%values(ik, 1)%list => vals
821 :
822 0 : CALL timestop(handle)
823 :
824 0 : END SUBROUTINE update_quadrupoles_section
825 :
826 : ! **************************************************************************************************
827 : !> \brief Dump atomic, core, or shell coordinates to a file in CP2K &COORD
828 : !> section format
829 : !> \param particles ...
830 : !> \param molecules ...
831 : !> \param cell ...
832 : !> \param conv_factor ...
833 : !> \param output_unit ...
834 : !> \param core_or_shell ...
835 : !> \param scaled_coordinates ...
836 : !> \par History
837 : !> 07.02.2011 (Creation, MK)
838 : !> \author Matthias Krack (MK)
839 : !> \version 1.0
840 : ! **************************************************************************************************
841 0 : SUBROUTINE dump_coordinates_cp2k(particles, molecules, cell, conv_factor, &
842 : output_unit, core_or_shell, &
843 : scaled_coordinates)
844 :
845 : TYPE(particle_list_type), POINTER :: particles
846 : TYPE(molecule_list_type), POINTER :: molecules
847 : TYPE(cell_type), POINTER :: cell
848 : REAL(KIND=dp), INTENT(IN) :: conv_factor
849 : INTEGER, INTENT(IN) :: output_unit
850 : LOGICAL, INTENT(IN) :: core_or_shell, scaled_coordinates
851 :
852 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_coordinates_cp2k'
853 :
854 : CHARACTER(LEN=default_string_length) :: kind_name, tag
855 : INTEGER :: handle, imolecule, iparticle, last_atom
856 : LOGICAL :: ldum, molname_generated
857 : REAL(KIND=dp), DIMENSION(3) :: r, s
858 : TYPE(molecule_type), POINTER :: molecule
859 :
860 0 : CALL timeset(routineN, handle)
861 :
862 0 : CPASSERT(ASSOCIATED(particles))
863 0 : IF (.NOT. core_or_shell) THEN
864 0 : CPASSERT(ASSOCIATED(molecules))
865 : END IF
866 0 : IF (scaled_coordinates) THEN
867 0 : CPASSERT(ASSOCIATED(cell))
868 : END IF
869 :
870 0 : kind_name = ""
871 0 : tag = ""
872 0 : imolecule = 0
873 0 : last_atom = 0
874 0 : DO iparticle = 1, particles%n_els
875 0 : CALL get_atomic_kind(particles%els(iparticle)%atomic_kind, name=kind_name)
876 0 : IF (.NOT. core_or_shell) THEN
877 0 : IF (iparticle > last_atom) THEN
878 0 : imolecule = imolecule + 1
879 0 : molecule => molecules%els(imolecule)
880 0 : CALL get_molecule(molecule, last_atom=last_atom)
881 : CALL get_molecule_kind(molecule%molecule_kind, &
882 : molname_generated=molname_generated, &
883 0 : name=tag)
884 0 : IF (molname_generated) tag = ""
885 : END IF
886 0 : ldum = qmmm_ff_precond_only_qm(tag)
887 0 : ldum = qmmm_ff_precond_only_qm(kind_name)
888 : END IF
889 0 : IF (scaled_coordinates) THEN
890 0 : CALL real_to_scaled(s, particles%els(iparticle)%r, cell)
891 0 : r(1:3) = s(1:3)
892 : ELSE
893 0 : r(1:3) = particles%els(iparticle)%r(1:3)*conv_factor
894 : END IF
895 0 : IF (core_or_shell) THEN
896 : WRITE (UNIT=output_unit, FMT="(A,3ES25.16E3,1X,I0)") &
897 0 : TRIM(ADJUSTL(kind_name)), r(1:3), particles%els(iparticle)%atom_index
898 : ELSE
899 0 : IF (LEN_TRIM(tag) > 0) THEN
900 : WRITE (UNIT=output_unit, FMT="(A,3ES25.16E3,1X,A,1X,I0)") &
901 0 : TRIM(ADJUSTL(kind_name)), r(1:3), TRIM(tag), imolecule
902 : ELSE
903 : WRITE (UNIT=output_unit, FMT="(A,3ES25.16E3)") &
904 0 : TRIM(ADJUSTL(kind_name)), r(1:3)
905 : END IF
906 : END IF
907 : END DO
908 :
909 0 : CALL timestop(handle)
910 :
911 0 : END SUBROUTINE dump_coordinates_cp2k
912 :
913 : END MODULE input_restart_force_eval
|
