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 : 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 14200 : 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 14200 : INTEGER, DIMENSION(:), POINTER :: i_force_eval
95 : LOGICAL :: is_present, multiple_subsys, &
96 : skip_vel_section
97 14200 : 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 14200 : 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 14200 : 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 14200 : (myid /= ehrenfest))
114 :
115 : ! Go on updatig the force_env_section
116 14200 : force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
117 14200 : 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 14200 : iforce_eval = 1
121 : subsys_section => section_vals_get_subs_vals3(force_env_sections, "SUBSYS", &
122 14200 : i_rep_section=i_force_eval(iforce_eval))
123 : CALL update_subsys(subsys_section, force_env, skip_vel_section, &
124 14200 : write_binary_restart_file)
125 :
126 : ! For RESPA we need to update all subsystems
127 14200 : IF (PRESENT(respa)) THEN
128 14178 : 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 14200 : rng_section => section_vals_get_subs_vals(subsys_section, "RNG_INIT")
140 14200 : CALL update_rng_particle(rng_section, force_env)
141 :
142 : qmmm_section => section_vals_get_subs_vals3(force_env_sections, "QMMM", &
143 14200 : i_rep_section=i_force_eval(iforce_eval))
144 14200 : CALL update_qmmm(qmmm_section, force_env)
145 :
146 : ! Non-mixed CDFT calculation: update constraint Lagrangian and counter
147 14200 : IF (nforce_eval == 1) THEN
148 14136 : IF (ASSOCIATED(force_env%qs_env)) THEN
149 2356 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
150 2356 : 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 14200 : 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 14200 : IF (ASSOCIATED(force_env%qs_env)) THEN
194 2360 : 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 2360 : i_rep_section=i_force_eval(iforce_eval))
197 2360 : tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION")
198 2360 : CALL section_vals_get(tmp_section, explicit=is_present)
199 2360 : IF (is_present) THEN
200 : CALL section_vals_val_set(root_section, "FORCE_EVAL%DFT%REAL_TIME_PROPAGATION%INITIAL_WFN", &
201 156 : i_val=use_rt_restart)
202 156 : 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 14200 : DEALLOCATE (i_force_eval)
214 :
215 14200 : 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 14200 : 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 14200 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
232 :
233 14454 : 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 14200 : 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 14200 : 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 14200 : 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 14200 : 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 14200 : particles=particles)
273 :
274 14200 : 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 14200 : 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 14238 : 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 14238 : 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 14238 : NULLIFY (work_section, core_particles, particles, shell_particles, &
338 14238 : subsys, cell, para_env, multipoles)
339 14238 : CALL timeset(routineN, handle)
340 14238 : 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 14238 : 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 14238 : ALLOCATE (multiple_unit_cell(3))
349 56952 : multiple_unit_cell = 1
350 : CALL section_vals_val_set(subsys_section, "TOPOLOGY%MULTIPLE_UNIT_CELL", &
351 14238 : i_vals_ptr=multiple_unit_cell)
352 :
353 : ! Coordinates and Velocities
354 14238 : work_section => section_vals_get_subs_vals(subsys_section, "COORD")
355 14238 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
356 14238 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
357 14238 : conv_factor = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
358 14238 : IF (.NOT. write_binary_restart_file) THEN
359 : CALL section_vals_val_get(subsys_section, "TOPOLOGY%COORD_FILE_FORMAT", &
360 14102 : i_val=coord_file_format)
361 14102 : 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 : CALL section_coord_val_set(work_section, particles, molecules, conv_factor, scale, &
380 14102 : cell)
381 : END IF
382 : END IF
383 14238 : CALL section_vals_val_set(subsys_section, "TOPOLOGY%NUMBER_OF_ATOMS", i_val=particles%n_els)
384 14238 : work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
385 14238 : IF (.NOT. skip_vel_section) THEN
386 5570 : IF (.NOT. write_binary_restart_file) THEN
387 5434 : CALL section_velocity_val_set(work_section, particles, conv_factor=1.0_dp)
388 : END IF
389 : ELSE
390 8668 : CALL section_vals_remove_values(work_section)
391 : END IF
392 :
393 : ! Write restart input for core-shell model
394 14238 : IF (.NOT. write_binary_restart_file) THEN
395 14102 : IF (ASSOCIATED(shell_particles)) THEN
396 3302 : work_section => section_vals_get_subs_vals(subsys_section, "SHELL_COORD")
397 3302 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
398 3302 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
399 3302 : conv_factor = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
400 : CALL section_coord_val_set(work_section, shell_particles, molecules, &
401 3302 : conv_factor, scale, cell, shell=.TRUE.)
402 3302 : IF (.NOT. skip_vel_section) THEN
403 362 : work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
404 362 : CALL section_velocity_val_set(work_section, shell_particles, conv_factor=1.0_dp)
405 : END IF
406 : END IF
407 14102 : IF (ASSOCIATED(core_particles)) THEN
408 3302 : work_section => section_vals_get_subs_vals(subsys_section, "CORE_COORD")
409 3302 : CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str)
410 3302 : CALL section_vals_val_get(work_section, "SCALED", l_val=scale)
411 3302 : conv_factor = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
412 : CALL section_coord_val_set(work_section, core_particles, molecules, &
413 3302 : conv_factor, scale, cell, shell=.TRUE.)
414 3302 : IF (.NOT. skip_vel_section) THEN
415 362 : work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
416 362 : CALL section_velocity_val_set(work_section, core_particles, conv_factor=1.0_dp)
417 : END IF
418 : END IF
419 : END IF
420 :
421 : ! Updating cell info
422 14238 : CALL force_env_get(force_env, cell=cell)
423 14238 : work_section => section_vals_get_subs_vals(subsys_section, "CELL")
424 14238 : CALL update_cell_section(cell, cell_section=work_section)
425 :
426 : ! Updating cell_ref info
427 14238 : use_ref_cell = .FALSE.
428 16606 : SELECT CASE (force_env%in_use)
429 : CASE (use_qs_force)
430 2368 : CALL get_qs_env(force_env%qs_env, cell_ref=cell, use_ref_cell=use_ref_cell)
431 : CASE (use_eip_force)
432 16 : CALL eip_env_get(force_env%eip_env, cell_ref=cell, use_ref_cell=use_ref_cell)
433 : CASE (use_nnp_force)
434 14238 : CALL nnp_env_get(force_env%nnp_env, cell_ref=cell, use_ref_cell=use_ref_cell)
435 : END SELECT
436 14238 : IF (use_ref_cell) THEN
437 34 : work_section => section_vals_get_subs_vals(subsys_section, "CELL%CELL_REF")
438 34 : CALL update_cell_section(cell, cell_section=work_section)
439 : END IF
440 :
441 : ! Updating multipoles
442 14238 : IF (ASSOCIATED(multipoles)) THEN
443 24 : work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES")
444 16 : DO
445 40 : IF (SIZE(work_section%values, 2) == 1) EXIT
446 16 : CALL section_vals_add_values(work_section)
447 : END DO
448 24 : IF (ASSOCIATED(multipoles%dipoles)) THEN
449 24 : work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES%DIPOLES")
450 24 : CALL update_dipoles_section(multipoles%dipoles, work_section)
451 : END IF
452 24 : IF (ASSOCIATED(multipoles%quadrupoles)) THEN
453 0 : work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES%QUADRUPOLES")
454 0 : CALL update_quadrupoles_section(multipoles%quadrupoles, work_section)
455 : END IF
456 : END IF
457 :
458 14238 : CALL timestop(handle)
459 :
460 28476 : END SUBROUTINE update_subsys
461 :
462 : ! **************************************************************************************************
463 : !> \brief Routine to update a cell section
464 : !> \param cell ...
465 : !> \param cell_section ...
466 : !> \author Ole Schuett
467 : ! **************************************************************************************************
468 14292 : SUBROUTINE update_cell_section(cell, cell_section)
469 :
470 : TYPE(cell_type), POINTER :: cell
471 : TYPE(section_vals_type), POINTER :: cell_section
472 :
473 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_cell_section'
474 :
475 : INTEGER :: handle
476 14292 : INTEGER, DIMENSION(:), POINTER :: iwork
477 14292 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
478 :
479 14292 : NULLIFY (work, iwork)
480 14292 : CALL timeset(routineN, handle)
481 :
482 : ! CELL VECTORS - A
483 14292 : ALLOCATE (work(3))
484 114336 : work(1:3) = cell%hmat(1:3, 1)
485 14292 : CALL section_vals_val_set(cell_section, "A", r_vals_ptr=work)
486 :
487 : ! CELL VECTORS - B
488 14292 : ALLOCATE (work(3))
489 114336 : work(1:3) = cell%hmat(1:3, 2)
490 14292 : CALL section_vals_val_set(cell_section, "B", r_vals_ptr=work)
491 :
492 : ! CELL VECTORS - C
493 14292 : ALLOCATE (work(3))
494 114336 : work(1:3) = cell%hmat(1:3, 3)
495 14292 : CALL section_vals_val_set(cell_section, "C", r_vals_ptr=work)
496 :
497 : ! MULTIPLE_UNIT_CELL
498 14292 : ALLOCATE (iwork(3))
499 57168 : iwork = 1
500 14292 : CALL section_vals_val_set(cell_section, "MULTIPLE_UNIT_CELL", i_vals_ptr=iwork)
501 :
502 : ! Unset unused or misleading information
503 14292 : CALL section_vals_val_unset(cell_section, "ABC")
504 14292 : CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA")
505 :
506 14292 : CALL timestop(handle)
507 :
508 14292 : END SUBROUTINE update_cell_section
509 :
510 : ! **************************************************************************************************
511 : !> \brief routine to dump coordinates.. fast implementation
512 : !> \param coord_section ...
513 : !> \param particles ...
514 : !> \param molecules ...
515 : !> \param conv_factor ...
516 : !> \param scale ...
517 : !> \param cell ...
518 : !> \param shell ...
519 : !> \par History
520 : !> 02.2006 created [teo]
521 : !> \author Teodoro Laino
522 : ! **************************************************************************************************
523 20706 : SUBROUTINE section_coord_val_set(coord_section, particles, molecules, conv_factor, &
524 : scale, cell, shell)
525 :
526 : TYPE(section_vals_type), POINTER :: coord_section
527 : TYPE(particle_list_type), POINTER :: particles
528 : TYPE(molecule_list_type), POINTER :: molecules
529 : REAL(KIND=dp) :: conv_factor
530 : LOGICAL, INTENT(IN) :: scale
531 : TYPE(cell_type), POINTER :: cell
532 : LOGICAL, INTENT(IN), OPTIONAL :: shell
533 :
534 : CHARACTER(LEN=*), PARAMETER :: routineN = 'section_coord_val_set'
535 :
536 : CHARACTER(LEN=2*default_string_length) :: line
537 : CHARACTER(LEN=default_string_length) :: my_tag, name
538 : INTEGER :: handle, ik, imol, irk, last_atom, Nlist
539 : LOGICAL :: ldum, molname_generated, my_shell
540 : REAL(KIND=dp), DIMENSION(3) :: r, s
541 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
542 : TYPE(molecule_type), POINTER :: molecule_now
543 : TYPE(section_type), POINTER :: section
544 : TYPE(val_type), POINTER :: my_val, old_val
545 :
546 20706 : CALL timeset(routineN, handle)
547 :
548 20706 : NULLIFY (my_val, old_val, section, vals)
549 20706 : my_shell = .FALSE.
550 20706 : IF (PRESENT(shell)) my_shell = shell
551 20706 : CPASSERT(ASSOCIATED(coord_section))
552 20706 : CPASSERT(coord_section%ref_count > 0)
553 20706 : section => coord_section%section
554 20706 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
555 20706 : IF (ik == -2) &
556 : CALL cp_abort(__LOCATION__, &
557 : "section "//TRIM(section%name)//" does not contain keyword "// &
558 0 : "_DEFAULT_KEYWORD_")
559 :
560 896 : DO
561 21602 : IF (SIZE(coord_section%values, 2) == 1) EXIT
562 896 : CALL section_vals_add_values(coord_section)
563 : END DO
564 20706 : vals => coord_section%values(ik, 1)%list
565 20706 : Nlist = 0
566 20706 : IF (ASSOCIATED(vals)) THEN
567 19670 : Nlist = cp_sll_val_get_length(vals)
568 : END IF
569 20706 : imol = 0
570 20706 : last_atom = 0
571 1587812 : DO irk = 1, particles%n_els
572 1567106 : CALL get_atomic_kind(particles%els(irk)%atomic_kind, name=name)
573 1587812 : IF (my_shell) THEN
574 2766288 : s = particles%els(irk)%r
575 691572 : IF (scale) THEN
576 0 : r = s
577 0 : CALL real_to_scaled(s, r, cell)
578 : ELSE
579 2766288 : s = s*conv_factor
580 : END IF
581 : WRITE (UNIT=line, FMT="(T7,A,3ES25.16E3,1X,I0)") &
582 691572 : TRIM(ADJUSTL(name)), s(1:3), particles%els(irk)%atom_index
583 691572 : CALL val_create(my_val, lc_val=line)
584 691572 : IF (Nlist /= 0) THEN
585 651672 : IF (irk == 1) THEN
586 6348 : new_pos => vals
587 : ELSE
588 645324 : new_pos => new_pos%rest
589 : END IF
590 651672 : old_val => new_pos%first_el
591 651672 : CALL val_release(old_val)
592 651672 : new_pos%first_el => my_val
593 : ELSE
594 39900 : IF (irk == 1) THEN
595 256 : NULLIFY (new_pos)
596 256 : CALL cp_sll_val_create(new_pos, first_el=my_val)
597 256 : vals => new_pos
598 : ELSE
599 39644 : NULLIFY (new_pos%rest)
600 39644 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
601 39644 : new_pos => new_pos%rest
602 : END IF
603 : END IF
604 691572 : NULLIFY (my_val)
605 : ELSE
606 875534 : IF (last_atom < irk) THEN
607 377498 : imol = imol + 1
608 377498 : molecule_now => molecules%els(imol)
609 377498 : CALL get_molecule(molecule_now, last_atom=last_atom)
610 : CALL get_molecule_kind(molecule_now%molecule_kind, molname_generated=molname_generated, &
611 377498 : name=my_tag)
612 377498 : IF (molname_generated) my_tag = ""
613 : END IF
614 875534 : ldum = qmmm_ff_precond_only_qm(my_tag)
615 875534 : ldum = qmmm_ff_precond_only_qm(name)
616 3502136 : s = particles%els(irk)%r
617 875534 : IF (scale) THEN
618 400710 : r = s
619 400710 : CALL real_to_scaled(s, r, cell)
620 : ELSE
621 1899296 : s = s*conv_factor
622 : END IF
623 875534 : IF (LEN_TRIM(my_tag) > 0) THEN
624 : WRITE (UNIT=line, FMT="(T7,A,3ES25.16E3,1X,A,1X,I0)") &
625 447594 : TRIM(ADJUSTL(name)), s(1:3), TRIM(my_tag), imol
626 : ELSE
627 : WRITE (UNIT=line, FMT="(T7,A,3ES25.16E3)") &
628 427940 : TRIM(ADJUSTL(name)), s(1:3)
629 : END IF
630 875534 : CALL val_create(my_val, lc_val=line)
631 :
632 875534 : IF (Nlist /= 0) THEN
633 676546 : IF (irk == 1) THEN
634 13322 : new_pos => vals
635 : ELSE
636 663224 : new_pos => new_pos%rest
637 : END IF
638 676546 : old_val => new_pos%first_el
639 676546 : CALL val_release(old_val)
640 676546 : new_pos%first_el => my_val
641 : ELSE
642 198988 : IF (irk == 1) THEN
643 780 : NULLIFY (new_pos)
644 780 : CALL cp_sll_val_create(new_pos, first_el=my_val)
645 780 : vals => new_pos
646 : ELSE
647 198208 : NULLIFY (new_pos%rest)
648 198208 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
649 198208 : new_pos => new_pos%rest
650 : END IF
651 : END IF
652 875534 : NULLIFY (my_val)
653 : END IF
654 : END DO
655 :
656 20706 : coord_section%values(ik, 1)%list => vals
657 :
658 20706 : CALL timestop(handle)
659 :
660 20706 : END SUBROUTINE section_coord_val_set
661 :
662 : ! **************************************************************************************************
663 : !> \brief routine to dump dipoles.. fast implementation
664 : !> \param dipoles ...
665 : !> \param dipoles_section ...
666 : !> \par History
667 : !> 12.2007 created [teo]
668 : !> \author Teodoro Laino
669 : ! **************************************************************************************************
670 24 : SUBROUTINE update_dipoles_section(dipoles, dipoles_section)
671 :
672 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dipoles
673 : TYPE(section_vals_type), POINTER :: dipoles_section
674 :
675 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_dipoles_section'
676 :
677 : INTEGER :: handle, ik, irk, Nlist, nloop
678 24 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
679 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
680 : TYPE(section_type), POINTER :: section
681 : TYPE(val_type), POINTER :: my_val, old_val
682 :
683 24 : CALL timeset(routineN, handle)
684 24 : NULLIFY (my_val, old_val, section, vals)
685 24 : CPASSERT(ASSOCIATED(dipoles_section))
686 24 : CPASSERT(dipoles_section%ref_count > 0)
687 24 : section => dipoles_section%section
688 24 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
689 24 : IF (ik == -2) &
690 : CALL cp_abort(__LOCATION__, &
691 : "section "//TRIM(section%name)//" does not contain keyword "// &
692 0 : "_DEFAULT_KEYWORD_")
693 :
694 : ! At least one of the two arguments must be present..
695 24 : nloop = SIZE(dipoles, 2)
696 16 : DO
697 40 : IF (SIZE(dipoles_section%values, 2) == 1) EXIT
698 16 : CALL section_vals_add_values(dipoles_section)
699 : END DO
700 24 : vals => dipoles_section%values(ik, 1)%list
701 24 : Nlist = 0
702 24 : IF (ASSOCIATED(vals)) THEN
703 8 : Nlist = cp_sll_val_get_length(vals)
704 : END IF
705 4656 : DO irk = 1, nloop
706 4632 : ALLOCATE (work(3))
707 : ! Always stored in A.U.
708 37056 : work = dipoles(1:3, irk)
709 4632 : CALL val_create(my_val, r_vals_ptr=work)
710 :
711 4632 : IF (Nlist /= 0) THEN
712 2154 : IF (irk == 1) THEN
713 8 : new_pos => vals
714 : ELSE
715 2146 : new_pos => new_pos%rest
716 : END IF
717 2154 : old_val => new_pos%first_el
718 2154 : CALL val_release(old_val)
719 2154 : new_pos%first_el => my_val
720 : ELSE
721 2478 : IF (irk == 1) THEN
722 16 : NULLIFY (new_pos)
723 16 : CALL cp_sll_val_create(new_pos, first_el=my_val)
724 16 : vals => new_pos
725 : ELSE
726 2462 : NULLIFY (new_pos%rest)
727 2462 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
728 2462 : new_pos => new_pos%rest
729 : END IF
730 : END IF
731 4656 : NULLIFY (my_val)
732 : END DO
733 24 : dipoles_section%values(ik, 1)%list => vals
734 :
735 24 : CALL timestop(handle)
736 :
737 24 : END SUBROUTINE update_dipoles_section
738 :
739 : ! **************************************************************************************************
740 : !> \brief routine to dump quadrupoles.. fast implementation
741 : !> \param quadrupoles ...
742 : !> \param quadrupoles_section ...
743 : !> \par History
744 : !> 12.2007 created [teo]
745 : !> \author Teodoro Laino
746 : ! **************************************************************************************************
747 0 : SUBROUTINE update_quadrupoles_section(quadrupoles, quadrupoles_section)
748 :
749 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: quadrupoles
750 : TYPE(section_vals_type), POINTER :: quadrupoles_section
751 :
752 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_quadrupoles_section'
753 :
754 : INTEGER :: handle, i, ik, ind, irk, j, Nlist, nloop
755 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
756 : TYPE(cp_sll_val_type), POINTER :: new_pos, vals
757 : TYPE(section_type), POINTER :: section
758 : TYPE(val_type), POINTER :: my_val, old_val
759 :
760 0 : CALL timeset(routineN, handle)
761 0 : NULLIFY (my_val, old_val, section, vals)
762 0 : CPASSERT(ASSOCIATED(quadrupoles_section))
763 0 : CPASSERT(quadrupoles_section%ref_count > 0)
764 0 : section => quadrupoles_section%section
765 0 : ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_")
766 0 : IF (ik == -2) &
767 : CALL cp_abort(__LOCATION__, &
768 : "section "//TRIM(section%name)//" does not contain keyword "// &
769 0 : "_DEFAULT_KEYWORD_")
770 :
771 : ! At least one of the two arguments must be present..
772 0 : nloop = SIZE(quadrupoles, 2)
773 0 : DO
774 0 : IF (SIZE(quadrupoles_section%values, 2) == 1) EXIT
775 0 : CALL section_vals_add_values(quadrupoles_section)
776 : END DO
777 0 : vals => quadrupoles_section%values(ik, 1)%list
778 0 : Nlist = 0
779 0 : IF (ASSOCIATED(vals)) THEN
780 0 : Nlist = cp_sll_val_get_length(vals)
781 : END IF
782 0 : DO irk = 1, nloop
783 0 : ALLOCATE (work(6))
784 : ! Always stored in A.U.
785 0 : ind = 0
786 0 : DO i = 1, 3
787 0 : DO j = i, 3
788 0 : ind = ind + 1
789 0 : work(ind) = quadrupoles(j, i, irk)
790 : END DO
791 : END DO
792 0 : CALL val_create(my_val, r_vals_ptr=work)
793 :
794 0 : IF (Nlist /= 0) THEN
795 0 : IF (irk == 1) THEN
796 0 : new_pos => vals
797 : ELSE
798 0 : new_pos => new_pos%rest
799 : END IF
800 0 : old_val => new_pos%first_el
801 0 : CALL val_release(old_val)
802 0 : new_pos%first_el => my_val
803 : ELSE
804 0 : IF (irk == 1) THEN
805 0 : NULLIFY (new_pos)
806 0 : CALL cp_sll_val_create(new_pos, first_el=my_val)
807 0 : vals => new_pos
808 : ELSE
809 0 : NULLIFY (new_pos%rest)
810 0 : CALL cp_sll_val_create(new_pos%rest, first_el=my_val)
811 0 : new_pos => new_pos%rest
812 : END IF
813 : END IF
814 0 : NULLIFY (my_val)
815 : END DO
816 :
817 0 : quadrupoles_section%values(ik, 1)%list => vals
818 :
819 0 : CALL timestop(handle)
820 :
821 0 : END SUBROUTINE update_quadrupoles_section
822 :
823 : ! **************************************************************************************************
824 : !> \brief Dump atomic, core, or shell coordinates to a file in CP2K &COORD
825 : !> section format
826 : !> \param particles ...
827 : !> \param molecules ...
828 : !> \param cell ...
829 : !> \param conv_factor ...
830 : !> \param output_unit ...
831 : !> \param core_or_shell ...
832 : !> \param scaled_coordinates ...
833 : !> \par History
834 : !> 07.02.2011 (Creation, MK)
835 : !> \author Matthias Krack (MK)
836 : !> \version 1.0
837 : ! **************************************************************************************************
838 0 : SUBROUTINE dump_coordinates_cp2k(particles, molecules, cell, conv_factor, &
839 : output_unit, core_or_shell, &
840 : scaled_coordinates)
841 :
842 : TYPE(particle_list_type), POINTER :: particles
843 : TYPE(molecule_list_type), POINTER :: molecules
844 : TYPE(cell_type), POINTER :: cell
845 : REAL(KIND=dp), INTENT(IN) :: conv_factor
846 : INTEGER, INTENT(IN) :: output_unit
847 : LOGICAL, INTENT(IN) :: core_or_shell, scaled_coordinates
848 :
849 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_coordinates_cp2k'
850 :
851 : CHARACTER(LEN=default_string_length) :: kind_name, tag
852 : INTEGER :: handle, imolecule, iparticle, last_atom
853 : LOGICAL :: ldum, molname_generated
854 : REAL(KIND=dp), DIMENSION(3) :: r, s
855 : TYPE(molecule_type), POINTER :: molecule
856 :
857 0 : CALL timeset(routineN, handle)
858 :
859 0 : CPASSERT(ASSOCIATED(particles))
860 0 : IF (.NOT. core_or_shell) THEN
861 0 : CPASSERT(ASSOCIATED(molecules))
862 : END IF
863 0 : IF (scaled_coordinates) THEN
864 0 : CPASSERT(ASSOCIATED(cell))
865 : END IF
866 :
867 0 : kind_name = ""
868 0 : tag = ""
869 0 : imolecule = 0
870 0 : last_atom = 0
871 0 : DO iparticle = 1, particles%n_els
872 0 : CALL get_atomic_kind(particles%els(iparticle)%atomic_kind, name=kind_name)
873 0 : IF (.NOT. core_or_shell) THEN
874 0 : IF (iparticle > last_atom) THEN
875 0 : imolecule = imolecule + 1
876 0 : molecule => molecules%els(imolecule)
877 0 : CALL get_molecule(molecule, last_atom=last_atom)
878 : CALL get_molecule_kind(molecule%molecule_kind, &
879 : molname_generated=molname_generated, &
880 0 : name=tag)
881 0 : IF (molname_generated) tag = ""
882 : END IF
883 0 : ldum = qmmm_ff_precond_only_qm(tag)
884 0 : ldum = qmmm_ff_precond_only_qm(kind_name)
885 : END IF
886 0 : IF (scaled_coordinates) THEN
887 0 : CALL real_to_scaled(s, particles%els(iparticle)%r, cell)
888 0 : r(1:3) = s(1:3)
889 : ELSE
890 0 : r(1:3) = particles%els(iparticle)%r(1:3)*conv_factor
891 : END IF
892 0 : IF (core_or_shell) THEN
893 : WRITE (UNIT=output_unit, FMT="(A,3ES25.16E3,1X,I0)") &
894 0 : TRIM(ADJUSTL(kind_name)), r(1:3), particles%els(iparticle)%atom_index
895 : ELSE
896 0 : IF (LEN_TRIM(tag) > 0) THEN
897 : WRITE (UNIT=output_unit, FMT="(A,3ES25.16E3,1X,A,1X,I0)") &
898 0 : TRIM(ADJUSTL(kind_name)), r(1:3), TRIM(tag), imolecule
899 : ELSE
900 : WRITE (UNIT=output_unit, FMT="(A,3ES25.16E3)") &
901 0 : TRIM(ADJUSTL(kind_name)), r(1:3)
902 : END IF
903 : END IF
904 : END DO
905 :
906 0 : CALL timestop(handle)
907 :
908 0 : END SUBROUTINE dump_coordinates_cp2k
909 :
910 : END MODULE input_restart_force_eval
|
