(git:c641773)
Loading...
Searching...
No Matches
force_env_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
9!> \brief Interface for the force calculations
10!> \par History
11!> cjm, FEB-20-2001: pass variable box_ref
12!> cjm, SEPT-12-2002: major reorganization
13!> fawzi, APR-12-2003: introduced force_env (based on the work by CJM&JGH)
14!> fawzi, NOV-3-2004: reorganized interface for f77 interface
15!> \author fawzi
16! **************************************************************************************************
17MODULE force_env_methods
18 USE atprop_types, ONLY: atprop_init,&
19 atprop_type
20 USE bibliography, ONLY: heaton_burgess2007,&
21 huang2011,&
22 cite_reference
23 USE cell_methods, ONLY: cell_create,&
24 init_cell
25 USE cell_types, ONLY: cell_clone,&
26 cell_release,&
27 cell_sym_triclinic,&
28 cell_type,&
29 real_to_scaled,&
30 scaled_to_real
31 USE constraint_fxd, ONLY: fix_atom_control
32 USE constraint_vsite, ONLY: vsite_force_control
33 USE cp_blacs_env, ONLY: cp_blacs_env_type
34 USE cp_control_types, ONLY: dft_control_type
35 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
36 USE cp_fm_types, ONLY: cp_fm_copy_general
37 USE cp_iter_types, ONLY: cp_iteration_info_copy_iter
38 USE cp_log_handling, ONLY: cp_add_default_logger,&
39 cp_get_default_logger,&
40 cp_logger_type,&
41 cp_rm_default_logger,&
42 cp_to_string
43 USE cp_output_handling, ONLY: cp_p_file,&
44 cp_print_key_finished_output,&
45 cp_print_key_should_output,&
46 cp_print_key_unit_nr,&
47 low_print_level
48 USE cp_result_methods, ONLY: cp_results_erase,&
49 cp_results_mp_bcast,&
50 get_results,&
51 test_for_result
52 USE cp_result_types, ONLY: cp_result_copy,&
53 cp_result_create,&
54 cp_result_p_type,&
55 cp_result_release,&
56 cp_result_type
57 USE cp_subsys_types, ONLY: cp_subsys_get,&
58 cp_subsys_p_type,&
59 cp_subsys_set,&
60 cp_subsys_type
61 USE cp_units, ONLY: cp_unit_from_cp2k
62 USE eip_environment_types, ONLY: eip_environment_type
63 USE eip_silicon, ONLY: eip_bazant,&
64 eip_lenosky,&
65 eip_stillinger_weber,&
66 eip_tersoff
67 USE embed_types, ONLY: embed_env_type,&
68 opt_dmfet_pot_type,&
69 opt_embed_pot_type
70 USE external_potential_methods, ONLY: add_external_potential
71 USE fist_environment_types, ONLY: fist_environment_type
72 USE fist_force, ONLY: fist_calc_energy_force
73 USE force_env_types, ONLY: &
74 force_env_get, force_env_get_natom, force_env_p_type, force_env_set, force_env_type, &
75 use_eip_force, use_embed, use_fist_force, use_ipi, use_mixed_force, use_nnp_force, &
76 use_prog_name, use_pwdft_force, use_qmmm, use_qmmmx, use_qs_force
77 USE force_env_utils, ONLY: rescale_forces,&
78 write_atener,&
79 write_forces
80 USE fp_methods, ONLY: fp_eval
81 USE fparser, ONLY: evalerrtype,&
82 evalf,&
83 evalfd,&
84 finalizef,&
85 initf,&
86 parsef
87 USE generic_function_input, ONLY: parse_function_section
88 USE global_types, ONLY: global_environment_type,&
89 globenv_retain
90 USE grrm_utils, ONLY: write_grrm
91 USE input_constants, ONLY: &
92 cell_opt_run, debug_run, dfet, dmfet, do_method_gapw, do_method_gapw_xc, do_method_gpw, &
93 do_method_lrigpw, do_method_ofgpw, do_method_rigpw, driver_run, ehrenfest, geo_opt_run, &
94 mix_cdft, mix_coupled, mix_generic, mix_linear_combination, mix_minimum, mix_restrained, &
95 mixed_cdft_serial, mol_dyn_run, use_bazant_eip, use_lenosky_eip, use_stillinger_weber_eip, &
96 use_tersoff_eip
97 USE input_section_types, ONLY: section_vals_get_subs_vals,&
98 section_vals_retain,&
99 section_vals_type,&
100 section_vals_val_get
101 USE ipi_environment_types, ONLY: ipi_environment_type
102 USE ipi_server, ONLY: request_forces
103 USE kahan_sum, ONLY: accurate_sum
104 USE kinds, ONLY: default_path_length,&
105 default_string_length,&
106 dp
107 USE kpoint_methods, ONLY: kpoint_env_initialize,&
108 kpoint_initialize,&
109 kpoint_initialize_mos
110 USE kpoint_types, ONLY: get_kpoint_info,&
111 kpoint_reset_initialization,&
112 kpoint_sym_type,&
113 kpoint_type,&
114 set_kpoint_info
115 USE machine, ONLY: m_memory
116 USE mathlib, ONLY: abnormal_value
117 USE message_passing, ONLY: mp_para_env_type
118 USE metadynamics_types, ONLY: meta_env_type
119 USE mixed_cdft_methods, ONLY: mixed_cdft_build_weight,&
120 mixed_cdft_calculate_coupling,&
121 mixed_cdft_init
122 USE mixed_energy_types, ONLY: mixed_energy_type,&
123 mixed_force_type
124 USE mixed_environment_types, ONLY: get_mixed_env,&
125 mixed_environment_type
126 USE mixed_environment_utils, ONLY: get_subsys_map_index,&
127 mixed_map_forces
128 USE molecule_kind_list_types, ONLY: molecule_kind_list_type
129 USE molecule_kind_types, ONLY: get_molecule_kind,&
130 molecule_kind_type
131 USE nnp_environment_types, ONLY: nnp_type
132 USE nnp_force, ONLY: nnp_calc_energy_force
133 USE optimize_dmfet_potential, ONLY: build_full_dm,&
134 check_dmfet,&
135 prepare_dmfet_opt,&
136 release_dmfet_opt,&
137 subsys_spin
138 USE optimize_embedding_potential, ONLY: &
139 coulomb_guess, calculate_embed_pot_grad, conv_check_embed, get_max_subsys_diff, &
140 get_prev_density, init_embed_pot, make_subsys_embed_pot, opt_embed_step, &
141 prepare_embed_opt, print_emb_opt_info, print_embed_restart, print_pot_simple_grid, &
142 print_rho_diff, print_rho_spin_diff, read_embed_pot, release_opt_embed, step_control, &
143 understand_spin_states
144 USE particle_list_types, ONLY: particle_list_p_type,&
145 particle_list_type
146 USE particle_types, ONLY: particle_type
147 USE physcon, ONLY: debye
148 USE pw_env_types, ONLY: pw_env_get,&
149 pw_env_type
150 USE pw_methods, ONLY: pw_axpy,&
151 pw_copy,&
152 pw_integral_ab,&
153 pw_zero
154 USE pw_pool_types, ONLY: pw_pool_type
155 USE pw_types, ONLY: pw_r3d_rs_type
156 USE pwdft_environment, ONLY: pwdft_calc_energy_force
157 USE pwdft_environment_types, ONLY: pwdft_environment_type
158 USE qmmm_force, ONLY: qmmm_calc_energy_force
159 USE qmmm_types, ONLY: qmmm_env_type
160 USE qmmm_util, ONLY: apply_qmmm_translate
161 USE qmmmx_force, ONLY: qmmmx_calc_energy_force
162 USE qmmmx_types, ONLY: qmmmx_env_type
163 USE qs_apt_fdiff_methods, ONLY: apt_fdiff
164 USE qs_basis_rotation_methods, ONLY: qs_basis_rotation
165 USE qs_energy_types, ONLY: qs_energy_type
166 USE qs_environment_types, ONLY: get_qs_env,&
167 qs_environment_type,&
168 set_qs_env
169 USE qs_force, ONLY: qs_calc_energy_force
170 USE qs_mo_types, ONLY: mo_set_type
171 USE qs_rho_types, ONLY: qs_rho_get,&
172 qs_rho_type
173 USE qs_wf_history_types, ONLY: qs_wf_history_type,&
174 wfi_clear
175 USE restraint, ONLY: restraint_control
176 USE scine_utils, ONLY: write_scine
177 USE string_utilities, ONLY: compress
178 USE virial_methods, ONLY: write_stress_tensor,&
179 write_stress_tensor_components
180 USE virial_types, ONLY: symmetrize_virial,&
181 virial_p_type,&
182 virial_type,&
183 zero_virial
184#include "./base/base_uses.f90"
185
186 IMPLICIT NONE
187
188 PRIVATE
189
190 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_env_methods'
191
192 PUBLIC :: force_env_create, &
193 force_env_calc_energy_force, &
194 force_env_calc_num_pressure
195
196 INTEGER, SAVE, PRIVATE :: last_force_env_id = 0
197
198CONTAINS
199
200! **************************************************************************************************
201!> \brief Interface routine for force and energy calculations
202!> \param force_env the force_env of which you want the energy and forces
203!> \param calc_force if false the forces *might* be left unchanged
204!> or be invalid, no guarantees can be given. Defaults to true
205!> \param consistent_energies Performs an additional qs_ks_update_qs_env, so
206!> that the energies are appropriate to the forces, they are in the
207!> non-selfconsistent case not consistent to each other! [08.2005, TdK]
208!> \param skip_external_control ...
209!> \param eval_energy_forces ...
210!> \param require_consistent_energy_force ...
211!> \param linres ...
212!> \param calc_stress_tensor ...
213!> \author CJM & fawzi
214! **************************************************************************************************
215 RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
216 consistent_energies, skip_external_control, eval_energy_forces, &
217 require_consistent_energy_force, linres, calc_stress_tensor)
218
219 TYPE(force_env_type), POINTER :: force_env
220 LOGICAL, INTENT(IN), OPTIONAL :: calc_force, consistent_energies, skip_external_control, &
221 eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor
222
223 REAL(kind=dp), PARAMETER :: ateps = 1.0e-6_dp
224
225 CHARACTER(LEN=default_string_length) :: unit_string
226 INTEGER :: ikind, nat, ndigits, nfixed_atoms, &
227 nfixed_atoms_total, nkind, &
228 output_unit, print_forces, print_grrm, &
229 print_scine
230 LOGICAL :: calculate_forces, calculate_stress_tensor, do_apt_fd, energy_consistency, &
231 eval_ef, linres_run, my_skip, print_components
232 REAL(kind=dp) :: checksum, e_entropy, e_gap, e_pot, &
233 fconv, sum_energy
234 REAL(kind=dp), DIMENSION(3) :: grand_total_force, total_force
235 TYPE(atprop_type), POINTER :: atprop_env
236 TYPE(cell_type), POINTER :: cell
237 TYPE(cp_logger_type), POINTER :: logger
238 TYPE(cp_subsys_type), POINTER :: subsys
239 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
240 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
241 TYPE(molecule_kind_type), POINTER :: molecule_kind
242 TYPE(particle_list_type), POINTER :: core_particles, particles, &
243 shell_particles
244 TYPE(section_vals_type), POINTER :: print_key
245 TYPE(virial_type), POINTER :: virial
246
247 NULLIFY (logger, virial, subsys, atprop_env, cell)
248 logger => cp_get_default_logger()
249 eval_ef = .true.
250 my_skip = .false.
251 calculate_forces = .true.
252 energy_consistency = .false.
253 linres_run = .false.
254 e_gap = -1.0_dp
255 e_entropy = -1.0_dp
256 unit_string = ""
257
258 IF (PRESENT(eval_energy_forces)) eval_ef = eval_energy_forces
259 IF (PRESENT(skip_external_control)) my_skip = skip_external_control
260 IF (PRESENT(calc_force)) calculate_forces = calc_force
261 IF (PRESENT(calc_stress_tensor)) THEN
262 calculate_stress_tensor = calc_stress_tensor
263 ELSE
264 calculate_stress_tensor = calculate_forces
265 END IF
266 IF (PRESENT(consistent_energies)) energy_consistency = consistent_energies
267 IF (PRESENT(linres)) linres_run = linres
268
269 cpassert(ASSOCIATED(force_env))
270 cpassert(force_env%ref_count > 0)
271 CALL force_env_get(force_env, subsys=subsys)
272 CALL force_env_set(force_env, additional_potential=0.0_dp)
273 CALL cp_subsys_get(subsys, virial=virial, atprop=atprop_env, cell=cell)
274 IF (virial%pv_availability) CALL zero_virial(virial, reset=.false.)
275
276 nat = force_env_get_natom(force_env)
277 CALL atprop_init(atprop_env, nat)
278 IF (eval_ef) THEN
279 SELECT CASE (force_env%in_use)
280 CASE (use_fist_force)
281 CALL fist_calc_energy_force(force_env%fist_env)
282 CASE (use_qs_force)
283 CALL force_env_refresh_kpoint_symmetry(force_env, fd_energy=.NOT. calculate_forces)
284 CALL qs_calc_energy_force(force_env%qs_env, calculate_forces, energy_consistency, linres_run)
285 CASE (use_pwdft_force)
286 IF (virial%pv_availability .AND. calculate_stress_tensor) THEN
287 CALL pwdft_calc_energy_force(force_env%pwdft_env, calculate_forces,.NOT. virial%pv_numer)
288 ELSE
289 CALL pwdft_calc_energy_force(force_env%pwdft_env, calculate_forces, .false.)
290 END IF
291 e_gap = force_env%pwdft_env%energy%band_gap
292 e_entropy = force_env%pwdft_env%energy%entropy
293 CASE (use_eip_force)
294 SELECT CASE (force_env%eip_env%eip_model)
295 CASE (use_lenosky_eip)
296 CALL eip_lenosky(force_env%eip_env)
297 CASE (use_bazant_eip)
298 CALL eip_bazant(force_env%eip_env)
299 CASE (use_stillinger_weber_eip)
300 CALL eip_stillinger_weber(force_env%eip_env)
301 CASE (use_tersoff_eip)
302 CALL eip_tersoff(force_env%eip_env)
303 CASE DEFAULT
304 cpabort("Unknown EIP model.")
305 END SELECT
306 CASE (use_qmmm)
307 CALL qmmm_calc_energy_force(force_env%qmmm_env, &
308 calculate_forces, energy_consistency, linres=linres_run)
309 CASE (use_qmmmx)
310 CALL qmmmx_calc_energy_force(force_env%qmmmx_env, &
311 calculate_forces, energy_consistency, linres=linres_run, &
312 require_consistent_energy_force=require_consistent_energy_force)
313 CASE (use_mixed_force)
314 CALL mixed_energy_forces(force_env, calculate_forces)
315 CASE (use_nnp_force)
316 CALL nnp_calc_energy_force(force_env%nnp_env, &
317 calculate_forces)
318 CASE (use_embed)
319 CALL embed_energy(force_env)
320 CASE (use_ipi)
321 CALL request_forces(force_env%ipi_env)
322 CASE default
323 cpabort("Unknown force environment; cannot evaluate energy or force")
324 END SELECT
325 END IF
326 ! In case it is requested, we evaluate the stress tensor numerically
327 IF (virial%pv_availability) THEN
328 IF (virial%pv_numer .AND. calculate_stress_tensor) THEN
329 ! Compute the numerical stress tensor
330 CALL force_env_calc_num_pressure(force_env)
331 ELSE
332 IF (calculate_forces) THEN
333 ! Symmetrize analytical stress tensor
334 CALL symmetrize_virial(virial)
335 ELSE
336 IF (calculate_stress_tensor) THEN
337 CALL cp_warn(__location__, "The calculation of the stress tensor "// &
338 "requires the calculation of the forces")
339 END IF
340 END IF
341 END IF
342 END IF
343
344 ! In case requested, compute the APT numerically
345 do_apt_fd = .false.
346 IF (force_env%in_use == use_qs_force) THEN
347 CALL section_vals_val_get(force_env%qs_env%input, "PROPERTIES%LINRES%DCDR%APT_FD", l_val=do_apt_fd)
348 IF (do_apt_fd) THEN
349 print_key => section_vals_get_subs_vals(force_env%qs_env%input, &
350 subsection_name="PROPERTIES%LINRES%DCDR%PRINT%APT")
351 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
352 CALL apt_fdiff(force_env)
353 END IF
354 END IF
355 END IF
356
357 !sample peak memory
358 CALL m_memory()
359
360 ! Some additional tasks..
361 IF (.NOT. my_skip) THEN
362 ! Flexible Partitioning
363 IF (ASSOCIATED(force_env%fp_env)) THEN
364 IF (force_env%fp_env%use_fp) THEN
365 CALL fp_eval(force_env%fp_env, subsys, cell)
366 END IF
367 END IF
368 ! Constraints ONLY of Fixed Atom type
369 CALL fix_atom_control(force_env)
370 ! All Restraints
371 CALL restraint_control(force_env)
372 ! Virtual Sites
373 CALL vsite_force_control(force_env)
374 ! External Potential
375 CALL add_external_potential(force_env)
376 ! Rescale forces if requested
377 CALL rescale_forces(force_env)
378 END IF
379
380 CALL force_env_get(force_env, potential_energy=e_pot)
381
382 ! Print energy always in the same format for all methods
383 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
384 extension=".Log")
385 IF (output_unit > 0) THEN
386 CALL section_vals_val_get(force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO%ENERGY_UNIT", &
387 c_val=unit_string)
388 fconv = cp_unit_from_cp2k(1.0_dp, trim(adjustl(unit_string)))
389 WRITE (unit=output_unit, fmt="(/,T2,A,T55,F26.15)") &
390 "ENERGY| Total FORCE_EVAL ( "//trim(adjustl(use_prog_name(force_env%in_use)))// &
391 " ) energy ["//trim(adjustl(unit_string))//"]", e_pot*fconv
392 IF (e_gap > -0.1_dp) THEN
393 WRITE (unit=output_unit, fmt="(/,T2,A,T55,F26.15)") &
394 "ENERGY| Total FORCE_EVAL ( "//trim(adjustl(use_prog_name(force_env%in_use)))// &
395 " ) gap ["//trim(adjustl(unit_string))//"]", e_gap*fconv
396 END IF
397 IF (e_entropy > -0.1_dp) THEN
398 WRITE (unit=output_unit, fmt="(/,T2,A,T55,F26.15)") &
399 "ENERGY| Total FORCE_EVAL ( "//trim(adjustl(use_prog_name(force_env%in_use)))// &
400 " ) free energy ["//trim(adjustl(unit_string))//"]", (e_pot - e_entropy)*fconv
401 END IF
402 END IF
403 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
404 "PRINT%PROGRAM_RUN_INFO")
405
406 ! terminate the run if the value of the potential is abnormal
407 IF (abnormal_value(e_pot)) THEN
408 cpabort("Potential energy is an abnormal value (NaN/Inf).")
409 END IF
410
411 ! Print forces, if requested
412 print_forces = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%FORCES", &
413 extension=".xyz")
414 IF ((print_forces > 0) .AND. calculate_forces) THEN
415 CALL force_env_get(force_env, subsys=subsys)
416 CALL cp_subsys_get(subsys, &
417 core_particles=core_particles, &
418 particles=particles, &
419 shell_particles=shell_particles)
420 ! Variable precision output of the forces
421 CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%NDIGITS", &
422 i_val=ndigits)
423 CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%FORCE_UNIT", &
424 c_val=unit_string)
425 IF (ASSOCIATED(core_particles) .OR. ASSOCIATED(shell_particles)) THEN
426 CALL write_forces(particles, print_forces, "Atomic", ndigits, unit_string, &
427 total_force, zero_force_core_shell_atom=.true.)
428 grand_total_force(1:3) = total_force(1:3)
429 IF (ASSOCIATED(core_particles)) THEN
430 CALL write_forces(core_particles, print_forces, "Core particle", ndigits, &
431 unit_string, total_force, zero_force_core_shell_atom=.false.)
432 grand_total_force(:) = grand_total_force(:) + total_force(:)
433 END IF
434 IF (ASSOCIATED(shell_particles)) THEN
435 CALL write_forces(shell_particles, print_forces, "Shell particle", ndigits, &
436 unit_string, total_force, zero_force_core_shell_atom=.false., &
437 grand_total_force=grand_total_force)
438 END IF
439 ELSE
440 CALL write_forces(particles, print_forces, "Atomic", ndigits, unit_string, total_force)
441 END IF
442 END IF
443 CALL cp_print_key_finished_output(print_forces, logger, force_env%force_env_section, "PRINT%FORCES")
444
445 ! Write stress tensor
446 IF (virial%pv_availability) THEN
447 ! If the virial is defined but we are not computing forces let's zero the
448 ! virial for consistency
449 IF (calculate_forces .AND. calculate_stress_tensor) THEN
450 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
451 extension=".stress_tensor")
452 IF (output_unit > 0) THEN
453 CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%COMPONENTS", &
454 l_val=print_components)
455 CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%STRESS_UNIT", &
456 c_val=unit_string)
457 IF (print_components) THEN
458 IF ((.NOT. virial%pv_numer) .AND. (force_env%in_use == use_qs_force)) THEN
459 CALL write_stress_tensor_components(virial, output_unit, cell, unit_string)
460 END IF
461 END IF
462 CALL write_stress_tensor(virial%pv_virial, output_unit, cell, unit_string, virial%pv_numer)
463 END IF
464 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
465 "PRINT%STRESS_TENSOR")
466 ELSE
467 CALL zero_virial(virial, reset=.false.)
468 END IF
469 ELSE
470 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
471 extension=".stress_tensor")
472 IF (output_unit > 0) THEN
473 CALL cp_warn(__location__, "To print the stress tensor switch on the "// &
474 "virial evaluation with the keyword: STRESS_TENSOR")
475 END IF
476 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
477 "PRINT%STRESS_TENSOR")
478 END IF
479
480 ! Atomic energy
481 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
482 extension=".Log")
483 IF (atprop_env%energy) THEN
484 CALL force_env%para_env%sum(atprop_env%atener)
485 CALL force_env_get(force_env, potential_energy=e_pot)
486 IF (output_unit > 0) THEN
487 IF (logger%iter_info%print_level >= low_print_level) THEN
488 CALL cp_subsys_get(subsys=subsys, particles=particles)
489 CALL write_atener(output_unit, particles, atprop_env%atener, "Mulliken Atomic Energies")
490 END IF
491 sum_energy = accurate_sum(atprop_env%atener(:))
492 checksum = abs(e_pot - sum_energy)
493 WRITE (unit=output_unit, fmt="(/,(T2,A,T56,F25.13))") &
494 "Potential energy (Atomic):", sum_energy, &
495 "Potential energy (Total) :", e_pot, &
496 "Difference :", checksum
497 cpassert((checksum < ateps*abs(e_pot)))
498 END IF
499 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
500 "PRINT%PROGRAM_RUN_INFO")
501 END IF
502
503 ! Print GRMM interface file
504 print_grrm = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%GRRM", &
505 file_position="REWIND", extension=".rrm")
506 IF (print_grrm > 0) THEN
507 CALL force_env_get(force_env, subsys=subsys)
508 CALL cp_subsys_get(subsys=subsys, particles=particles, &
509 molecule_kinds=molecule_kinds)
510 ! Count the number of fixed atoms
511 nfixed_atoms_total = 0
512 nkind = molecule_kinds%n_els
513 molecule_kind_set => molecule_kinds%els
514 DO ikind = 1, nkind
515 molecule_kind => molecule_kind_set(ikind)
516 CALL get_molecule_kind(molecule_kind, nfixd=nfixed_atoms)
517 nfixed_atoms_total = nfixed_atoms_total + nfixed_atoms
518 END DO
519 !
520 CALL write_grrm(print_grrm, force_env, particles%els, e_pot, fixed_atoms=nfixed_atoms_total)
521 END IF
522 CALL cp_print_key_finished_output(print_grrm, logger, force_env%force_env_section, "PRINT%GRRM")
523
524 print_scine = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%SCINE", &
525 file_position="REWIND", extension=".scine")
526 IF (print_scine > 0) THEN
527 CALL force_env_get(force_env, subsys=subsys)
528 CALL cp_subsys_get(subsys=subsys, particles=particles)
529 !
530 CALL write_scine(print_scine, force_env, particles%els, e_pot)
531 END IF
532 CALL cp_print_key_finished_output(print_scine, logger, force_env%force_env_section, "PRINT%SCINE")
533
534 END SUBROUTINE force_env_calc_energy_force
535
536! **************************************************************************************************
537!> \brief Rebuild k-point data for geometries whose atomic symmetry can change.
538!> Atomic k-point symmetry may change when atoms or cell vectors move.
539!> \param force_env ...
540!> \param fd_energy ...
541! **************************************************************************************************
542 SUBROUTINE force_env_refresh_kpoint_symmetry(force_env, fd_energy)
543
544 TYPE(force_env_type), POINTER :: force_env
545 LOGICAL, INTENT(IN) :: fd_energy
546
547 REAL(kind=dp), PARAMETER :: eps_cell = 1.0e-14_dp
548
549 CHARACTER(LEN=default_string_length) :: kp_scheme
550 INTEGER :: run_type_id
551 LOGICAL :: debug_full_kpoint_symmetry, debug_full_kpoint_symmetry_explicit, &
552 debug_inversion_only, do_kpoints, dynamic_symmetry, force_full_debug_symmetry, full_grid, &
553 input_full_grid, input_inversion_symmetry_only, inversion_symmetry_only, kpoint_symmetry, &
554 moving_geometry, non_lower_triangular_cell, use_full_grid, use_inversion_symmetry_only
555 TYPE(cell_type), POINTER :: cell
556 TYPE(cp_blacs_env_type), POINTER :: blacs_env
557 TYPE(dft_control_type), POINTER :: dft_control
558 TYPE(global_environment_type), POINTER :: globenv
559 TYPE(kpoint_type), POINTER :: kpoints
560 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
561 TYPE(mp_para_env_type), POINTER :: para_env
562 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
563 TYPE(qs_wf_history_type), POINTER :: wf_history
564 TYPE(section_vals_type), POINTER :: input, kpoint_section
565
566 IF (.NOT. ASSOCIATED(force_env)) RETURN
567 IF (force_env%in_use /= use_qs_force) RETURN
568
569 NULLIFY (globenv)
570 CALL force_env_get(force_env, globenv=globenv)
571 IF (.NOT. ASSOCIATED(globenv)) RETURN
572 run_type_id = globenv%run_type_id
573 moving_geometry = .false.
574 SELECT CASE (run_type_id)
575 CASE (cell_opt_run, driver_run, ehrenfest, geo_opt_run, mol_dyn_run)
576 moving_geometry = .true.
577 CASE DEFAULT
578 moving_geometry = .false.
579 END SELECT
580 IF (run_type_id /= debug_run .AND. .NOT. moving_geometry) RETURN
581
582 NULLIFY (blacs_env, cell, dft_control, input, kpoint_section, kpoints, mos, para_env, &
583 particle_set, wf_history)
584 CALL get_qs_env(qs_env=force_env%qs_env, &
585 blacs_env=blacs_env, &
586 cell=cell, &
587 dft_control=dft_control, &
588 do_kpoints=do_kpoints, &
589 input=input, &
590 kpoints=kpoints, &
591 mos=mos, &
592 para_env=para_env, &
593 particle_set=particle_set, &
594 wf_history=wf_history)
595 IF (.NOT. do_kpoints) RETURN
596
597 CALL get_kpoint_info(kpoints, kp_scheme=kp_scheme, symmetry=kpoint_symmetry, full_grid=full_grid, &
598 inversion_symmetry_only=inversion_symmetry_only)
599 IF (.NOT. kpoint_symmetry) RETURN
600 IF (trim(kp_scheme) /= "MONKHORST-PACK" .AND. trim(kp_scheme) /= "MACDONALD" .AND. &
601 trim(kp_scheme) /= "GENERAL") RETURN
602
603 input_full_grid = full_grid
604 input_inversion_symmetry_only = inversion_symmetry_only
605 debug_full_kpoint_symmetry = .false.
606 debug_full_kpoint_symmetry_explicit = .false.
607 IF (ASSOCIATED(input)) THEN
608 kpoint_section => section_vals_get_subs_vals(input, "DFT%KPOINTS")
609 CALL section_vals_val_get(kpoint_section, "FULL_GRID", l_val=input_full_grid)
610 CALL section_vals_val_get(kpoint_section, "INVERSION_SYMMETRY_ONLY", &
611 l_val=input_inversion_symmetry_only)
612 CALL section_vals_val_get(kpoint_section, "DEBUG_FULL_KPOINT_SYMMETRY", &
613 l_val=debug_full_kpoint_symmetry, &
614 explicit=debug_full_kpoint_symmetry_explicit)
615 END IF
616 ! Moving geometries and DEBUG finite differences must not reuse atomic symmetry from
617 ! another geometry. Rebuild the k-point symmetry from the current cell and positions.
618 ! An explicit DEBUG_FULL_KPOINT_SYMMETRY OFF keeps numerical finite-difference
619 ! energies and DFTB DEBUG checks on inversion/time-reversal reduction.
620 debug_inversion_only = run_type_id == debug_run .AND. .NOT. debug_full_kpoint_symmetry .AND. &
621 (fd_energy .OR. dft_control%qs_control%dftb)
622 force_full_debug_symmetry = run_type_id == debug_run .AND. debug_full_kpoint_symmetry_explicit .AND. &
623 debug_full_kpoint_symmetry
624 use_full_grid = input_full_grid
625 use_inversion_symmetry_only = (input_inversion_symmetry_only .OR. debug_inversion_only) .AND. &
626 (.NOT. use_full_grid)
627 ! Preserve restrictions selected during initial setup. The explicit DEBUG expert option
628 ! remains available for symmetry diagnostics on otherwise unsupported cell matrices.
629 IF (inversion_symmetry_only .AND. .NOT. force_full_debug_symmetry .AND. .NOT. use_full_grid) THEN
630 use_inversion_symmetry_only = .true.
631 END IF
632 non_lower_triangular_cell = (abs(cell%hmat(2, 1)) > eps_cell) .OR. &
633 (abs(cell%hmat(3, 1)) > eps_cell) .OR. &
634 (abs(cell%hmat(3, 2)) > eps_cell)
635 IF (non_lower_triangular_cell .AND. .NOT. force_full_debug_symmetry .AND. .NOT. use_full_grid) THEN
636 use_inversion_symmetry_only = .true.
637 END IF
638 dynamic_symmetry = kpoint_symmetry .AND. .NOT. use_full_grid .AND. &
639 .NOT. use_inversion_symmetry_only
640 IF (run_type_id == debug_run .AND. .NOT. fd_energy .AND. .NOT. dynamic_symmetry .AND. &
641 (full_grid .EQV. use_full_grid) .AND. &
642 (inversion_symmetry_only .EQV. use_inversion_symmetry_only)) THEN
643 CALL qs_basis_rotation(force_env%qs_env, kpoints)
644 RETURN
645 END IF
646 IF (moving_geometry .AND. .NOT. dynamic_symmetry) RETURN
647 IF (moving_geometry .AND. .NOT. kpoint_has_nontrivial_atomic_symmetry(kpoints)) RETURN
648 CALL set_kpoint_info(kpoints, full_grid=use_full_grid, &
649 inversion_symmetry_only=use_inversion_symmetry_only)
650
651 CALL kpoint_reset_initialization(kpoints)
652 CALL kpoint_initialize(kpoints, particle_set, cell)
653 CALL kpoint_env_initialize(kpoints, para_env, blacs_env, with_aux_fit=dft_control%do_admm)
654 CALL kpoint_initialize_mos(kpoints, mos)
655 CALL wfi_clear(wf_history)
656 CALL qs_basis_rotation(force_env%qs_env, kpoints)
657
658 END SUBROUTINE force_env_refresh_kpoint_symmetry
659
660! **************************************************************************************************
661!> \brief Return whether the current reduced mesh uses nontrivial atomic symmetry operations.
662!> \param kpoints ...
663!> \return has_symmetry
664! **************************************************************************************************
665 FUNCTION kpoint_has_nontrivial_atomic_symmetry(kpoints) RESULT(has_symmetry)
666
667 TYPE(kpoint_type), POINTER :: kpoints
668 LOGICAL :: has_symmetry
669
670 INTEGER :: iatom, ik, isym, natom
671 REAL(kind=dp), DIMENSION(3, 3) :: eye3
672 TYPE(kpoint_sym_type), POINTER :: kpsym
673
674 has_symmetry = .false.
675 IF (.NOT. ASSOCIATED(kpoints)) RETURN
676 IF (.NOT. ASSOCIATED(kpoints%kp_sym)) RETURN
677
678 eye3 = 0.0_dp
679 eye3(1, 1) = 1.0_dp
680 eye3(2, 2) = 1.0_dp
681 eye3(3, 3) = 1.0_dp
682
683 DO ik = 1, kpoints%nkp
684 kpsym => kpoints%kp_sym(ik)%kpoint_sym
685 IF (.NOT. ASSOCIATED(kpsym)) cycle
686 IF (.NOT. kpsym%apply_symmetry) cycle
687 IF (.NOT. ASSOCIATED(kpsym%rot)) cycle
688 IF (.NOT. ASSOCIATED(kpsym%f0)) cycle
689 IF (.NOT. ASSOCIATED(kpsym%fcell)) cycle
690
691 natom = SIZE(kpsym%f0, 1)
692 DO isym = 1, SIZE(kpsym%rot, 3)
693 IF (maxval(abs(kpsym%rot(1:3, 1:3, isym) - eye3(1:3, 1:3))) > 1.e-12_dp .OR. &
694 any(kpsym%fcell(1:3, 1:natom, isym) /= 0)) THEN
695 has_symmetry = .true.
696 RETURN
697 END IF
698 DO iatom = 1, natom
699 IF (kpsym%f0(iatom, isym) /= iatom) THEN
700 has_symmetry = .true.
701 RETURN
702 END IF
703 END DO
704 END DO
705 END DO
706
707 END FUNCTION kpoint_has_nontrivial_atomic_symmetry
708
709! **************************************************************************************************
710!> \brief Evaluates the stress tensor and pressure numerically
711!> \param force_env ...
712!> \param dx ...
713!> \par History
714!> 10.2005 created [JCS]
715!> 05.2009 Teodoro Laino [tlaino] - rewriting for general force_env
716!>
717!> \author JCS
718! **************************************************************************************************
719 SUBROUTINE force_env_calc_num_pressure(force_env, dx)
720
721 TYPE(force_env_type), POINTER :: force_env
722 REAL(kind=dp), INTENT(IN), OPTIONAL :: dx
723
724 REAL(kind=dp), PARAMETER :: default_dx = 0.001_dp
725
726 CHARACTER(LEN=default_string_length) :: unit_string
727 INTEGER :: i, ip, iq, j, k, method_id, natom, &
728 ncore, nshell, output_unit, symmetry_id
729 LOGICAL :: use_sym_strain_2d
730 REAL(kind=dp) :: dx_w, eps_w
731 REAL(kind=dp), DIMENSION(2) :: numer_energy
732 REAL(kind=dp), DIMENSION(3) :: s
733 REAL(kind=dp), DIMENSION(3, 3) :: hmat_deformed, numer_pv_2d, &
734 numer_stress, strain
735 REAL(kind=dp), DIMENSION(:, :), POINTER :: ref_pos_atom, ref_pos_core, ref_pos_shell
736 TYPE(cell_type), POINTER :: cell, cell_local
737 TYPE(cp_logger_type), POINTER :: logger
738 TYPE(cp_subsys_type), POINTER :: subsys
739 TYPE(dft_control_type), POINTER :: dft_control
740 TYPE(global_environment_type), POINTER :: globenv
741 TYPE(particle_list_type), POINTER :: core_particles, particles, &
742 shell_particles
743 TYPE(virial_type), POINTER :: virial
744
745 NULLIFY (cell_local)
746 NULLIFY (dft_control)
747 NULLIFY (core_particles)
748 NULLIFY (particles)
749 NULLIFY (shell_particles)
750 NULLIFY (ref_pos_atom)
751 NULLIFY (ref_pos_core)
752 NULLIFY (ref_pos_shell)
753 natom = 0
754 method_id = 0
755 ncore = 0
756 nshell = 0
757 numer_pv_2d = 0.0_dp
758 numer_stress = 0.0_dp
759 use_sym_strain_2d = .false.
760
761 logger => cp_get_default_logger()
762
763 dx_w = default_dx
764 IF (PRESENT(dx)) dx_w = dx
765 CALL force_env_get(force_env, subsys=subsys, globenv=globenv, in_use=method_id)
766 CALL cp_subsys_get(subsys, &
767 core_particles=core_particles, &
768 particles=particles, &
769 shell_particles=shell_particles, &
770 virial=virial)
771 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
772 extension=".stress_tensor")
773 IF (output_unit > 0) THEN
774 WRITE (output_unit, "(/A,A/)") " **************************** ", &
775 "NUMERICAL STRESS ********************************"
776 END IF
777
778 ! Save all original particle positions
779 natom = particles%n_els
780 ALLOCATE (ref_pos_atom(natom, 3))
781 DO i = 1, natom
782 ref_pos_atom(i, :) = particles%els(i)%r
783 END DO
784 IF (ASSOCIATED(core_particles)) THEN
785 ncore = core_particles%n_els
786 ALLOCATE (ref_pos_core(ncore, 3))
787 DO i = 1, ncore
788 ref_pos_core(i, :) = core_particles%els(i)%r
789 END DO
790 END IF
791 IF (ASSOCIATED(shell_particles)) THEN
792 nshell = shell_particles%n_els
793 ALLOCATE (ref_pos_shell(nshell, 3))
794 DO i = 1, nshell
795 ref_pos_shell(i, :) = shell_particles%els(i)%r
796 END DO
797 END IF
798 CALL force_env_get(force_env, cell=cell)
799 ! Save cell symmetry (distorted cell has no symmetry)
800 symmetry_id = cell%symmetry_id
801 cell%symmetry_id = cell_sym_triclinic
802 !
803 CALL cell_create(cell_local)
804 CALL cell_clone(cell, cell_local)
805 IF (count(cell_local%perd /= 0) == 2 .AND. method_id == use_qs_force) THEN
806 CALL get_qs_env(qs_env=force_env%qs_env, dft_control=dft_control)
807 SELECT CASE (dft_control%qs_control%method_id)
808 CASE (do_method_gapw, do_method_gapw_xc, do_method_gpw, &
809 do_method_lrigpw, do_method_ofgpw, do_method_rigpw)
810 use_sym_strain_2d = .true.
811 END SELECT
812 END IF
813 ! First change box
814 DO ip = 1, 3
815 DO iq = 1, 3
816 IF (use_sym_strain_2d) THEN
817 IF (cell_local%perd(ip) == 0 .OR. cell_local%perd(iq) == 0) cycle
818 IF (iq < ip) cycle
819 END IF
820 IF (virial%pv_diagonal .AND. (ip /= iq)) cycle
821 DO k = 1, 2
822 hmat_deformed = cell_local%hmat
823 IF (use_sym_strain_2d) THEN
824 eps_w = -(-1.0_dp)**k*dx_w
825 strain = 0.0_dp
826 DO i = 1, 3
827 strain(i, i) = 1.0_dp
828 END DO
829 IF (ip == iq) THEN
830 strain(ip, ip) = strain(ip, ip) + eps_w
831 ELSE
832 strain(ip, iq) = strain(ip, iq) + 0.5_dp*eps_w
833 strain(iq, ip) = strain(iq, ip) + 0.5_dp*eps_w
834 END IF
835 hmat_deformed = matmul(strain, cell_local%hmat)
836 ELSE
837 hmat_deformed(ip, iq) = hmat_deformed(ip, iq) - (-1.0_dp)**k*dx_w
838 END IF
839 cell%hmat = hmat_deformed
840 CALL init_cell(cell)
841 ! Scale positions
842 DO i = 1, natom
843 CALL real_to_scaled(s, ref_pos_atom(i, 1:3), cell_local)
844 CALL scaled_to_real(particles%els(i)%r, s, cell)
845 END DO
846 DO i = 1, ncore
847 CALL real_to_scaled(s, ref_pos_core(i, 1:3), cell_local)
848 CALL scaled_to_real(core_particles%els(i)%r, s, cell)
849 END DO
850 DO i = 1, nshell
851 CALL real_to_scaled(s, ref_pos_shell(i, 1:3), cell_local)
852 CALL scaled_to_real(shell_particles%els(i)%r, s, cell)
853 END DO
854 ! Compute energies
855 CALL force_env_calc_energy_force(force_env, &
856 calc_force=.false., &
857 consistent_energies=.true., &
858 calc_stress_tensor=.false.)
859 CALL force_env_get(force_env, potential_energy=numer_energy(k))
860 ! Reset cell
861 cell%hmat = cell_local%hmat
862 END DO
863 CALL init_cell(cell)
864 IF (use_sym_strain_2d) THEN
865 numer_pv_2d(ip, iq) = -0.5_dp*(numer_energy(1) - numer_energy(2))/dx_w
866 numer_pv_2d(iq, ip) = numer_pv_2d(ip, iq)
867 IF (output_unit > 0) THEN
868 IF (globenv%run_type_id == debug_run) THEN
869 WRITE (unit=output_unit, fmt="(/,T2,A,T19,A,F7.4,A,T44,A,F7.4,A,T69,A)") &
870 "DEBUG|", "E(e"//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
871 "E(e"//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
872 "pv(numerical)"
873 WRITE (unit=output_unit, fmt="(T2,A,2(1X,F24.8),1X,F22.8)") &
874 "DEBUG|", numer_energy(1:2), numer_pv_2d(ip, iq)
875 ELSE
876 WRITE (unit=output_unit, fmt="(/,T7,A,F7.4,A,T27,A,F7.4,A,T49,A)") &
877 "E(e"//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
878 "E(e"//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
879 "pv(numerical)"
880 WRITE (unit=output_unit, fmt="(3(1X,F19.8))") &
881 numer_energy(1:2), numer_pv_2d(ip, iq)
882 END IF
883 END IF
884 ELSE
885 numer_stress(ip, iq) = 0.5_dp*(numer_energy(1) - numer_energy(2))/dx_w
886 IF (output_unit > 0) THEN
887 IF (globenv%run_type_id == debug_run) THEN
888 WRITE (unit=output_unit, fmt="(/,T2,A,T19,A,F7.4,A,T44,A,F7.4,A,T69,A)") &
889 "DEBUG|", "E("//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
890 "E("//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
891 "f(numerical)"
892 WRITE (unit=output_unit, fmt="(T2,A,2(1X,F24.8),1X,F22.8)") &
893 "DEBUG|", numer_energy(1:2), numer_stress(ip, iq)
894 ELSE
895 WRITE (unit=output_unit, fmt="(/,T7,A,F7.4,A,T27,A,F7.4,A,T49,A)") &
896 "E("//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
897 "E("//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
898 "f(numerical)"
899 WRITE (unit=output_unit, fmt="(3(1X,F19.8))") &
900 numer_energy(1:2), numer_stress(ip, iq)
901 END IF
902 END IF
903 END IF
904 END DO
905 END DO
906
907 ! Reset positions and rebuild original environment
908 cell%symmetry_id = symmetry_id
909 CALL init_cell(cell)
910 DO i = 1, natom
911 particles%els(i)%r = ref_pos_atom(i, :)
912 END DO
913 DO i = 1, ncore
914 core_particles%els(i)%r = ref_pos_core(i, :)
915 END DO
916 DO i = 1, nshell
917 shell_particles%els(i)%r = ref_pos_shell(i, :)
918 END DO
919 CALL force_env_calc_energy_force(force_env, &
920 calc_force=.false., &
921 consistent_energies=.true., &
922 calc_stress_tensor=.false.)
923
924 ! Computing pv_test
925 virial%pv_virial = 0.0_dp
926 IF (use_sym_strain_2d) THEN
927 virial%pv_virial = numer_pv_2d
928 ELSE
929 DO i = 1, 3
930 DO j = 1, 3
931 DO k = 1, 3
932 virial%pv_virial(i, j) = virial%pv_virial(i, j) - &
933 0.5_dp*(numer_stress(i, k)*cell_local%hmat(j, k) + &
934 numer_stress(j, k)*cell_local%hmat(i, k))
935 END DO
936 END DO
937 END DO
938 END IF
939 IF (output_unit > 0) THEN
940 IF (globenv%run_type_id == debug_run) THEN
941 CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%FORCE_UNIT", &
942 c_val=unit_string)
943 CALL write_stress_tensor(virial%pv_virial, output_unit, cell, unit_string, virial%pv_numer)
944 END IF
945 WRITE (output_unit, "(/,A,/)") " **************************** "// &
946 "NUMERICAL STRESS END *****************************"
947 END IF
948
949 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
950 "PRINT%STRESS_TENSOR")
951
952 ! Release storage
953 IF (ASSOCIATED(ref_pos_atom)) THEN
954 DEALLOCATE (ref_pos_atom)
955 END IF
956 IF (ASSOCIATED(ref_pos_core)) THEN
957 DEALLOCATE (ref_pos_core)
958 END IF
959 IF (ASSOCIATED(ref_pos_shell)) THEN
960 DEALLOCATE (ref_pos_shell)
961 END IF
962 IF (ASSOCIATED(cell_local)) CALL cell_release(cell_local)
963
964 END SUBROUTINE force_env_calc_num_pressure
965
966! **************************************************************************************************
967!> \brief creates and initializes a force environment
968!> \param force_env the force env to create
969!> \param root_section ...
970!> \param para_env ...
971!> \param globenv ...
972!> \param fist_env , qs_env: exactly one of these should be
973!> associated, the one that is active
974!> \param qs_env ...
975!> \param meta_env ...
976!> \param sub_force_env ...
977!> \param qmmm_env ...
978!> \param qmmmx_env ...
979!> \param eip_env ...
980!> \param pwdft_env ...
981!> \param force_env_section ...
982!> \param mixed_env ...
983!> \param embed_env ...
984!> \param nnp_env ...
985!> \param ipi_env ...
986!> \par History
987!> 04.2003 created [fawzi]
988!> \author fawzi
989! **************************************************************************************************
990 SUBROUTINE force_env_create(force_env, root_section, para_env, globenv, fist_env, &
991 qs_env, meta_env, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, force_env_section, &
992 mixed_env, embed_env, nnp_env, ipi_env)
993
994 TYPE(force_env_type), POINTER :: force_env
995 TYPE(section_vals_type), POINTER :: root_section
996 TYPE(mp_para_env_type), POINTER :: para_env
997 TYPE(global_environment_type), POINTER :: globenv
998 TYPE(fist_environment_type), OPTIONAL, POINTER :: fist_env
999 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
1000 TYPE(meta_env_type), OPTIONAL, POINTER :: meta_env
1001 TYPE(force_env_p_type), DIMENSION(:), OPTIONAL, &
1002 POINTER :: sub_force_env
1003 TYPE(qmmm_env_type), OPTIONAL, POINTER :: qmmm_env
1004 TYPE(qmmmx_env_type), OPTIONAL, POINTER :: qmmmx_env
1005 TYPE(eip_environment_type), OPTIONAL, POINTER :: eip_env
1006 TYPE(pwdft_environment_type), OPTIONAL, POINTER :: pwdft_env
1007 TYPE(section_vals_type), POINTER :: force_env_section
1008 TYPE(mixed_environment_type), OPTIONAL, POINTER :: mixed_env
1009 TYPE(embed_env_type), OPTIONAL, POINTER :: embed_env
1010 TYPE(nnp_type), OPTIONAL, POINTER :: nnp_env
1011 TYPE(ipi_environment_type), OPTIONAL, POINTER :: ipi_env
1012
1013 ALLOCATE (force_env)
1014 NULLIFY (force_env%fist_env, force_env%qs_env, &
1015 force_env%para_env, force_env%globenv, &
1016 force_env%meta_env, force_env%sub_force_env, &
1017 force_env%qmmm_env, force_env%qmmmx_env, force_env%fp_env, &
1018 force_env%force_env_section, force_env%eip_env, force_env%mixed_env, &
1019 force_env%embed_env, force_env%pwdft_env, force_env%nnp_env, &
1020 force_env%root_section)
1021 last_force_env_id = last_force_env_id + 1
1022 force_env%ref_count = 1
1023 force_env%in_use = 0
1024 force_env%additional_potential = 0.0_dp
1025
1026 force_env%globenv => globenv
1027 CALL globenv_retain(force_env%globenv)
1028
1029 force_env%root_section => root_section
1030 CALL section_vals_retain(root_section)
1031
1032 force_env%para_env => para_env
1033 CALL force_env%para_env%retain()
1034
1035 CALL section_vals_retain(force_env_section)
1036 force_env%force_env_section => force_env_section
1037
1038 IF (PRESENT(fist_env)) THEN
1039 cpassert(ASSOCIATED(fist_env))
1040 cpassert(force_env%in_use == 0)
1041 force_env%in_use = use_fist_force
1042 force_env%fist_env => fist_env
1043 END IF
1044 IF (PRESENT(eip_env)) THEN
1045 cpassert(ASSOCIATED(eip_env))
1046 cpassert(force_env%in_use == 0)
1047 force_env%in_use = use_eip_force
1048 force_env%eip_env => eip_env
1049 END IF
1050 IF (PRESENT(pwdft_env)) THEN
1051 cpassert(ASSOCIATED(pwdft_env))
1052 cpassert(force_env%in_use == 0)
1053 force_env%in_use = use_pwdft_force
1054 force_env%pwdft_env => pwdft_env
1055 END IF
1056 IF (PRESENT(qs_env)) THEN
1057 cpassert(ASSOCIATED(qs_env))
1058 cpassert(force_env%in_use == 0)
1059 force_env%in_use = use_qs_force
1060 force_env%qs_env => qs_env
1061 END IF
1062 IF (PRESENT(qmmm_env)) THEN
1063 cpassert(ASSOCIATED(qmmm_env))
1064 cpassert(force_env%in_use == 0)
1065 force_env%in_use = use_qmmm
1066 force_env%qmmm_env => qmmm_env
1067 END IF
1068 IF (PRESENT(qmmmx_env)) THEN
1069 cpassert(ASSOCIATED(qmmmx_env))
1070 cpassert(force_env%in_use == 0)
1071 force_env%in_use = use_qmmmx
1072 force_env%qmmmx_env => qmmmx_env
1073 END IF
1074 IF (PRESENT(mixed_env)) THEN
1075 cpassert(ASSOCIATED(mixed_env))
1076 cpassert(force_env%in_use == 0)
1077 force_env%in_use = use_mixed_force
1078 force_env%mixed_env => mixed_env
1079 END IF
1080 IF (PRESENT(embed_env)) THEN
1081 cpassert(ASSOCIATED(embed_env))
1082 cpassert(force_env%in_use == 0)
1083 force_env%in_use = use_embed
1084 force_env%embed_env => embed_env
1085 END IF
1086 IF (PRESENT(nnp_env)) THEN
1087 cpassert(ASSOCIATED(nnp_env))
1088 cpassert(force_env%in_use == 0)
1089 force_env%in_use = use_nnp_force
1090 force_env%nnp_env => nnp_env
1091 END IF
1092 IF (PRESENT(ipi_env)) THEN
1093 cpassert(ASSOCIATED(ipi_env))
1094 cpassert(force_env%in_use == 0)
1095 force_env%in_use = use_ipi
1096 force_env%ipi_env => ipi_env
1097 END IF
1098 cpassert(force_env%in_use /= 0)
1099
1100 IF (PRESENT(sub_force_env)) THEN
1101 force_env%sub_force_env => sub_force_env
1102 END IF
1103
1104 IF (PRESENT(meta_env)) THEN
1105 force_env%meta_env => meta_env
1106 ELSE
1107 NULLIFY (force_env%meta_env)
1108 END IF
1109
1110 END SUBROUTINE force_env_create
1111
1112! **************************************************************************************************
1113!> \brief ****f* force_env_methods/mixed_energy_forces [1.0]
1114!>
1115!> Computes energy and forces for a mixed force_env type
1116!> \param force_env the force_env that holds the mixed_env type
1117!> \param calculate_forces decides if forces should be calculated
1118!> \par History
1119!> 11.06 created [fschiff]
1120!> 04.07 generalization to an illimited number of force_eval [tlaino]
1121!> 04.07 further generalization to force_eval with different geometrical
1122!> structures [tlaino]
1123!> 04.08 reorganizing the genmix structure (collecting common code)
1124!> 01.16 added CDFT [Nico Holmberg]
1125!> 08.17 added DFT embedding [Vladimir Rybkin]
1126!> \author Florian Schiffmann
1127! **************************************************************************************************
1128 SUBROUTINE mixed_energy_forces(force_env, calculate_forces)
1129
1130 TYPE(force_env_type), POINTER :: force_env
1131 LOGICAL, INTENT(IN) :: calculate_forces
1132
1133 CHARACTER(LEN=default_path_length) :: coupling_function
1134 CHARACTER(LEN=default_string_length) :: def_error, description, this_error
1135 INTEGER :: iforce_eval, iparticle, istate(2), &
1136 jparticle, mixing_type, my_group, &
1137 natom, nforce_eval, source, unit_nr
1138 INTEGER, DIMENSION(:), POINTER :: glob_natoms, itmplist, map_index
1139 LOGICAL :: dip_exists
1140 REAL(kind=dp) :: coupling_parameter, dedf, der_1, der_2, &
1141 dx, energy, err, lambda, lerr, &
1142 restraint_strength, restraint_target, &
1143 sd
1144 REAL(kind=dp), DIMENSION(3) :: dip_mix
1145 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1146 TYPE(cell_type), POINTER :: cell_mix
1147 TYPE(cp_logger_type), POINTER :: logger, my_logger
1148 TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1149 TYPE(cp_result_type), POINTER :: loc_results, results_mix
1150 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1151 TYPE(cp_subsys_type), POINTER :: subsys_mix
1152 TYPE(mixed_energy_type), POINTER :: mixed_energy
1153 TYPE(mixed_force_type), DIMENSION(:), POINTER :: global_forces
1154 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1155 TYPE(particle_list_type), POINTER :: particles_mix
1156 TYPE(section_vals_type), POINTER :: force_env_section, gen_section, &
1157 mapping_section, mixed_section, &
1158 root_section
1159 TYPE(virial_p_type), DIMENSION(:), POINTER :: virials
1160 TYPE(virial_type), POINTER :: loc_virial, virial_mix
1161
1162 logger => cp_get_default_logger()
1163 cpassert(ASSOCIATED(force_env))
1164 ! Get infos about the mixed subsys
1165 CALL force_env_get(force_env=force_env, &
1166 subsys=subsys_mix, &
1167 force_env_section=force_env_section, &
1168 root_section=root_section, &
1169 cell=cell_mix)
1170 CALL cp_subsys_get(subsys=subsys_mix, &
1171 particles=particles_mix, &
1172 virial=virial_mix, &
1173 results=results_mix)
1174 NULLIFY (map_index, glob_natoms, global_forces, itmplist)
1175
1176 nforce_eval = SIZE(force_env%sub_force_env)
1177 mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
1178 mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
1179 ! Global Info
1180 ALLOCATE (subsystems(nforce_eval))
1181 ALLOCATE (particles(nforce_eval))
1182 ! Local Info to sync
1183 ALLOCATE (global_forces(nforce_eval))
1184 ALLOCATE (energies(nforce_eval))
1185 ALLOCATE (glob_natoms(nforce_eval))
1186 ALLOCATE (virials(nforce_eval))
1187 ALLOCATE (results(nforce_eval))
1188 energies = 0.0_dp
1189 glob_natoms = 0
1190 ! Check if mixed CDFT calculation is requested and initialize
1191 CALL mixed_cdft_init(force_env, calculate_forces)
1192
1193 !
1194 IF (.NOT. force_env%mixed_env%do_mixed_cdft) THEN
1195 DO iforce_eval = 1, nforce_eval
1196 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1197 NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
1198 ALLOCATE (virials(iforce_eval)%virial)
1199 CALL cp_result_create(results(iforce_eval)%results)
1200 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1201 ! From this point on the error is the sub_error
1202 my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
1203 my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
1204 ! Copy iterations info (they are updated only in the main mixed_env)
1205 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1206 CALL cp_add_default_logger(my_logger)
1207
1208 ! Get all available subsys
1209 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1210 subsys=subsystems(iforce_eval)%subsys)
1211
1212 ! all force_env share the same cell
1213 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
1214
1215 ! Get available particles
1216 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1217 particles=particles(iforce_eval)%list)
1218
1219 ! Get Mapping index array
1220 natom = SIZE(particles(iforce_eval)%list%els)
1221
1222 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1223 map_index)
1224
1225 ! Mapping particles from iforce_eval environment to the mixed env
1226 DO iparticle = 1, natom
1227 jparticle = map_index(iparticle)
1228 particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
1229 END DO
1230
1231 ! Calculate energy and forces for each sub_force_env
1232 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1233 calc_force=calculate_forces, &
1234 skip_external_control=.true.)
1235
1236 ! Only the rank 0 process collect info for each computation
1237 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1238 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1239 potential_energy=energy)
1240 CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1241 virial=loc_virial, results=loc_results)
1242 energies(iforce_eval) = energy
1243 glob_natoms(iforce_eval) = natom
1244 virials(iforce_eval)%virial = loc_virial
1245 CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1246 END IF
1247 ! Deallocate map_index array
1248 IF (ASSOCIATED(map_index)) THEN
1249 DEALLOCATE (map_index)
1250 END IF
1251 CALL cp_rm_default_logger()
1252 END DO
1253 ELSE
1254 CALL mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
1255 glob_natoms, virials, results)
1256 END IF
1257 ! Handling Parallel execution
1258 CALL force_env%para_env%sync()
1259 ! Post CDFT operations
1260 CALL mixed_cdft_post_energy_forces(force_env)
1261 ! Let's transfer energy, natom, forces, virials
1262 CALL force_env%para_env%sum(energies)
1263 CALL force_env%para_env%sum(glob_natoms)
1264 ! Transfer forces
1265 DO iforce_eval = 1, nforce_eval
1266 ALLOCATE (global_forces(iforce_eval)%forces(3, glob_natoms(iforce_eval)))
1267 global_forces(iforce_eval)%forces = 0.0_dp
1268 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN
1269 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1270 ! Forces
1271 DO iparticle = 1, glob_natoms(iforce_eval)
1272 global_forces(iforce_eval)%forces(:, iparticle) = &
1273 particles(iforce_eval)%list%els(iparticle)%f
1274 END DO
1275 END IF
1276 END IF
1277 CALL force_env%para_env%sum(global_forces(iforce_eval)%forces)
1278 !Transfer only the relevant part of the virial..
1279 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_total)
1280 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_kinetic)
1281 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_virial)
1282 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_xc)
1283 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_fock_4c)
1284 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_constraint)
1285 !Transfer results
1286 source = 0
1287 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN
1288 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1289 source = force_env%para_env%mepos
1290 END IF
1291 END IF
1292 CALL force_env%para_env%sum(source)
1293 CALL cp_results_mp_bcast(results(iforce_eval)%results, source, force_env%para_env)
1294 END DO
1295
1296 IF (ALLOCATED(force_env%mixed_env%energies)) THEN
1297 DEALLOCATE (force_env%mixed_env%energies)
1298 END IF
1299 ALLOCATE (force_env%mixed_env%energies(nforce_eval))
1300 force_env%mixed_env%energies(:) = energies(:)
1301 ! Start combining the different sub_force_env
1302 CALL get_mixed_env(mixed_env=force_env%mixed_env, &
1303 mixed_energy=mixed_energy)
1304
1305 !NB: do this for all MIXING_TYPE values, since some need it (e.g. linear mixing
1306 !NB if the first system has fewer atoms than the second)
1307 DO iparticle = 1, SIZE(particles_mix%els)
1308 particles_mix%els(iparticle)%f(:) = 0.0_dp
1309 END DO
1310
1311 CALL section_vals_val_get(mixed_section, "MIXING_TYPE", i_val=mixing_type)
1312 SELECT CASE (mixing_type)
1313 CASE (mix_linear_combination)
1314 ! Support offered only 2 force_eval
1315 cpassert(nforce_eval == 2)
1316 CALL section_vals_val_get(mixed_section, "LINEAR%LAMBDA", r_val=lambda)
1317 mixed_energy%pot = lambda*energies(1) + (1.0_dp - lambda)*energies(2)
1318 ! General Mapping of forces...
1319 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1320 lambda, 1, nforce_eval, map_index, mapping_section, .true.)
1321 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1322 (1.0_dp - lambda), 2, nforce_eval, map_index, mapping_section, .false.)
1323 CASE (mix_minimum)
1324 ! Support offered only 2 force_eval
1325 cpassert(nforce_eval == 2)
1326 IF (energies(1) < energies(2)) THEN
1327 mixed_energy%pot = energies(1)
1328 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1329 1.0_dp, 1, nforce_eval, map_index, mapping_section, .true.)
1330 ELSE
1331 mixed_energy%pot = energies(2)
1332 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1333 1.0_dp, 2, nforce_eval, map_index, mapping_section, .true.)
1334 END IF
1335 CASE (mix_coupled)
1336 ! Support offered only 2 force_eval
1337 cpassert(nforce_eval == 2)
1338 CALL section_vals_val_get(mixed_section, "COUPLING%COUPLING_PARAMETER", &
1339 r_val=coupling_parameter)
1340 sd = sqrt((energies(1) - energies(2))**2 + 4.0_dp*coupling_parameter**2)
1341 der_1 = (1.0_dp - (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1342 der_2 = (1.0_dp + (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1343 mixed_energy%pot = (energies(1) + energies(2) - sd)/2.0_dp
1344 ! General Mapping of forces...
1345 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1346 der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1347 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1348 der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1349 CASE (mix_restrained)
1350 ! Support offered only 2 force_eval
1351 cpassert(nforce_eval == 2)
1352 CALL section_vals_val_get(mixed_section, "RESTRAINT%RESTRAINT_TARGET", &
1353 r_val=restraint_target)
1354 CALL section_vals_val_get(mixed_section, "RESTRAINT%RESTRAINT_STRENGTH", &
1355 r_val=restraint_strength)
1356 mixed_energy%pot = energies(1) + restraint_strength*(energies(1) - energies(2) - restraint_target)**2
1357 der_2 = -2.0_dp*restraint_strength*(energies(1) - energies(2) - restraint_target)
1358 der_1 = 1.0_dp - der_2
1359 ! General Mapping of forces...
1360 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1361 der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1362 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1363 der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1364 CASE (mix_generic)
1365 ! Support any number of force_eval sections
1366 gen_section => section_vals_get_subs_vals(mixed_section, "GENERIC")
1367 CALL parse_function_section(gen_section, "MIXING_FUNCTION", coupling_function, &
1368 force_env%mixed_env%par, force_env%mixed_env%val, &
1369 var_values=energies, dx=dx, error_limit=lerr)
1370 CALL initf(1)
1371 CALL parsef(1, trim(coupling_function), force_env%mixed_env%par)
1372 ! Now the hardest part.. map energy with corresponding force_eval
1373 mixed_energy%pot = evalf(1, force_env%mixed_env%val)
1374 cpassert(evalerrtype <= 0)
1375 CALL zero_virial(virial_mix, reset=.false.)
1376 CALL cp_results_erase(results_mix)
1377 DO iforce_eval = 1, nforce_eval
1378 dedf = evalfd(1, iforce_eval, force_env%mixed_env%val, dx, err)
1379 IF (abs(err) > lerr) THEN
1380 WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
1381 WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
1382 CALL compress(this_error, .true.)
1383 CALL compress(def_error, .true.)
1384 CALL cp_warn(__location__, &
1385 'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
1386 ' Error '//trim(this_error)//' in computing numerical derivatives larger than'// &
1387 trim(def_error)//' .')
1388 END IF
1389 ! General Mapping of forces...
1390 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1391 dedf, iforce_eval, nforce_eval, map_index, mapping_section, .false.)
1392 force_env%mixed_env%val(iforce_eval) = energies(iforce_eval)
1393 END DO
1394 ! Let's store the needed information..
1395 force_env%mixed_env%dx = dx
1396 force_env%mixed_env%lerr = lerr
1397 force_env%mixed_env%coupling_function = trim(coupling_function)
1398 CALL finalizef()
1399 CASE (mix_cdft)
1400 ! Supports any number of force_evals for calculation of CDFT properties, but forces only from two
1401 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LAMBDA", r_val=lambda)
1402 ! Get the states which determine the forces
1403 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%FORCE_STATES", i_vals=itmplist)
1404 IF (SIZE(itmplist) /= 2) THEN
1405 CALL cp_abort(__location__, &
1406 "Keyword FORCE_STATES takes exactly two input values.")
1407 END IF
1408 IF (any(itmplist < 0)) THEN
1409 cpabort("Invalid force_eval index.")
1410 END IF
1411 istate = itmplist
1412 IF (istate(1) > nforce_eval .OR. istate(2) > nforce_eval) THEN
1413 cpabort("Invalid force_eval index.")
1414 END IF
1415 mixed_energy%pot = lambda*energies(istate(1)) + (1.0_dp - lambda)*energies(istate(2))
1416 ! General Mapping of forces...
1417 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1418 lambda, istate(1), nforce_eval, map_index, mapping_section, .true.)
1419 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1420 (1.0_dp - lambda), istate(2), nforce_eval, map_index, mapping_section, .false.)
1421 CASE DEFAULT
1422 cpabort("Unknown mixing type for mixed_energy_forces")
1423 END SELECT
1424 !Simply deallocate and loose the pointer references..
1425 DO iforce_eval = 1, nforce_eval
1426 DEALLOCATE (global_forces(iforce_eval)%forces)
1427 IF (ASSOCIATED(virials(iforce_eval)%virial)) DEALLOCATE (virials(iforce_eval)%virial)
1428 CALL cp_result_release(results(iforce_eval)%results)
1429 END DO
1430 DEALLOCATE (global_forces)
1431 DEALLOCATE (subsystems)
1432 DEALLOCATE (particles)
1433 DEALLOCATE (energies)
1434 DEALLOCATE (glob_natoms)
1435 DEALLOCATE (virials)
1436 DEALLOCATE (results)
1437 ! Print Section
1438 unit_nr = cp_print_key_unit_nr(logger, mixed_section, "PRINT%DIPOLE", &
1439 extension=".data", middle_name="MIXED_DIPOLE", log_filename=.false.)
1440 IF (unit_nr > 0) THEN
1441 description = '[DIPOLE]'
1442 dip_exists = test_for_result(results=results_mix, description=description)
1443 IF (dip_exists) THEN
1444 CALL get_results(results=results_mix, description=description, values=dip_mix)
1445 WRITE (unit_nr, '(/,1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE ( A.U.)|", dip_mix
1446 WRITE (unit_nr, '( 1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE (Debye)|", dip_mix*debye
1447 ELSE
1448 WRITE (unit_nr, *) "NO FORCE_EVAL section calculated the dipole"
1449 END IF
1450 END IF
1451 CALL cp_print_key_finished_output(unit_nr, logger, mixed_section, "PRINT%DIPOLE")
1452 END SUBROUTINE mixed_energy_forces
1453
1454! **************************************************************************************************
1455!> \brief Driver routine for mixed CDFT energy and force calculations
1456!> \param force_env the force_env that holds the mixed_env
1457!> \param calculate_forces if forces should be calculated
1458!> \param particles system particles
1459!> \param energies the energies of the CDFT states
1460!> \param glob_natoms the total number of particles
1461!> \param virials the virials stored in subsys
1462!> \param results results stored in subsys
1463!> \par History
1464!> 01.17 created [Nico Holmberg]
1465!> \author Nico Holmberg
1466! **************************************************************************************************
1467 SUBROUTINE mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
1468 glob_natoms, virials, results)
1469 TYPE(force_env_type), POINTER :: force_env
1470 LOGICAL, INTENT(IN) :: calculate_forces
1471 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1472 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1473 INTEGER, DIMENSION(:), POINTER :: glob_natoms
1474 TYPE(virial_p_type), DIMENSION(:), POINTER :: virials
1475 TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1476
1477 INTEGER :: iforce_eval, iparticle, jparticle, &
1478 my_group, natom, nforce_eval
1479 INTEGER, DIMENSION(:), POINTER :: map_index
1480 REAL(kind=dp) :: energy
1481 TYPE(cell_type), POINTER :: cell_mix
1482 TYPE(cp_logger_type), POINTER :: logger, my_logger
1483 TYPE(cp_result_type), POINTER :: loc_results, results_mix
1484 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1485 TYPE(cp_subsys_type), POINTER :: subsys_mix
1486 TYPE(particle_list_type), POINTER :: particles_mix
1487 TYPE(section_vals_type), POINTER :: force_env_section, mapping_section, &
1488 mixed_section, root_section
1489 TYPE(virial_type), POINTER :: loc_virial, virial_mix
1490
1491 logger => cp_get_default_logger()
1492 cpassert(ASSOCIATED(force_env))
1493 ! Get infos about the mixed subsys
1494 CALL force_env_get(force_env=force_env, &
1495 subsys=subsys_mix, &
1496 force_env_section=force_env_section, &
1497 root_section=root_section, &
1498 cell=cell_mix)
1499 CALL cp_subsys_get(subsys=subsys_mix, &
1500 particles=particles_mix, &
1501 virial=virial_mix, &
1502 results=results_mix)
1503 NULLIFY (map_index)
1504 nforce_eval = SIZE(force_env%sub_force_env)
1505 mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
1506 mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
1507 ALLOCATE (subsystems(nforce_eval))
1508 DO iforce_eval = 1, nforce_eval
1509 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1510 NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
1511 ALLOCATE (virials(iforce_eval)%virial)
1512 CALL cp_result_create(results(iforce_eval)%results)
1513 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1514 ! Get all available subsys
1515 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1516 subsys=subsystems(iforce_eval)%subsys)
1517
1518 ! all force_env share the same cell
1519 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
1520
1521 ! Get available particles
1522 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1523 particles=particles(iforce_eval)%list)
1524
1525 ! Get Mapping index array
1526 natom = SIZE(particles(iforce_eval)%list%els)
1527 ! Serial mode need to deallocate first
1528 IF (ASSOCIATED(map_index)) THEN
1529 DEALLOCATE (map_index)
1530 END IF
1531 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1532 map_index)
1533
1534 ! Mapping particles from iforce_eval environment to the mixed env
1535 DO iparticle = 1, natom
1536 jparticle = map_index(iparticle)
1537 particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
1538 END DO
1539 ! Mixed CDFT + QMMM: Need to translate now
1540 IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1541 CALL apply_qmmm_translate(force_env%sub_force_env(iforce_eval)%force_env%qmmm_env)
1542 END IF
1543 END DO
1544 ! For mixed CDFT calculations parallelized over CDFT states
1545 ! build weight and gradient on all processors before splitting into groups and
1546 ! starting energy calculation
1547 CALL mixed_cdft_build_weight(force_env, calculate_forces)
1548 DO iforce_eval = 1, nforce_eval
1549 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1550 ! From this point on the error is the sub_error
1551 IF (force_env%mixed_env%cdft_control%run_type == mixed_cdft_serial .AND. iforce_eval >= 2) THEN
1552 my_logger => force_env%mixed_env%cdft_control%sub_logger(iforce_eval - 1)%p
1553 ELSE
1554 my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
1555 my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
1556 END IF
1557 ! Copy iterations info (they are updated only in the main mixed_env)
1558 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1559 CALL cp_add_default_logger(my_logger)
1560 ! Serial CDFT calculation: transfer weight/gradient
1561 CALL mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
1562 ! Calculate energy and forces for each sub_force_env
1563 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1564 calc_force=calculate_forces, &
1565 skip_external_control=.true.)
1566 ! Only the rank 0 process collect info for each computation
1567 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1568 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1569 potential_energy=energy)
1570 CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1571 virial=loc_virial, results=loc_results)
1572 energies(iforce_eval) = energy
1573 glob_natoms(iforce_eval) = natom
1574 virials(iforce_eval)%virial = loc_virial
1575 CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1576 END IF
1577 ! Deallocate map_index array
1578 IF (ASSOCIATED(map_index)) THEN
1579 DEALLOCATE (map_index)
1580 END IF
1581 CALL cp_rm_default_logger()
1582 END DO
1583 DEALLOCATE (subsystems)
1584
1585 END SUBROUTINE mixed_cdft_energy_forces
1586
1587! **************************************************************************************************
1588!> \brief Perform additional tasks for mixed CDFT calculations after solving the electronic structure
1589!> of both CDFT states
1590!> \param force_env the force_env that holds the CDFT states
1591!> \par History
1592!> 01.17 created [Nico Holmberg]
1593!> \author Nico Holmberg
1594! **************************************************************************************************
1595 SUBROUTINE mixed_cdft_post_energy_forces(force_env)
1596 TYPE(force_env_type), POINTER :: force_env
1597
1598 INTEGER :: iforce_eval, nforce_eval, nvar
1599 TYPE(dft_control_type), POINTER :: dft_control
1600 TYPE(qs_environment_type), POINTER :: qs_env
1601
1602 cpassert(ASSOCIATED(force_env))
1603 NULLIFY (qs_env, dft_control)
1604 IF (force_env%mixed_env%do_mixed_cdft) THEN
1605 nforce_eval = SIZE(force_env%sub_force_env)
1606 nvar = force_env%mixed_env%cdft_control%nconstraint
1607 ! Transfer cdft strengths for writing restart
1608 IF (.NOT. ASSOCIATED(force_env%mixed_env%strength)) THEN
1609 ALLOCATE (force_env%mixed_env%strength(nforce_eval, nvar))
1610 END IF
1611 force_env%mixed_env%strength = 0.0_dp
1612 DO iforce_eval = 1, nforce_eval
1613 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1614 IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1615 qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
1616 ELSE
1617 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1618 END IF
1619 CALL get_qs_env(qs_env, dft_control=dft_control)
1620 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1621 force_env%mixed_env%strength(iforce_eval, :) = dft_control%qs_control%cdft_control%strength(:)
1622 END IF
1623 END DO
1624 CALL force_env%para_env%sum(force_env%mixed_env%strength)
1625 ! Mixed CDFT: calculate ET coupling
1626 IF (force_env%mixed_env%do_mixed_et) THEN
1627 IF (modulo(force_env%mixed_env%cdft_control%sim_step, force_env%mixed_env%et_freq) == 0) THEN
1628 CALL mixed_cdft_calculate_coupling(force_env)
1629 END IF
1630 END IF
1631 END IF
1632
1633 END SUBROUTINE mixed_cdft_post_energy_forces
1634
1635! **************************************************************************************************
1636!> \brief Computes the total energy for an embedded calculation
1637!> \param force_env ...
1638!> \author Vladimir Rybkin
1639! **************************************************************************************************
1640 SUBROUTINE embed_energy(force_env)
1641
1642 TYPE(force_env_type), POINTER :: force_env
1643
1644 INTEGER :: iforce_eval, iparticle, jparticle, &
1645 my_group, natom, nforce_eval
1646 INTEGER, DIMENSION(:), POINTER :: glob_natoms, map_index
1647 LOGICAL :: converged_embed
1648 REAL(kind=dp) :: energy
1649 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1650 TYPE(cell_type), POINTER :: cell_embed
1651 TYPE(cp_logger_type), POINTER :: logger, my_logger
1652 TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1653 TYPE(cp_result_type), POINTER :: loc_results, results_embed
1654 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1655 TYPE(cp_subsys_type), POINTER :: subsys_embed
1656 TYPE(dft_control_type), POINTER :: dft_control
1657 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1658 TYPE(particle_list_type), POINTER :: particles_embed
1659 TYPE(pw_env_type), POINTER :: pw_env
1660 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1661 TYPE(pw_r3d_rs_type), POINTER :: embed_pot, spin_embed_pot
1662 TYPE(section_vals_type), POINTER :: embed_section, force_env_section, &
1663 mapping_section, root_section
1664
1665 logger => cp_get_default_logger()
1666 cpassert(ASSOCIATED(force_env))
1667 ! Get infos about the embedding subsys
1668 CALL force_env_get(force_env=force_env, &
1669 subsys=subsys_embed, &
1670 force_env_section=force_env_section, &
1671 root_section=root_section, &
1672 cell=cell_embed)
1673 CALL cp_subsys_get(subsys=subsys_embed, &
1674 particles=particles_embed, &
1675 results=results_embed)
1676 NULLIFY (map_index, glob_natoms)
1677
1678 nforce_eval = SIZE(force_env%sub_force_env)
1679 embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1680 mapping_section => section_vals_get_subs_vals(embed_section, "MAPPING")
1681 ! Global Info
1682 ALLOCATE (subsystems(nforce_eval))
1683 ALLOCATE (particles(nforce_eval))
1684 ! Local Info to sync
1685 ALLOCATE (energies(nforce_eval))
1686 ALLOCATE (glob_natoms(nforce_eval))
1687 ALLOCATE (results(nforce_eval))
1688 energies = 0.0_dp
1689 glob_natoms = 0
1690
1691 DO iforce_eval = 1, nforce_eval
1692 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1693 NULLIFY (results(iforce_eval)%results)
1694 CALL cp_result_create(results(iforce_eval)%results)
1695 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1696 ! From this point on the error is the sub_error
1697 my_group = force_env%embed_env%group_distribution(force_env%para_env%mepos)
1698 my_logger => force_env%embed_env%sub_logger(my_group + 1)%p
1699 ! Copy iterations info (they are updated only in the main embed_env)
1700 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1701 CALL cp_add_default_logger(my_logger)
1702
1703 ! Get all available subsys
1704 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1705 subsys=subsystems(iforce_eval)%subsys)
1706
1707 ! Check if we import density from previous force calculations
1708 ! Only for QUICKSTEP
1709 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env%qs_env)) THEN
1710 NULLIFY (dft_control)
1711 CALL get_qs_env(force_env%sub_force_env(iforce_eval)%force_env%qs_env, dft_control=dft_control)
1712 IF (dft_control%qs_control%ref_embed_subsys) THEN
1713 IF (iforce_eval == 2) cpabort("Density importing force_eval can't be the first.")
1714 END IF
1715 END IF
1716
1717 ! all force_env share the same cell
1718 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_embed)
1719
1720 ! Get available particles
1721 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1722 particles=particles(iforce_eval)%list)
1723
1724 ! Get Mapping index array
1725 natom = SIZE(particles(iforce_eval)%list%els)
1726
1727 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1728 map_index, .true.)
1729
1730 ! Mapping particles from iforce_eval environment to the embed env
1731 DO iparticle = 1, natom
1732 jparticle = map_index(iparticle)
1733 particles(iforce_eval)%list%els(iparticle)%r = particles_embed%els(jparticle)%r
1734 END DO
1735
1736 ! Calculate energy and forces for each sub_force_env
1737 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1738 calc_force=.false., &
1739 skip_external_control=.true.)
1740
1741 ! Call DFT embedding
1742 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env%qs_env)) THEN
1743 NULLIFY (dft_control)
1744 CALL get_qs_env(force_env%sub_force_env(iforce_eval)%force_env%qs_env, dft_control=dft_control)
1745 IF (dft_control%qs_control%ref_embed_subsys) THEN
1746 ! Now we can optimize the embedding potential
1747 CALL dft_embedding(force_env, iforce_eval, energies, converged_embed)
1748 IF (.NOT. converged_embed) cpabort("Embedding potential optimization not converged.")
1749 END IF
1750 ! Deallocate embedding potential on the high-level subsystem
1751 IF (dft_control%qs_control%high_level_embed_subsys) THEN
1752 CALL get_qs_env(qs_env=force_env%sub_force_env(iforce_eval)%force_env%qs_env, &
1753 embed_pot=embed_pot, spin_embed_pot=spin_embed_pot, pw_env=pw_env)
1754 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1755 CALL auxbas_pw_pool%give_back_pw(embed_pot)
1756 IF (ASSOCIATED(embed_pot)) THEN
1757 CALL embed_pot%release()
1758 DEALLOCATE (embed_pot)
1759 END IF
1760 IF (ASSOCIATED(spin_embed_pot)) THEN
1761 CALL auxbas_pw_pool%give_back_pw(spin_embed_pot)
1762 CALL spin_embed_pot%release()
1763 DEALLOCATE (spin_embed_pot)
1764 END IF
1765 END IF
1766 END IF
1767
1768 ! Only the rank 0 process collect info for each computation
1769 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1770 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1771 potential_energy=energy)
1772 CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1773 results=loc_results)
1774 energies(iforce_eval) = energy
1775 glob_natoms(iforce_eval) = natom
1776 CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1777 END IF
1778 ! Deallocate map_index array
1779 IF (ASSOCIATED(map_index)) THEN
1780 DEALLOCATE (map_index)
1781 END IF
1782 CALL cp_rm_default_logger()
1783 END DO
1784
1785 ! Handling Parallel execution
1786 CALL force_env%para_env%sync()
1787 ! Let's transfer energy, natom
1788 CALL force_env%para_env%sum(energies)
1789 CALL force_env%para_env%sum(glob_natoms)
1790
1791 force_env%embed_env%energies = energies
1792
1793 !NB if the first system has fewer atoms than the second)
1794 DO iparticle = 1, SIZE(particles_embed%els)
1795 particles_embed%els(iparticle)%f(:) = 0.0_dp
1796 END DO
1797
1798 ! ONIOM type of mixing in embedding: E = E_total + E_cluster_high - E_cluster
1799 force_env%embed_env%pot_energy = energies(3) + energies(4) - energies(2)
1800
1801 !Simply deallocate and loose the pointer references..
1802 DO iforce_eval = 1, nforce_eval
1803 CALL cp_result_release(results(iforce_eval)%results)
1804 END DO
1805 DEALLOCATE (subsystems)
1806 DEALLOCATE (particles)
1807 DEALLOCATE (energies)
1808 DEALLOCATE (glob_natoms)
1809 DEALLOCATE (results)
1810
1811 END SUBROUTINE embed_energy
1812
1813! **************************************************************************************************
1814!> \brief ...
1815!> \param force_env ...
1816!> \param ref_subsys_number ...
1817!> \param energies ...
1818!> \param converged_embed ...
1819! **************************************************************************************************
1820 SUBROUTINE dft_embedding(force_env, ref_subsys_number, energies, converged_embed)
1821 TYPE(force_env_type), POINTER :: force_env
1822 INTEGER :: ref_subsys_number
1823 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1824 LOGICAL :: converged_embed
1825
1826 INTEGER :: embed_method
1827 TYPE(section_vals_type), POINTER :: embed_section, force_env_section
1828
1829 ! Find out which embedding scheme is used
1830 CALL force_env_get(force_env=force_env, &
1831 force_env_section=force_env_section)
1832 embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1833
1834 CALL section_vals_val_get(embed_section, "EMBED_METHOD", i_val=embed_method)
1835 SELECT CASE (embed_method)
1836 CASE (dfet)
1837 ! Density functional embedding
1838 CALL dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1839 CASE (dmfet)
1840 ! Density matrix embedding theory
1841 CALL dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1842 END SELECT
1843
1844 END SUBROUTINE dft_embedding
1845! **************************************************************************************************
1846!> \brief ... Main driver for DFT embedding
1847!> \param force_env ...
1848!> \param ref_subsys_number ...
1849!> \param energies ...
1850!> \param converged_embed ...
1851!> \author Vladimir Rybkin
1852! **************************************************************************************************
1853 SUBROUTINE dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1854 TYPE(force_env_type), POINTER :: force_env
1855 INTEGER :: ref_subsys_number
1856 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1857 LOGICAL :: converged_embed
1858
1859 CHARACTER(LEN=*), PARAMETER :: routinen = 'dfet_embedding'
1860
1861 INTEGER :: cluster_subsys_num, handle, &
1862 i_force_eval, i_iter, i_spin, &
1863 nforce_eval, nspins, nspins_subsys, &
1864 output_unit
1865 REAL(kind=dp) :: cluster_energy
1866 REAL(kind=dp), DIMENSION(:), POINTER :: rhs
1867 TYPE(cp_logger_type), POINTER :: logger
1868 TYPE(dft_control_type), POINTER :: dft_control
1869 TYPE(opt_embed_pot_type) :: opt_embed
1870 TYPE(pw_env_type), POINTER :: pw_env
1871 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1872 TYPE(pw_r3d_rs_type) :: diff_rho_r, diff_rho_spin
1873 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_ref, rho_r_subsys
1874 TYPE(pw_r3d_rs_type), POINTER :: embed_pot, embed_pot_subsys, &
1875 spin_embed_pot, spin_embed_pot_subsys
1876 TYPE(qs_energy_type), POINTER :: energy
1877 TYPE(qs_rho_type), POINTER :: rho, subsys_rho
1878 TYPE(section_vals_type), POINTER :: dft_section, embed_section, &
1879 force_env_section, input, &
1880 mapping_section, opt_embed_section
1881
1882 CALL timeset(routinen, handle)
1883
1884 CALL cite_reference(huang2011)
1885 CALL cite_reference(heaton_burgess2007)
1886
1887 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
1888
1889 ! Reveal input file
1890 NULLIFY (logger)
1891 logger => cp_get_default_logger()
1892 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
1893 extension=".Log")
1894
1895 NULLIFY (dft_section, input, opt_embed_section)
1896 NULLIFY (energy, dft_control)
1897
1898 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1899 pw_env=pw_env, dft_control=dft_control, rho=rho, energy=energy, &
1900 input=input)
1901 nspins = dft_control%nspins
1902
1903 dft_section => section_vals_get_subs_vals(input, "DFT")
1904 opt_embed_section => section_vals_get_subs_vals(input, &
1905 "DFT%QS%OPT_EMBED")
1906 ! Rho_r is the reference
1907 CALL qs_rho_get(rho_struct=rho, rho_r=rho_r_ref)
1908
1909 ! We need to understand how to treat spins states
1910 CALL understand_spin_states(force_env, ref_subsys_number, opt_embed%change_spin, opt_embed%open_shell_embed, &
1911 opt_embed%all_nspins)
1912
1913 ! Prepare everything for the potential maximization
1914 CALL prepare_embed_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, opt_embed, &
1915 opt_embed_section)
1916
1917 ! Initialize embedding potential
1918 CALL init_embed_pot(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, &
1919 opt_embed%add_const_pot, opt_embed%Fermi_Amaldi, opt_embed%const_pot, &
1920 opt_embed%open_shell_embed, spin_embed_pot, &
1921 opt_embed%pot_diff, opt_embed%Coulomb_guess, opt_embed%grid_opt)
1922
1923 ! Read embedding potential vector from the file
1924 IF (opt_embed%read_embed_pot .OR. opt_embed%read_embed_pot_cube) CALL read_embed_pot( &
1925 force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, spin_embed_pot, &
1926 opt_embed_section, opt_embed)
1927
1928 ! Prepare the pw object to store density differences
1929 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1930 CALL auxbas_pw_pool%create_pw(diff_rho_r)
1931 CALL pw_zero(diff_rho_r)
1932 IF (opt_embed%open_shell_embed) THEN
1933 CALL auxbas_pw_pool%create_pw(diff_rho_spin)
1934 CALL pw_zero(diff_rho_spin)
1935 END IF
1936
1937 ! Check the preliminary density differences
1938 DO i_spin = 1, nspins
1939 CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
1940 END DO
1941 IF (opt_embed%open_shell_embed) THEN ! Spin part
1942 IF (nspins == 2) THEN ! Reference systems has an open shell, else the reference diff_rho_spin is zero
1943 CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
1944 CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
1945 END IF
1946 END IF
1947
1948 DO i_force_eval = 1, ref_subsys_number - 1
1949 NULLIFY (subsys_rho, rho_r_subsys, dft_control)
1950 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, energy=energy, &
1951 dft_control=dft_control)
1952 nspins_subsys = dft_control%nspins
1953 ! Add subsystem densities
1954 CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
1955 DO i_spin = 1, nspins_subsys
1956 CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
1957 END DO
1958 IF (opt_embed%open_shell_embed) THEN ! Spin part
1959 IF (nspins_subsys == 2) THEN ! The subsystem makes contribution if it is spin-polarized
1960 ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
1961 IF ((i_force_eval == 2) .AND. (opt_embed%change_spin)) THEN
1962 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1963 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1964 ELSE
1965 ! First subsystem (always) and second subsystem (without spin change)
1966 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1967 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1968 END IF
1969 END IF
1970 END IF
1971 END DO
1972
1973 ! Print density difference
1974 CALL print_rho_diff(diff_rho_r, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1975 IF (opt_embed%open_shell_embed) THEN ! Spin part
1976 CALL print_rho_spin_diff(diff_rho_spin, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1977 END IF
1978
1979 ! Construct electrostatic guess if needed
1980 IF (opt_embed%Coulomb_guess) THEN
1981 ! Reveal resp charges for total system
1982 nforce_eval = SIZE(force_env%sub_force_env)
1983 NULLIFY (rhs)
1984 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, rhs=rhs)
1985 ! Get the mapping
1986 CALL force_env_get(force_env=force_env, &
1987 force_env_section=force_env_section)
1988 embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1989 mapping_section => section_vals_get_subs_vals(embed_section, "MAPPING")
1990
1991 DO i_force_eval = 1, ref_subsys_number - 1
1992 IF (i_force_eval == 1) THEN
1993 CALL coulomb_guess(embed_pot, rhs, mapping_section, &
1994 force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1995 ELSE
1996 CALL coulomb_guess(opt_embed%pot_diff, rhs, mapping_section, &
1997 force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1998 END IF
1999 END DO
2000 CALL pw_axpy(opt_embed%pot_diff, embed_pot)
2001 IF (.NOT. opt_embed%grid_opt) CALL pw_copy(embed_pot, opt_embed%const_pot)
2002
2003 END IF
2004
2005 ! Difference guess
2006 IF (opt_embed%diff_guess) THEN
2007 CALL pw_copy(diff_rho_r, embed_pot)
2008 IF (.NOT. opt_embed%grid_opt) CALL pw_copy(embed_pot, opt_embed%const_pot)
2009 ! Open shell
2010 IF (opt_embed%open_shell_embed) CALL pw_copy(diff_rho_spin, spin_embed_pot)
2011 END IF
2012
2013 ! Calculate subsystems with trial embedding potential
2014 DO i_iter = 1, opt_embed%n_iter
2015 opt_embed%i_iter = i_iter
2016
2017 ! Set the density difference as the negative reference one
2018 CALL pw_zero(diff_rho_r)
2019 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, dft_control=dft_control)
2020 nspins = dft_control%nspins
2021 DO i_spin = 1, nspins
2022 CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
2023 END DO
2024 IF (opt_embed%open_shell_embed) THEN ! Spin part
2025 CALL pw_zero(diff_rho_spin)
2026 IF (nspins == 2) THEN ! Reference systems has an open shell, else the reference diff_rho_spin is zero
2027 CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
2028 CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
2029 END IF
2030 END IF
2031
2032 DO i_force_eval = 1, ref_subsys_number - 1
2033 NULLIFY (dft_control)
2034 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
2035 nspins_subsys = dft_control%nspins
2036
2037 IF ((i_force_eval == 2) .AND. (opt_embed%change_spin)) THEN
2038 ! Here we change the sign of the spin embedding potential due to spin change:
2039 ! only in spin_embed_subsys
2040 CALL make_subsys_embed_pot(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2041 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
2042 opt_embed%open_shell_embed, .true.)
2043 ELSE ! Regular case
2044 CALL make_subsys_embed_pot(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2045 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
2046 opt_embed%open_shell_embed, .false.)
2047 END IF
2048
2049 ! Switch on external potential in the subsystems
2050 dft_control%apply_embed_pot = .true.
2051
2052 ! Add the embedding potential
2053 CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, embed_pot=embed_pot_subsys)
2054 IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys == 2)) THEN ! Spin part
2055 CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2056 spin_embed_pot=spin_embed_pot_subsys)
2057 END IF
2058
2059 ! Get the previous subsystem densities
2060 CALL get_prev_density(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
2061
2062 ! Calculate the new density
2063 CALL force_env_calc_energy_force(force_env=force_env%sub_force_env(i_force_eval)%force_env, &
2064 calc_force=.false., &
2065 skip_external_control=.true.)
2066
2067 CALL get_max_subsys_diff(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
2068
2069 ! Extract subsystem density and energy
2070 NULLIFY (rho_r_subsys, energy)
2071
2072 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, &
2073 energy=energy)
2074 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + energy%total
2075
2076 ! Find out which subsystem is the cluster
2077 IF (dft_control%qs_control%cluster_embed_subsys) THEN
2078 cluster_subsys_num = i_force_eval
2079 cluster_energy = energy%total
2080 END IF
2081
2082 ! Add subsystem densities
2083 CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
2084 DO i_spin = 1, nspins_subsys
2085 CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
2086 END DO
2087 IF (opt_embed%open_shell_embed) THEN ! Spin part
2088 IF (nspins_subsys == 2) THEN ! The subsystem makes contribution if it is spin-polarized
2089 ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
2090 IF ((i_force_eval == 2) .AND. (opt_embed%change_spin)) THEN
2091 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
2092 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
2093 ELSE
2094 ! First subsystem (always) and second subsystem (without spin change)
2095 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
2096 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
2097 END IF
2098 END IF
2099 END IF
2100
2101 ! Release embedding potential for subsystem
2102 CALL embed_pot_subsys%release()
2103 DEALLOCATE (embed_pot_subsys)
2104 IF (opt_embed%open_shell_embed) THEN
2105 CALL spin_embed_pot_subsys%release()
2106 DEALLOCATE (spin_embed_pot_subsys)
2107 END IF
2108
2109 END DO ! i_force_eval
2110
2111 ! Print embedding potential for restart
2112 CALL print_embed_restart(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2113 opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
2114 spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .false.)
2115 CALL print_pot_simple_grid(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2116 embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .false., &
2117 force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
2118
2119 ! Integrate the potential over density differences and add to w functional; also add regularization contribution
2120 DO i_spin = 1, nspins ! Sum over nspins for the reference system, not subsystem!
2121 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) - pw_integral_ab(embed_pot, rho_r_ref(i_spin))
2122 END DO
2123 ! Spin part
2124 IF (opt_embed%open_shell_embed) THEN
2125 ! If reference system is not spin-polarized then it does not make a contribution to W functional
2126 IF (nspins == 2) THEN
2127 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) &
2128 - pw_integral_ab(spin_embed_pot, rho_r_ref(1)) &
2129 + pw_integral_ab(spin_embed_pot, rho_r_ref(2))
2130 END IF
2131 END IF
2132 ! Finally, add the regularization term
2133 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + opt_embed%reg_term
2134
2135 ! Print information and check convergence
2136 CALL print_emb_opt_info(output_unit, i_iter, opt_embed)
2137 CALL conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit)
2138 IF (opt_embed%converged) EXIT
2139
2140 ! Update the trust radius and control the step
2141 IF ((i_iter > 1) .AND. (.NOT. opt_embed%steep_desc)) CALL step_control(opt_embed)
2142
2143 ! Print density difference
2144 CALL print_rho_diff(diff_rho_r, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
2145 IF (opt_embed%open_shell_embed) THEN ! Spin part
2146 CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
2147 END IF
2148
2149 ! Calculate potential gradient if the step has been accepted. Otherwise, we reuse the previous one
2150
2151 IF (opt_embed%accept_step .AND. (.NOT. opt_embed%grid_opt)) THEN
2152 CALL calculate_embed_pot_grad(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2153 diff_rho_r, diff_rho_spin, opt_embed)
2154 END IF
2155 ! Take the embedding step
2156 CALL opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, &
2157 force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2158
2159 END DO ! i_iter
2160
2161 ! Print final embedding potential for restart
2162 CALL print_embed_restart(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2163 opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
2164 spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .true.)
2165 CALL print_pot_simple_grid(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2166 embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .true., &
2167 force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
2168
2169 ! Print final density difference
2170 !CALL print_rho_diff(diff_rho_r, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .TRUE.)
2171 IF (opt_embed%open_shell_embed) THEN ! Spin part
2172 CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .true.)
2173 END IF
2174
2175 ! Give away plane waves pools
2176 CALL diff_rho_r%release()
2177 IF (opt_embed%open_shell_embed) THEN
2178 CALL diff_rho_spin%release()
2179 END IF
2180
2181 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
2182 "PRINT%PROGRAM_RUN_INFO")
2183
2184 ! If converged send the embedding potential to the higher-level calculation.
2185 IF (opt_embed%converged) THEN
2186 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, dft_control=dft_control, &
2187 pw_env=pw_env)
2188 nspins_subsys = dft_control%nspins
2189 dft_control%apply_embed_pot = .true.
2190 ! The embedded subsystem corresponds to subsystem #2, where spin change is possible
2191 CALL make_subsys_embed_pot(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, &
2192 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
2193 opt_embed%open_shell_embed, opt_embed%change_spin)
2194
2195 IF (opt_embed%Coulomb_guess) THEN
2196 CALL pw_axpy(opt_embed%pot_diff, embed_pot_subsys, -1.0_dp, allow_noncompatible_grids=.true.)
2197 END IF
2198
2199 CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, embed_pot=embed_pot_subsys)
2200
2201 IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys == 2)) THEN
2202 CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, &
2203 spin_embed_pot=spin_embed_pot_subsys)
2204 END IF
2205
2206 ! Substitute the correct energy in energies: only on rank 0
2207 IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source()) THEN
2208 energies(cluster_subsys_num) = cluster_energy
2209 END IF
2210 END IF
2211
2212 ! Deallocate and release opt_embed content
2213 CALL release_opt_embed(opt_embed)
2214
2215 ! Deallocate embedding potential
2216 CALL embed_pot%release()
2217 DEALLOCATE (embed_pot)
2218 IF (opt_embed%open_shell_embed) THEN
2219 CALL spin_embed_pot%release()
2220 DEALLOCATE (spin_embed_pot)
2221 END IF
2222
2223 converged_embed = opt_embed%converged
2224
2225 CALL timestop(handle)
2226
2227 END SUBROUTINE dfet_embedding
2228
2229! **************************************************************************************************
2230!> \brief Main driver for the DMFET embedding
2231!> \param force_env ...
2232!> \param ref_subsys_number ...
2233!> \param energies ...
2234!> \param converged_embed ...
2235!> \author Vladimir Rybkin
2236! **************************************************************************************************
2237 SUBROUTINE dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
2238 TYPE(force_env_type), POINTER :: force_env
2239 INTEGER :: ref_subsys_number
2240 REAL(kind=dp), DIMENSION(:), POINTER :: energies
2241 LOGICAL :: converged_embed
2242
2243 CHARACTER(LEN=*), PARAMETER :: routinen = 'dmfet_embedding'
2244
2245 INTEGER :: cluster_subsys_num, handle, &
2246 i_force_eval, i_iter, output_unit
2247 LOGICAL :: subsys_open_shell
2248 REAL(kind=dp) :: cluster_energy
2249 TYPE(cp_logger_type), POINTER :: logger
2250 TYPE(dft_control_type), POINTER :: dft_control
2251 TYPE(mp_para_env_type), POINTER :: para_env
2252 TYPE(opt_dmfet_pot_type) :: opt_dmfet
2253 TYPE(qs_energy_type), POINTER :: energy
2254 TYPE(section_vals_type), POINTER :: dft_section, input, opt_dmfet_section
2255
2256 CALL timeset(routinen, handle)
2257
2258 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2259 para_env=para_env)
2260
2261 ! Reveal input file
2262 NULLIFY (logger)
2263 logger => cp_get_default_logger()
2264 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
2265 extension=".Log")
2266
2267 NULLIFY (dft_section, input, opt_dmfet_section)
2268 NULLIFY (energy)
2269
2270 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2271 energy=energy, input=input)
2272
2273 dft_section => section_vals_get_subs_vals(input, "DFT")
2274 opt_dmfet_section => section_vals_get_subs_vals(input, &
2275 "DFT%QS%OPT_DMFET")
2276
2277 ! We need to understand how to treat spins states
2278 CALL understand_spin_states(force_env, ref_subsys_number, opt_dmfet%change_spin, opt_dmfet%open_shell_embed, &
2279 opt_dmfet%all_nspins)
2280
2281 ! Prepare for the potential optimization
2282 CALL prepare_dmfet_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2283 opt_dmfet, opt_dmfet_section)
2284
2285 ! Get the reference density matrix/matrices
2286 subsys_open_shell = subsys_spin(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2287 CALL build_full_dm(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2288 opt_dmfet%dm_total, subsys_open_shell, opt_dmfet%open_shell_embed, opt_dmfet%dm_total_beta)
2289
2290 ! Check the preliminary DM difference
2291 CALL cp_fm_copy_general(opt_dmfet%dm_total, opt_dmfet%dm_diff, para_env)
2292 IF (opt_dmfet%open_shell_embed) CALL cp_fm_copy_general(opt_dmfet%dm_total_beta, &
2293 opt_dmfet%dm_diff_beta, para_env)
2294
2295 DO i_force_eval = 1, ref_subsys_number - 1
2296
2297 ! Get the subsystem density matrix/matrices
2298 subsys_open_shell = subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2299
2300 CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2301 opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2302 opt_dmfet%dm_subsys_beta)
2303
2304 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2305
2306 IF (opt_dmfet%open_shell_embed) CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, &
2307 1.0_dp, opt_dmfet%dm_subsys_beta)
2308
2309 END DO
2310
2311 ! Main loop of iterative matrix potential optimization
2312 DO i_iter = 1, opt_dmfet%n_iter
2313
2314 opt_dmfet%i_iter = i_iter
2315
2316 ! Set the dm difference as the reference one
2317 CALL cp_fm_copy_general(opt_dmfet%dm_total, opt_dmfet%dm_diff, para_env)
2318
2319 IF (opt_dmfet%open_shell_embed) CALL cp_fm_copy_general(opt_dmfet%dm_total_beta, &
2320 opt_dmfet%dm_diff_beta, para_env)
2321
2322 ! Loop over force evaluations
2323 DO i_force_eval = 1, ref_subsys_number - 1
2324
2325 ! Switch on external potential in the subsystems
2326 NULLIFY (dft_control)
2327 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
2328 dft_control%apply_dmfet_pot = .true.
2329
2330 ! Calculate the new density
2331 CALL force_env_calc_energy_force(force_env=force_env%sub_force_env(i_force_eval)%force_env, &
2332 calc_force=.false., &
2333 skip_external_control=.true.)
2334
2335 ! Extract subsystem density matrix and energy
2336 NULLIFY (energy)
2337
2338 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, energy=energy)
2339 opt_dmfet%w_func(i_iter) = opt_dmfet%w_func(i_iter) + energy%total
2340
2341 ! Find out which subsystem is the cluster
2342 IF (dft_control%qs_control%cluster_embed_subsys) THEN
2343 cluster_subsys_num = i_force_eval
2344 cluster_energy = energy%total
2345 END IF
2346
2347 ! Add subsystem density matrices
2348 subsys_open_shell = subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2349
2350 CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2351 opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2352 opt_dmfet%dm_subsys_beta)
2353
2354 IF (opt_dmfet%open_shell_embed) THEN ! Open-shell embedding
2355 ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
2356 IF ((i_force_eval == 2) .AND. (opt_dmfet%change_spin)) THEN
2357 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, 1.0_dp, opt_dmfet%dm_subsys)
2358 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys_beta)
2359 ELSE
2360 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2361 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, 1.0_dp, opt_dmfet%dm_subsys_beta)
2362 END IF
2363 ELSE ! Closed-shell embedding
2364 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2365 END IF
2366
2367 END DO ! i_force_eval
2368
2369 CALL check_dmfet(opt_dmfet, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2370
2371 END DO ! i_iter
2372
2373 ! Substitute the correct energy in energies: only on rank 0
2374 IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source()) THEN
2375 energies(cluster_subsys_num) = cluster_energy
2376 END IF
2377
2378 CALL release_dmfet_opt(opt_dmfet)
2379
2380 converged_embed = .false.
2381
2382 END SUBROUTINE dmfet_embedding
2383
2384END MODULE force_env_methods
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
cp_result_methods::get_results
Definition cp_result_methods.F:43
kahan_sum::accurate_sum
Definition kahan_sum.F:41
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_zero
Definition pw_methods.F:83
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
atprop_types::atprop_init
subroutine, public atprop_init(atprop_env, natom)
...
Definition atprop_types.F:64
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::huang2011
integer, save, public huang2011
Definition bibliography.F:36
bibliography::heaton_burgess2007
integer, save, public heaton_burgess2007
Definition bibliography.F:36
cell_methods
Handles all functions related to the CELL.
Definition cell_methods.F:15
cell_methods::init_cell
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
Definition cell_methods.F:185
cell_methods::cell_create
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition cell_methods.F:95
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::scaled_to_real
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:661
cell_types::cell_sym_triclinic
integer, parameter, public cell_sym_triclinic
Definition cell_types.F:29
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:631
cell_types::cell_release
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:704
cell_types::cell_clone
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition cell_types.F:142
constraint_fxd
Definition constraint_fxd.F:12
constraint_fxd::fix_atom_control
subroutine, public fix_atom_control(force_env, w)
allows for fix atom constraints
Definition constraint_fxd.F:65
constraint_vsite
Routines to handle the virtual site constraint/restraint.
Definition constraint_vsite.F:14
constraint_vsite::vsite_force_control
subroutine, public vsite_force_control(force_env)
control force distribution for virtual sites
Definition constraint_vsite.F:53
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:171
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition cp_fm_types.F:1587
cp_iter_types
Collection of routines to handle the iteration info.
Definition cp_iter_types.F:11
cp_iter_types::cp_iteration_info_copy_iter
subroutine, public cp_iteration_info_copy_iter(iteration_info_in, iteration_info_out)
Copies iterations info of an iteration info into another iteration info.
Definition cp_iter_types.F:176
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_rm_default_logger
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
Definition cp_log_handling.F:213
cp_log_handling::cp_add_default_logger
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
Definition cp_log_handling.F:188
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:77
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1156
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:107
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout, file_access)
...
Definition cp_output_handling.F:935
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition cp_result_methods.F:17
cp_result_methods::cp_results_mp_bcast
subroutine, public cp_results_mp_bcast(results, source, para_env)
broadcast results type
Definition cp_result_methods.F:437
cp_result_methods::cp_results_erase
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
Definition cp_result_methods.F:365
cp_result_methods::test_for_result
logical function, public test_for_result(results, description)
test for a certain result in the result_list
Definition cp_result_methods.F:127
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition cp_result_types.F:18
cp_result_types::cp_result_copy
subroutine, public cp_result_copy(results_in, results_out)
Copies the cp_result type.
Definition cp_result_types.F:284
cp_result_types::cp_result_release
subroutine, public cp_result_release(results)
Releases cp_result type.
Definition cp_result_types.F:120
cp_result_types::cp_result_create
subroutine, public cp_result_create(results)
Allocates and intitializes the cp_result.
Definition cp_result_types.F:96
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_set
subroutine, public cp_subsys_set(subsys, atomic_kinds, particles, local_particles, molecules, molecule_kinds, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, results, cell, cell_ref, use_ref_cell)
sets various propreties of the subsys
Definition cp_subsys_types.F:210
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1311
eip_environment_types
The environment for the empirical interatomic potential methods.
Definition eip_environment_types.F:14
eip_silicon
Empirical interatomic potentials for Silicon.
Definition eip_silicon.F:17
eip_silicon::eip_stillinger_weber
subroutine, public eip_stillinger_weber(eip_env)
Interface routine of the Stillinger-Weber force field to CP2K.
Definition eip_silicon.F:407
eip_silicon::eip_lenosky
subroutine, public eip_lenosky(eip_env)
Interface routine of Goedecker's Lenosky force field to CP2K.
Definition eip_silicon.F:243
eip_silicon::eip_tersoff
subroutine, public eip_tersoff(eip_env)
Interface routine of the Tersoff force field to CP2K.
Definition eip_silicon.F:569
eip_silicon::eip_bazant
subroutine, public eip_bazant(eip_env)
Interface routine of Goedecker's Bazant EDIP to CP2K.
Definition eip_silicon.F:77
embed_types
Definition embed_types.F:8
external_potential_methods
Methods to include the effect of an external potential during an MD or energy calculation.
Definition external_potential_methods.F:13
external_potential_methods::add_external_potential
subroutine, public add_external_potential(force_env)
...
Definition external_potential_methods.F:55
fist_environment_types
Definition fist_environment_types.F:15
fist_force
Definition fist_force.F:16
fist_force::fist_calc_energy_force
subroutine, public fist_calc_energy_force(fist_env, debug)
Calculates the total potential energy, total force, and the total pressure tensor from the potentials...
Definition fist_force.F:112
force_env_methods
Interface for the force calculations.
Definition force_env_methods.F:17
force_env_methods::force_env_calc_energy_force
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Definition force_env_methods.F:218
force_env_methods::force_env_calc_num_pressure
subroutine, public force_env_calc_num_pressure(force_env, dx)
Evaluates the stress tensor and pressure numerically.
Definition force_env_methods.F:720
force_env_methods::force_env_create
subroutine, public force_env_create(force_env, root_section, para_env, globenv, fist_env, qs_env, meta_env, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, force_env_section, mixed_env, embed_env, nnp_env, ipi_env)
creates and initializes a force environment
Definition force_env_methods.F:993
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get_natom
integer function, public force_env_get_natom(force_env)
returns the number of atoms
Definition force_env_types.F:540
force_env_types::use_qmmm
integer, parameter, public use_qmmm
Definition force_env_types.F:83
force_env_types::use_mixed_force
integer, parameter, public use_mixed_force
Definition force_env_types.F:83
force_env_types::use_prog_name
character(len=10), dimension(501:510), parameter, public use_prog_name
Definition force_env_types.F:94
force_env_types::use_eip_force
integer, parameter, public use_eip_force
Definition force_env_types.F:83
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:344
force_env_types::use_embed
integer, parameter, public use_embed
Definition force_env_types.F:83
force_env_types::use_qmmmx
integer, parameter, public use_qmmmx
Definition force_env_types.F:83
force_env_types::force_env_set
subroutine, public force_env_set(force_env, meta_env, fp_env, force_env_section, method_name_id, additional_potential)
changes some attributes of the force_env
Definition force_env_types.F:669
force_env_types::use_qs_force
integer, parameter, public use_qs_force
Definition force_env_types.F:83
force_env_types::use_pwdft_force
integer, parameter, public use_pwdft_force
Definition force_env_types.F:83
force_env_types::use_nnp_force
integer, parameter, public use_nnp_force
Definition force_env_types.F:83
force_env_types::use_ipi
integer, parameter, public use_ipi
Definition force_env_types.F:83
force_env_types::use_fist_force
integer, parameter, public use_fist_force
Definition force_env_types.F:83
force_env_utils
Util force_env module.
Definition force_env_utils.F:12
force_env_utils::write_forces
subroutine, public write_forces(particles, iw, label, ndigits, unit_string, total_force, grand_total_force, zero_force_core_shell_atom)
Write forces either to the screen or to a file.
Definition force_env_utils.F:420
force_env_utils::write_atener
subroutine, public write_atener(iounit, particles, atener, label)
Write the atomic coordinates, types, and energies.
Definition force_env_utils.F:633
force_env_utils::rescale_forces
subroutine, public rescale_forces(force_env)
Rescale forces if requested.
Definition force_env_utils.F:374
fp_methods
methods used in the flexible partitioning scheme
Definition fp_methods.F:14
fp_methods::fp_eval
subroutine, public fp_eval(fp_env, subsys, cell)
Computes the forces and the energy due to the flexible potential & bias, and writes the weights file.
Definition fp_methods.F:55
fparser
This public domain function parser module is intended for applications where a set of mathematical ex...
Definition fparser.F:17
fparser::parsef
subroutine, public parsef(i, funcstr, var)
Parse ith function string FuncStr and compile it into bytecode.
Definition fparser.F:174
fparser::evalf
real(rn) function, public evalf(i, val)
...
Definition fparser.F:206
fparser::evalerrtype
integer, public evalerrtype
Definition fparser.F:33
fparser::evalfd
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
Definition fparser.F:1097
fparser::finalizef
subroutine, public finalizef()
...
Definition fparser.F:127
fparser::initf
subroutine, public initf(n)
...
Definition fparser.F:156
generic_function_input
Creates and parses an input section to define a generic function.
Definition generic_function_input.F:76
generic_function_input::parse_function_section
subroutine, public parse_function_section(gen_section, func_name, xfunction, parameters, values, var_values, i_rep_sec, input_variables, dx, error_limit)
Reads from the input structure all information for generic functions.
Definition generic_function_input.F:262
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition global_types.F:21
global_types::globenv_retain
subroutine, public globenv_retain(globenv)
Retains the global environment globenv.
Definition global_types.F:123
grrm_utils
GRRM interface.
Definition grrm_utils.F:12
grrm_utils::write_grrm
subroutine, public write_grrm(iounit, force_env, particles, energy, dipole, hessian, dipder, polar, fixed_atoms)
Write GRRM interface file.
Definition grrm_utils.F:50
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::use_bazant_eip
integer, parameter, public use_bazant_eip
Definition input_constants.F:775
input_constants::driver_run
integer, parameter, public driver_run
Definition input_constants.F:143
input_constants::mixed_cdft_serial
integer, parameter, public mixed_cdft_serial
Definition input_constants.F:839
input_constants::mix_cdft
integer, parameter, public mix_cdft
Definition input_constants.F:848
input_constants::do_method_ofgpw
integer, parameter, public do_method_ofgpw
Definition input_constants.F:282
input_constants::mix_linear_combination
integer, parameter, public mix_linear_combination
Definition input_constants.F:848
input_constants::do_method_rigpw
integer, parameter, public do_method_rigpw
Definition input_constants.F:282
input_constants::do_method_gpw
integer, parameter, public do_method_gpw
Definition input_constants.F:282
input_constants::dmfet
integer, parameter, public dmfet
Definition input_constants.F:1365
input_constants::use_tersoff_eip
integer, parameter, public use_tersoff_eip
Definition input_constants.F:775
input_constants::use_lenosky_eip
integer, parameter, public use_lenosky_eip
Definition input_constants.F:775
input_constants::ehrenfest
integer, parameter, public ehrenfest
Definition input_constants.F:143
input_constants::use_stillinger_weber_eip
integer, parameter, public use_stillinger_weber_eip
Definition input_constants.F:775
input_constants::mix_coupled
integer, parameter, public mix_coupled
Definition input_constants.F:848
input_constants::debug_run
integer, parameter, public debug_run
Definition input_constants.F:143
input_constants::mix_restrained
integer, parameter, public mix_restrained
Definition input_constants.F:848
input_constants::do_method_gapw
integer, parameter, public do_method_gapw
Definition input_constants.F:282
input_constants::mix_generic
integer, parameter, public mix_generic
Definition input_constants.F:848
input_constants::mol_dyn_run
integer, parameter, public mol_dyn_run
Definition input_constants.F:143
input_constants::dfet
integer, parameter, public dfet
Definition input_constants.F:1365
input_constants::cell_opt_run
integer, parameter, public cell_opt_run
Definition input_constants.F:143
input_constants::do_method_lrigpw
integer, parameter, public do_method_lrigpw
Definition input_constants.F:282
input_constants::mix_minimum
integer, parameter, public mix_minimum
Definition input_constants.F:848
input_constants::do_method_gapw_xc
integer, parameter, public do_method_gapw_xc
Definition input_constants.F:282
input_constants::geo_opt_run
integer, parameter, public geo_opt_run
Definition input_constants.F:143
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_retain
subroutine, public section_vals_retain(section_vals)
retains the given section values (see doc/ReferenceCounting.html)
Definition input_section_types.F:642
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
ipi_environment_types
The environment for the empirical interatomic potential methods.
Definition ipi_environment_types.F:14
ipi_server
i–PI server mode: Communication with i–PI clients
Definition ipi_server.F:14
ipi_server::request_forces
subroutine, public request_forces(ipi_env)
Send atomic positions to a client and retrieve forces.
Definition ipi_server.F:177
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::kpoint_initialize_mos
subroutine, public kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
Initialize a set of MOs and density matrix for each kpoint (kpoint group).
Definition kpoint_methods.F:834
kpoint_methods::kpoint_initialize
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
Definition kpoint_methods.F:379
kpoint_methods::kpoint_env_initialize
subroutine, public kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
Initialize the kpoint environment.
Definition kpoint_methods.F:670
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::set_kpoint_info
subroutine, public set_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Set information in a kpoint environment.
Definition kpoint_types.F:731
kpoint_types::kpoint_reset_initialization
subroutine, public kpoint_reset_initialization(kpoint)
Reset all data derived from a concrete k-point initialization. Input options such as the scheme,...
Definition kpoint_types.F:438
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:613
kpsym
Crystal-symmetry routines originating from the K290/ACMI code.
Definition kpsym.F:14
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_memory
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition machine.F:443
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::abnormal_value
logical function, public abnormal_value(a)
determines if a value is not normal (e.g. for Inf and Nan) based on IO to work also under optimizatio...
Definition mathlib.F:161
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
metadynamics_types
defines types for metadynamics calculation
Definition metadynamics_types.F:13
mixed_cdft_methods
Methods for mixed CDFT calculations.
Definition mixed_cdft_methods.F:14
mixed_cdft_methods::mixed_cdft_calculate_coupling
subroutine, public mixed_cdft_calculate_coupling(force_env)
Driver routine to calculate the electronic coupling(s) between CDFT states.
Definition mixed_cdft_methods.F:1130
mixed_cdft_methods::mixed_cdft_build_weight
subroutine, public mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
Driver routine to handle the build of CDFT weight/gradient in parallel and serial modes.
Definition mixed_cdft_methods.F:420
mixed_cdft_methods::mixed_cdft_init
subroutine, public mixed_cdft_init(force_env, calculate_forces)
Initialize a mixed CDFT calculation.
Definition mixed_cdft_methods.F:154
mixed_energy_types
Definition mixed_energy_types.F:12
mixed_environment_types
Definition mixed_environment_types.F:11
mixed_environment_types::get_mixed_env
subroutine, public get_mixed_env(mixed_env, atomic_kind_set, particle_set, local_particles, local_molecules, molecule_kind_set, molecule_set, cell, cell_ref, mixed_energy, para_env, sub_para_env, subsys, input, results, cdft_control)
Get the MIXED environment.
Definition mixed_environment_types.F:130
mixed_environment_utils
Util mixed_environment.
Definition mixed_environment_utils.F:12
mixed_environment_utils::get_subsys_map_index
subroutine, public get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, map_index, force_eval_embed)
performs mapping of the subsystems of different force_eval
Definition mixed_environment_utils.F:141
mixed_environment_utils::mixed_map_forces
subroutine, public mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, factor, iforce_eval, nforce_eval, map_index, mapping_section, overwrite)
Maps forces between the different force_eval sections/environments.
Definition mixed_environment_utils.F:63
molecule_kind_list_types
represent a simple array based list of the given type
Definition molecule_kind_list_types.F:15
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:518
nnp_environment_types
Data types for neural network potentials.
Definition nnp_environment_types.F:14
nnp_force
Methods dealing with Neural Network potentials.
Definition nnp_force.F:15
nnp_force::nnp_calc_energy_force
subroutine, public nnp_calc_energy_force(nnp, calc_forces)
Calculate the energy and force for a given configuration with the NNP.
Definition nnp_force.F:73
optimize_dmfet_potential
Definition optimize_dmfet_potential.F:8
optimize_dmfet_potential::check_dmfet
subroutine, public check_dmfet(opt_dmfet, qs_env)
...
Definition optimize_dmfet_potential.F:306
optimize_dmfet_potential::release_dmfet_opt
subroutine, public release_dmfet_opt(opt_dmfet)
...
Definition optimize_dmfet_potential.F:188
optimize_dmfet_potential::build_full_dm
subroutine, public build_full_dm(qs_env, dm, open_shell, open_shell_embed, dm_beta)
Builds density matrices from MO coefficients in full matrix format.
Definition optimize_dmfet_potential.F:250
optimize_dmfet_potential::prepare_dmfet_opt
subroutine, public prepare_dmfet_opt(qs_env, opt_dmfet, opt_dmfet_section)
...
Definition optimize_dmfet_potential.F:112
optimize_dmfet_potential::subsys_spin
logical function, public subsys_spin(qs_env)
...
Definition optimize_dmfet_potential.F:93
optimize_embedding_potential
Definition optimize_embedding_potential.F:8
optimize_embedding_potential::print_embed_restart
subroutine, public print_embed_restart(qs_env, dimen_aux, embed_pot_coef, embed_pot, i_iter, embed_pot_spin, open_shell_embed, grid_opt, final_one)
Print embedding potential as a cube and as a binary (for restarting).
Definition optimize_embedding_potential.F:3133
optimize_embedding_potential::make_subsys_embed_pot
subroutine, public make_subsys_embed_pot(qs_env, embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, open_shell_embed, change_spin_sign)
Creates a subsystem embedding potential.
Definition optimize_embedding_potential.F:1032
optimize_embedding_potential::print_rho_spin_diff
subroutine, public print_rho_spin_diff(spin_diff_rho_r, i_iter, qs_env, final_one)
Prints a cube for the (spin_rho_A + spin_rho_B - spin_rho_ref) to be minimized in embedding.
Definition optimize_embedding_potential.F:3078
optimize_embedding_potential::get_max_subsys_diff
subroutine, public get_max_subsys_diff(opt_embed, force_env, subsys_num)
...
Definition optimize_embedding_potential.F:3460
optimize_embedding_potential::print_emb_opt_info
subroutine, public print_emb_opt_info(output_unit, step_num, opt_embed)
...
Definition optimize_embedding_potential.F:3397
optimize_embedding_potential::init_embed_pot
subroutine, public init_embed_pot(qs_env, embed_pot, add_const_pot, fermi_amaldi, const_pot, open_shell_embed, spin_embed_pot, pot_diff, coulomb_guess, grid_opt)
...
Definition optimize_embedding_potential.F:208
optimize_embedding_potential::understand_spin_states
subroutine, public understand_spin_states(force_env, ref_subsys_number, change_spin, open_shell_embed, all_nspins)
Find out whether we need to swap alpha- and beta- spind densities in the second subsystem.
Definition optimize_embedding_potential.F:129
optimize_embedding_potential::read_embed_pot
subroutine, public read_embed_pot(qs_env, embed_pot, spin_embed_pot, section, opt_embed)
...
Definition optimize_embedding_potential.F:665
optimize_embedding_potential::get_prev_density
subroutine, public get_prev_density(opt_embed, force_env, subsys_num)
...
Definition optimize_embedding_potential.F:3430
optimize_embedding_potential::coulomb_guess
subroutine, public coulomb_guess(v_rspace, rhs, mapping_section, qs_env, nforce_eval, iforce_eval, eta)
Calculates subsystem Coulomb potential from the RESP charges of the total system.
Definition optimize_embedding_potential.F:941
optimize_embedding_potential::print_rho_diff
subroutine, public print_rho_diff(diff_rho_r, i_iter, qs_env, final_one)
Prints a cube for the (rho_A + rho_B - rho_ref) to be minimized in embedding.
Definition optimize_embedding_potential.F:3027
optimize_embedding_potential::conv_check_embed
subroutine, public conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit)
...
Definition optimize_embedding_potential.F:3493
optimize_embedding_potential::step_control
subroutine, public step_control(opt_embed)
Controls the step, changes the trust radius if needed in maximization of the V_emb.
Definition optimize_embedding_potential.F:2182
optimize_embedding_potential::opt_embed_step
subroutine, public opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, qs_env)
Takes maximization step in embedding potential optimization.
Definition optimize_embedding_potential.F:1475
optimize_embedding_potential::calculate_embed_pot_grad
subroutine, public calculate_embed_pot_grad(qs_env, diff_rho_r, diff_rho_spin, opt_embed)
Calculates the derivative of the embedding potential wrt to the expansion coefficients.
Definition optimize_embedding_potential.F:1080
optimize_embedding_potential::print_pot_simple_grid
subroutine, public print_pot_simple_grid(qs_env, embed_pot, embed_pot_spin, i_iter, open_shell_embed, final_one, qs_env_cluster)
Prints a volumetric file: X Y Z value for interfacing with external programs.
Definition optimize_embedding_potential.F:3231
optimize_embedding_potential::prepare_embed_opt
subroutine, public prepare_embed_opt(qs_env, opt_embed, opt_embed_section)
Creates and allocates objects for optimization of embedding potential.
Definition optimize_embedding_potential.F:293
optimize_embedding_potential::release_opt_embed
subroutine, public release_opt_embed(opt_embed)
Deallocate stuff for optimizing embedding potential.
Definition optimize_embedding_potential.F:864
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::debye
real(kind=dp), parameter, public debye
Definition physcon.F:201
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
pwdft_environment_types
The type definitions for the PWDFT environment.
Definition pwdft_environment_types.F:14
pwdft_environment
Methods and functions on the PWDFT environment.
Definition pwdft_environment.F:14
pwdft_environment::pwdft_calc_energy_force
subroutine, public pwdft_calc_energy_force(pwdft_env, calculate_forces, calculate_stress)
Calculate energy and forces within the PWDFT/SIRIUS code.
Definition pwdft_environment.F:210
qmmm_force
Calculates QM/MM energy and forces.
Definition qmmm_force.F:14
qmmm_force::qmmm_calc_energy_force
subroutine, public qmmm_calc_energy_force(qmmm_env, calc_force, consistent_energies, linres)
calculates the qm/mm energy and forces
Definition qmmm_force.F:76
qmmm_types
Basic container type for QM/MM.
Definition qmmm_types.F:12
qmmm_util
Definition qmmm_util.F:13
qmmm_util::apply_qmmm_translate
subroutine, public apply_qmmm_translate(qmmm_env)
Apply translation to the full system in order to center the QM system into the QM box.
Definition qmmm_util.F:375
qmmmx_force
Calculates QM/MM energy and forces with Force-Mixing.
Definition qmmmx_force.F:14
qmmmx_force::qmmmx_calc_energy_force
subroutine, public qmmmx_calc_energy_force(qmmmx_env, calc_force, consistent_energies, linres, require_consistent_energy_force)
calculates the qm/mm energy and forces
Definition qmmmx_force.F:66
qmmmx_types
Basic container type for QM/MM with force mixing.
Definition qmmmx_types.F:12
qs_apt_fdiff_methods
Atomic Polarization Tensor calculation by dF/d(E-field) finite differences.
Definition qs_apt_fdiff_methods.F:13
qs_apt_fdiff_methods::apt_fdiff
subroutine, public apt_fdiff(force_env)
Calculate Atomic Polarization Tensors by dF/d(E-field) finite differences.
Definition qs_apt_fdiff_methods.F:103
qs_basis_rotation_methods
Definition qs_basis_rotation_methods.F:9
qs_basis_rotation_methods::qs_basis_rotation
subroutine, public qs_basis_rotation(qs_env, kpoints, basis_type, storage_type)
Construct basis set rotation matrices.
Definition qs_basis_rotation_methods.F:54
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1134
qs_force
Quickstep force driver routine.
Definition qs_force.F:12
qs_force::qs_calc_energy_force
subroutine, public qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres)
...
Definition qs_force.F:113
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_wf_history_types
interpolate the wavefunctions to speed up the convergence when doing MD
Definition qs_wf_history_types.F:16
qs_wf_history_types::wfi_clear
subroutine, public wfi_clear(wf_history)
Clear stored wavefunction snapshots while preserving history settings.
Definition qs_wf_history_types.F:227
restraint
Handles all possible kinds of restraints in CP2K.
Definition restraint.F:14
restraint::restraint_control
subroutine, public restraint_control(force_env)
Computes restraints.
Definition restraint.F:67
scine_utils
SCINE interface.
Definition scine_utils.F:12
scine_utils::write_scine
subroutine, public write_scine(iounit, force_env, particles, energy, hessian)
Write SCINE interface file.
Definition scine_utils.F:46
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::compress
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Definition string_utilities.F:3191
virial_methods
Definition virial_methods.F:12
virial_methods::write_stress_tensor
subroutine, public write_stress_tensor(pv_virial, iw, cell, unit_string, numerical)
Print stress tensor to output file.
Definition virial_methods.F:261
virial_methods::write_stress_tensor_components
subroutine, public write_stress_tensor_components(virial, iw, cell, unit_string)
...
Definition virial_methods.F:167
virial_types
Definition virial_types.F:13
virial_types::zero_virial
subroutine, public zero_virial(virial, reset)
...
Definition virial_types.F:174
virial_types::symmetrize_virial
subroutine, public symmetrize_virial(virial)
Symmetrize the virial components.
Definition virial_types.F:118
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_result_types::cp_result_p_type
Definition cp_result_types.F:81
cp_result_types::cp_result_type
contains arbitrary information which need to be stored
Definition cp_result_types.F:73
cp_subsys_types::cp_subsys_p_type
represent a pointer to a subsys, to be able to create arrays of pointers
Definition cp_subsys_types.F:123
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
eip_environment_types::eip_environment_type
The empirical interatomic potential environment.
Definition eip_environment_types.F:83
embed_types::embed_env_type
Embedding environment type.
Definition embed_types.F:127
embed_types::opt_dmfet_pot_type
Type containing main data for matrix embedding potential optimization.
Definition embed_types.F:105
embed_types::opt_embed_pot_type
Type containing main data for embedding potential optimization.
Definition embed_types.F:67
fist_environment_types::fist_environment_type
Definition fist_environment_types.F:71
force_env_types::force_env_p_type
allows for the creation of an array of force_env
Definition force_env_types.F:174
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
ipi_environment_types::ipi_environment_type
The i–PI environment.
Definition ipi_environment_types.F:71
kpoint_types::kpoint_sym_type
Keeps symmetry information about a specific k-point.
Definition kpoint_types.F:144
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:196
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
metadynamics_types::meta_env_type
defines meta_env type
Definition metadynamics_types.F:95
mixed_energy_types::mixed_energy_type
Definition mixed_energy_types.F:22
mixed_energy_types::mixed_force_type
Definition mixed_energy_types.F:28
mixed_environment_types::mixed_environment_type
Definition mixed_environment_types.F:64
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:153
nnp_environment_types::nnp_type
Main data type collecting all relevant data for neural network potentials.
Definition nnp_environment_types.F:97
particle_list_types::particle_list_p_type
represents a pointer to a list
Definition particle_list_types.F:59
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
pwdft_environment_types::pwdft_environment_type
The PWDFT environment type.
Definition pwdft_environment_types.F:67
qmmm_types::qmmm_env_type
Definition qmmm_types.F:34
qmmmx_types::qmmmx_env_type
Definition qmmmx_types.F:27
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_wf_history_types::qs_wf_history_type
keeps track of the previous wavefunctions and can extrapolate them for the next step of md
Definition qs_wf_history_types.F:105
virial_types::virial_p_type
Definition virial_types.F:52
virial_types::virial_type
Definition virial_types.F:26