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