(git:871dbd5)
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 force_fields_util, ONLY: get_generic_info
81 USE fp_methods, ONLY: fp_eval
82 USE fparser, ONLY: evalerrtype,&
83 evalf,&
84 evalfd,&
85 finalizef,&
86 initf,&
87 parsef
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("")
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)) &
408 cpabort("Potential energy is an abnormal value (NaN/Inf).")
409
410 ! Print forces, if requested
411 print_forces = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%FORCES", &
412 extension=".xyz")
413 IF ((print_forces > 0) .AND. calculate_forces) THEN
414 CALL force_env_get(force_env, subsys=subsys)
415 CALL cp_subsys_get(subsys, &
416 core_particles=core_particles, &
417 particles=particles, &
418 shell_particles=shell_particles)
419 ! Variable precision output of the forces
420 CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%NDIGITS", &
421 i_val=ndigits)
422 CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%FORCE_UNIT", &
423 c_val=unit_string)
424 IF (ASSOCIATED(core_particles) .OR. ASSOCIATED(shell_particles)) THEN
425 CALL write_forces(particles, print_forces, "Atomic", ndigits, unit_string, &
426 total_force, zero_force_core_shell_atom=.true.)
427 grand_total_force(1:3) = total_force(1:3)
428 IF (ASSOCIATED(core_particles)) THEN
429 CALL write_forces(core_particles, print_forces, "Core particle", ndigits, &
430 unit_string, total_force, zero_force_core_shell_atom=.false.)
431 grand_total_force(:) = grand_total_force(:) + total_force(:)
432 END IF
433 IF (ASSOCIATED(shell_particles)) THEN
434 CALL write_forces(shell_particles, print_forces, "Shell particle", ndigits, &
435 unit_string, total_force, zero_force_core_shell_atom=.false., &
436 grand_total_force=grand_total_force)
437 END IF
438 ELSE
439 CALL write_forces(particles, print_forces, "Atomic", ndigits, unit_string, total_force)
440 END IF
441 END IF
442 CALL cp_print_key_finished_output(print_forces, logger, force_env%force_env_section, "PRINT%FORCES")
443
444 ! Write stress tensor
445 IF (virial%pv_availability) THEN
446 ! If the virial is defined but we are not computing forces let's zero the
447 ! virial for consistency
448 IF (calculate_forces .AND. calculate_stress_tensor) THEN
449 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
450 extension=".stress_tensor")
451 IF (output_unit > 0) THEN
452 CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%COMPONENTS", &
453 l_val=print_components)
454 CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%STRESS_UNIT", &
455 c_val=unit_string)
456 IF (print_components) THEN
457 IF ((.NOT. virial%pv_numer) .AND. (force_env%in_use == use_qs_force)) THEN
458 CALL write_stress_tensor_components(virial, output_unit, cell, unit_string)
459 END IF
460 END IF
461 CALL write_stress_tensor(virial%pv_virial, output_unit, cell, unit_string, virial%pv_numer)
462 END IF
463 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
464 "PRINT%STRESS_TENSOR")
465 ELSE
466 CALL zero_virial(virial, reset=.false.)
467 END IF
468 ELSE
469 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
470 extension=".stress_tensor")
471 IF (output_unit > 0) THEN
472 CALL cp_warn(__location__, "To print the stress tensor switch on the "// &
473 "virial evaluation with the keyword: STRESS_TENSOR")
474 END IF
475 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
476 "PRINT%STRESS_TENSOR")
477 END IF
478
479 ! Atomic energy
480 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
481 extension=".Log")
482 IF (atprop_env%energy) THEN
483 CALL force_env%para_env%sum(atprop_env%atener)
484 CALL force_env_get(force_env, potential_energy=e_pot)
485 IF (output_unit > 0) THEN
486 IF (logger%iter_info%print_level >= low_print_level) THEN
487 CALL cp_subsys_get(subsys=subsys, particles=particles)
488 CALL write_atener(output_unit, particles, atprop_env%atener, "Mulliken Atomic Energies")
489 END IF
490 sum_energy = accurate_sum(atprop_env%atener(:))
491 checksum = abs(e_pot - sum_energy)
492 WRITE (unit=output_unit, fmt="(/,(T2,A,T56,F25.13))") &
493 "Potential energy (Atomic):", sum_energy, &
494 "Potential energy (Total) :", e_pot, &
495 "Difference :", checksum
496 cpassert((checksum < ateps*abs(e_pot)))
497 END IF
498 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
499 "PRINT%PROGRAM_RUN_INFO")
500 END IF
501
502 ! Print GRMM interface file
503 print_grrm = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%GRRM", &
504 file_position="REWIND", extension=".rrm")
505 IF (print_grrm > 0) THEN
506 CALL force_env_get(force_env, subsys=subsys)
507 CALL cp_subsys_get(subsys=subsys, particles=particles, &
508 molecule_kinds=molecule_kinds)
509 ! Count the number of fixed atoms
510 nfixed_atoms_total = 0
511 nkind = molecule_kinds%n_els
512 molecule_kind_set => molecule_kinds%els
513 DO ikind = 1, nkind
514 molecule_kind => molecule_kind_set(ikind)
515 CALL get_molecule_kind(molecule_kind, nfixd=nfixed_atoms)
516 nfixed_atoms_total = nfixed_atoms_total + nfixed_atoms
517 END DO
518 !
519 CALL write_grrm(print_grrm, force_env, particles%els, e_pot, fixed_atoms=nfixed_atoms_total)
520 END IF
521 CALL cp_print_key_finished_output(print_grrm, logger, force_env%force_env_section, "PRINT%GRRM")
522
523 print_scine = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%SCINE", &
524 file_position="REWIND", extension=".scine")
525 IF (print_scine > 0) THEN
526 CALL force_env_get(force_env, subsys=subsys)
527 CALL cp_subsys_get(subsys=subsys, particles=particles)
528 !
529 CALL write_scine(print_scine, force_env, particles%els, e_pot)
530 END IF
531 CALL cp_print_key_finished_output(print_scine, logger, force_env%force_env_section, "PRINT%SCINE")
532
533 END SUBROUTINE force_env_calc_energy_force
534
535! **************************************************************************************************
536!> \brief Rebuild k-point data for geometries whose atomic symmetry can change.
537!> Atomic k-point symmetry may change when atoms or cell vectors move.
538!> \param force_env ...
539!> \param fd_energy ...
540! **************************************************************************************************
541 SUBROUTINE force_env_refresh_kpoint_symmetry(force_env, fd_energy)
542
543 TYPE(force_env_type), POINTER :: force_env
544 LOGICAL, INTENT(IN) :: fd_energy
545
546 CHARACTER(LEN=default_string_length) :: kp_scheme
547 INTEGER :: run_type_id
548 LOGICAL :: debug_full_kpoint_symmetry, debug_inversion_only, do_kpoints, dynamic_symmetry, &
549 full_grid, input_full_grid, input_inversion_symmetry_only, inversion_symmetry_only, &
550 kpoint_symmetry, moving_geometry, use_full_grid, use_inversion_symmetry_only
551 TYPE(cell_type), POINTER :: cell
552 TYPE(cp_blacs_env_type), POINTER :: blacs_env
553 TYPE(dft_control_type), POINTER :: dft_control
554 TYPE(global_environment_type), POINTER :: globenv
555 TYPE(kpoint_type), POINTER :: kpoints
556 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
557 TYPE(mp_para_env_type), POINTER :: para_env
558 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
559 TYPE(qs_wf_history_type), POINTER :: wf_history
560 TYPE(section_vals_type), POINTER :: input, kpoint_section
561
562 IF (.NOT. ASSOCIATED(force_env)) RETURN
563 IF (force_env%in_use /= use_qs_force) RETURN
564
565 NULLIFY (globenv)
566 CALL force_env_get(force_env, globenv=globenv)
567 IF (.NOT. ASSOCIATED(globenv)) RETURN
568 run_type_id = globenv%run_type_id
569 moving_geometry = .false.
570 SELECT CASE (run_type_id)
571 CASE (cell_opt_run, driver_run, ehrenfest, geo_opt_run, mol_dyn_run)
572 moving_geometry = .true.
573 CASE DEFAULT
574 moving_geometry = .false.
575 END SELECT
576 IF (run_type_id /= debug_run .AND. .NOT. moving_geometry) RETURN
577
578 NULLIFY (blacs_env, cell, dft_control, input, kpoint_section, kpoints, mos, para_env, &
579 particle_set, wf_history)
580 CALL get_qs_env(qs_env=force_env%qs_env, &
581 blacs_env=blacs_env, &
582 cell=cell, &
583 dft_control=dft_control, &
584 do_kpoints=do_kpoints, &
585 input=input, &
586 kpoints=kpoints, &
587 mos=mos, &
588 para_env=para_env, &
589 particle_set=particle_set, &
590 wf_history=wf_history)
591 IF (.NOT. do_kpoints) RETURN
592
593 CALL get_kpoint_info(kpoints, kp_scheme=kp_scheme, symmetry=kpoint_symmetry, full_grid=full_grid, &
594 inversion_symmetry_only=inversion_symmetry_only)
595 IF (.NOT. kpoint_symmetry) RETURN
596 IF (trim(kp_scheme) /= "MONKHORST-PACK" .AND. trim(kp_scheme) /= "MACDONALD" .AND. &
597 trim(kp_scheme) /= "GENERAL") RETURN
598
599 input_full_grid = full_grid
600 input_inversion_symmetry_only = inversion_symmetry_only
601 debug_full_kpoint_symmetry = .false.
602 IF (ASSOCIATED(input)) THEN
603 kpoint_section => section_vals_get_subs_vals(input, "DFT%KPOINTS")
604 CALL section_vals_val_get(kpoint_section, "FULL_GRID", l_val=input_full_grid)
605 CALL section_vals_val_get(kpoint_section, "INVERSION_SYMMETRY_ONLY", &
606 l_val=input_inversion_symmetry_only)
607 CALL section_vals_val_get(kpoint_section, "DEBUG_FULL_KPOINT_SYMMETRY", &
608 l_val=debug_full_kpoint_symmetry)
609 END IF
610 ! Moving geometries and DEBUG finite differences must not reuse atomic symmetry from
611 ! another geometry. Rebuild the k-point symmetry from the current cell and positions.
612 ! Keep numerical finite-difference energies and DFTB DEBUG checks on the established
613 ! inversion/time-reversal reduction.
614 debug_inversion_only = run_type_id == debug_run .AND. .NOT. debug_full_kpoint_symmetry .AND. &
615 (fd_energy .OR. dft_control%qs_control%dftb)
616 use_full_grid = input_full_grid
617 use_inversion_symmetry_only = (input_inversion_symmetry_only .OR. debug_inversion_only) .AND. &
618 (.NOT. use_full_grid)
619 dynamic_symmetry = kpoint_symmetry .AND. .NOT. use_full_grid .AND. &
620 .NOT. use_inversion_symmetry_only
621 IF (run_type_id == debug_run .AND. .NOT. fd_energy .AND. .NOT. dynamic_symmetry .AND. &
622 (full_grid .EQV. use_full_grid) .AND. &
623 (inversion_symmetry_only .EQV. use_inversion_symmetry_only)) THEN
624 CALL qs_basis_rotation(force_env%qs_env, kpoints)
625 RETURN
626 END IF
627 IF (moving_geometry .AND. .NOT. dynamic_symmetry) RETURN
628 IF (moving_geometry .AND. .NOT. kpoint_has_nontrivial_atomic_symmetry(kpoints)) RETURN
629 CALL set_kpoint_info(kpoints, full_grid=use_full_grid, &
630 inversion_symmetry_only=use_inversion_symmetry_only)
631
632 CALL kpoint_reset_initialization(kpoints)
633 CALL kpoint_initialize(kpoints, particle_set, cell)
634 CALL kpoint_env_initialize(kpoints, para_env, blacs_env, with_aux_fit=dft_control%do_admm)
635 CALL kpoint_initialize_mos(kpoints, mos)
636 CALL wfi_clear(wf_history)
637 CALL qs_basis_rotation(force_env%qs_env, kpoints)
638
639 END SUBROUTINE force_env_refresh_kpoint_symmetry
640
641! **************************************************************************************************
642!> \brief Return whether the current reduced mesh uses nontrivial atomic symmetry operations.
643!> \param kpoints ...
644!> \return has_symmetry
645! **************************************************************************************************
646 FUNCTION kpoint_has_nontrivial_atomic_symmetry(kpoints) RESULT(has_symmetry)
647
648 TYPE(kpoint_type), POINTER :: kpoints
649 LOGICAL :: has_symmetry
650
651 INTEGER :: iatom, ik, isym, natom
652 REAL(kind=dp), DIMENSION(3, 3) :: eye3
653 TYPE(kpoint_sym_type), POINTER :: kpsym
654
655 has_symmetry = .false.
656 IF (.NOT. ASSOCIATED(kpoints)) RETURN
657 IF (.NOT. ASSOCIATED(kpoints%kp_sym)) RETURN
658
659 eye3 = 0.0_dp
660 eye3(1, 1) = 1.0_dp
661 eye3(2, 2) = 1.0_dp
662 eye3(3, 3) = 1.0_dp
663
664 DO ik = 1, kpoints%nkp
665 kpsym => kpoints%kp_sym(ik)%kpoint_sym
666 IF (.NOT. ASSOCIATED(kpsym)) cycle
667 IF (.NOT. kpsym%apply_symmetry) cycle
668 IF (.NOT. ASSOCIATED(kpsym%rot)) cycle
669 IF (.NOT. ASSOCIATED(kpsym%f0)) cycle
670 IF (.NOT. ASSOCIATED(kpsym%fcell)) cycle
671
672 natom = SIZE(kpsym%f0, 1)
673 DO isym = 1, SIZE(kpsym%rot, 3)
674 IF (maxval(abs(kpsym%rot(1:3, 1:3, isym) - eye3(1:3, 1:3))) > 1.e-12_dp .OR. &
675 any(kpsym%fcell(1:3, 1:natom, isym) /= 0)) THEN
676 has_symmetry = .true.
677 RETURN
678 END IF
679 DO iatom = 1, natom
680 IF (kpsym%f0(iatom, isym) /= iatom) THEN
681 has_symmetry = .true.
682 RETURN
683 END IF
684 END DO
685 END DO
686 END DO
687
688 END FUNCTION kpoint_has_nontrivial_atomic_symmetry
689
690! **************************************************************************************************
691!> \brief Evaluates the stress tensor and pressure numerically
692!> \param force_env ...
693!> \param dx ...
694!> \par History
695!> 10.2005 created [JCS]
696!> 05.2009 Teodoro Laino [tlaino] - rewriting for general force_env
697!>
698!> \author JCS
699! **************************************************************************************************
700 SUBROUTINE force_env_calc_num_pressure(force_env, dx)
701
702 TYPE(force_env_type), POINTER :: force_env
703 REAL(kind=dp), INTENT(IN), OPTIONAL :: dx
704
705 REAL(kind=dp), PARAMETER :: default_dx = 0.001_dp
706
707 CHARACTER(LEN=default_string_length) :: unit_string
708 INTEGER :: i, ip, iq, j, k, method_id, natom, &
709 ncore, nshell, output_unit, symmetry_id
710 LOGICAL :: use_sym_strain_2d
711 REAL(kind=dp) :: dx_w, eps_w
712 REAL(kind=dp), DIMENSION(2) :: numer_energy
713 REAL(kind=dp), DIMENSION(3) :: s
714 REAL(kind=dp), DIMENSION(3, 3) :: hmat_deformed, numer_pv_2d, &
715 numer_stress, strain
716 REAL(kind=dp), DIMENSION(:, :), POINTER :: ref_pos_atom, ref_pos_core, ref_pos_shell
717 TYPE(cell_type), POINTER :: cell, cell_local
718 TYPE(cp_logger_type), POINTER :: logger
719 TYPE(cp_subsys_type), POINTER :: subsys
720 TYPE(dft_control_type), POINTER :: dft_control
721 TYPE(global_environment_type), POINTER :: globenv
722 TYPE(particle_list_type), POINTER :: core_particles, particles, &
723 shell_particles
724 TYPE(virial_type), POINTER :: virial
725
726 NULLIFY (cell_local)
727 NULLIFY (dft_control)
728 NULLIFY (core_particles)
729 NULLIFY (particles)
730 NULLIFY (shell_particles)
731 NULLIFY (ref_pos_atom)
732 NULLIFY (ref_pos_core)
733 NULLIFY (ref_pos_shell)
734 natom = 0
735 method_id = 0
736 ncore = 0
737 nshell = 0
738 numer_pv_2d = 0.0_dp
739 numer_stress = 0.0_dp
740 use_sym_strain_2d = .false.
741
742 logger => cp_get_default_logger()
743
744 dx_w = default_dx
745 IF (PRESENT(dx)) dx_w = dx
746 CALL force_env_get(force_env, subsys=subsys, globenv=globenv, in_use=method_id)
747 CALL cp_subsys_get(subsys, &
748 core_particles=core_particles, &
749 particles=particles, &
750 shell_particles=shell_particles, &
751 virial=virial)
752 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
753 extension=".stress_tensor")
754 IF (output_unit > 0) THEN
755 WRITE (output_unit, "(/A,A/)") " **************************** ", &
756 "NUMERICAL STRESS ********************************"
757 END IF
758
759 ! Save all original particle positions
760 natom = particles%n_els
761 ALLOCATE (ref_pos_atom(natom, 3))
762 DO i = 1, natom
763 ref_pos_atom(i, :) = particles%els(i)%r
764 END DO
765 IF (ASSOCIATED(core_particles)) THEN
766 ncore = core_particles%n_els
767 ALLOCATE (ref_pos_core(ncore, 3))
768 DO i = 1, ncore
769 ref_pos_core(i, :) = core_particles%els(i)%r
770 END DO
771 END IF
772 IF (ASSOCIATED(shell_particles)) THEN
773 nshell = shell_particles%n_els
774 ALLOCATE (ref_pos_shell(nshell, 3))
775 DO i = 1, nshell
776 ref_pos_shell(i, :) = shell_particles%els(i)%r
777 END DO
778 END IF
779 CALL force_env_get(force_env, cell=cell)
780 ! Save cell symmetry (distorted cell has no symmetry)
781 symmetry_id = cell%symmetry_id
782 cell%symmetry_id = cell_sym_triclinic
783 !
784 CALL cell_create(cell_local)
785 CALL cell_clone(cell, cell_local)
786 IF (count(cell_local%perd /= 0) == 2 .AND. method_id == use_qs_force) THEN
787 CALL get_qs_env(qs_env=force_env%qs_env, dft_control=dft_control)
788 SELECT CASE (dft_control%qs_control%method_id)
789 CASE (do_method_gapw, do_method_gapw_xc, do_method_gpw, &
790 do_method_lrigpw, do_method_ofgpw, do_method_rigpw)
791 use_sym_strain_2d = .true.
792 END SELECT
793 END IF
794 ! First change box
795 DO ip = 1, 3
796 DO iq = 1, 3
797 IF (use_sym_strain_2d) THEN
798 IF (cell_local%perd(ip) == 0 .OR. cell_local%perd(iq) == 0) cycle
799 IF (iq < ip) cycle
800 END IF
801 IF (virial%pv_diagonal .AND. (ip /= iq)) cycle
802 DO k = 1, 2
803 hmat_deformed = cell_local%hmat
804 IF (use_sym_strain_2d) THEN
805 eps_w = -(-1.0_dp)**k*dx_w
806 strain = 0.0_dp
807 DO i = 1, 3
808 strain(i, i) = 1.0_dp
809 END DO
810 IF (ip == iq) THEN
811 strain(ip, ip) = strain(ip, ip) + eps_w
812 ELSE
813 strain(ip, iq) = strain(ip, iq) + 0.5_dp*eps_w
814 strain(iq, ip) = strain(iq, ip) + 0.5_dp*eps_w
815 END IF
816 hmat_deformed = matmul(strain, cell_local%hmat)
817 ELSE
818 hmat_deformed(ip, iq) = hmat_deformed(ip, iq) - (-1.0_dp)**k*dx_w
819 END IF
820 cell%hmat = hmat_deformed
821 CALL init_cell(cell)
822 ! Scale positions
823 DO i = 1, natom
824 CALL real_to_scaled(s, ref_pos_atom(i, 1:3), cell_local)
825 CALL scaled_to_real(particles%els(i)%r, s, cell)
826 END DO
827 DO i = 1, ncore
828 CALL real_to_scaled(s, ref_pos_core(i, 1:3), cell_local)
829 CALL scaled_to_real(core_particles%els(i)%r, s, cell)
830 END DO
831 DO i = 1, nshell
832 CALL real_to_scaled(s, ref_pos_shell(i, 1:3), cell_local)
833 CALL scaled_to_real(shell_particles%els(i)%r, s, cell)
834 END DO
835 ! Compute energies
836 CALL force_env_calc_energy_force(force_env, &
837 calc_force=.false., &
838 consistent_energies=.true., &
839 calc_stress_tensor=.false.)
840 CALL force_env_get(force_env, potential_energy=numer_energy(k))
841 ! Reset cell
842 cell%hmat = cell_local%hmat
843 END DO
844 CALL init_cell(cell)
845 IF (use_sym_strain_2d) THEN
846 numer_pv_2d(ip, iq) = -0.5_dp*(numer_energy(1) - numer_energy(2))/dx_w
847 numer_pv_2d(iq, ip) = numer_pv_2d(ip, iq)
848 IF (output_unit > 0) THEN
849 IF (globenv%run_type_id == debug_run) THEN
850 WRITE (unit=output_unit, fmt="(/,T2,A,T19,A,F7.4,A,T44,A,F7.4,A,T69,A)") &
851 "DEBUG|", "E(e"//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
852 "E(e"//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
853 "pv(numerical)"
854 WRITE (unit=output_unit, fmt="(T2,A,2(1X,F24.8),1X,F22.8)") &
855 "DEBUG|", numer_energy(1:2), numer_pv_2d(ip, iq)
856 ELSE
857 WRITE (unit=output_unit, fmt="(/,T7,A,F7.4,A,T27,A,F7.4,A,T49,A)") &
858 "E(e"//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
859 "E(e"//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
860 "pv(numerical)"
861 WRITE (unit=output_unit, fmt="(3(1X,F19.8))") &
862 numer_energy(1:2), numer_pv_2d(ip, iq)
863 END IF
864 END IF
865 ELSE
866 numer_stress(ip, iq) = 0.5_dp*(numer_energy(1) - numer_energy(2))/dx_w
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("//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
871 "E("//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
872 "f(numerical)"
873 WRITE (unit=output_unit, fmt="(T2,A,2(1X,F24.8),1X,F22.8)") &
874 "DEBUG|", numer_energy(1:2), numer_stress(ip, iq)
875 ELSE
876 WRITE (unit=output_unit, fmt="(/,T7,A,F7.4,A,T27,A,F7.4,A,T49,A)") &
877 "E("//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
878 "E("//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
879 "f(numerical)"
880 WRITE (unit=output_unit, fmt="(3(1X,F19.8))") &
881 numer_energy(1:2), numer_stress(ip, iq)
882 END IF
883 END IF
884 END IF
885 END DO
886 END DO
887
888 ! Reset positions and rebuild original environment
889 cell%symmetry_id = symmetry_id
890 CALL init_cell(cell)
891 DO i = 1, natom
892 particles%els(i)%r = ref_pos_atom(i, :)
893 END DO
894 DO i = 1, ncore
895 core_particles%els(i)%r = ref_pos_core(i, :)
896 END DO
897 DO i = 1, nshell
898 shell_particles%els(i)%r = ref_pos_shell(i, :)
899 END DO
900 CALL force_env_calc_energy_force(force_env, &
901 calc_force=.false., &
902 consistent_energies=.true., &
903 calc_stress_tensor=.false.)
904
905 ! Computing pv_test
906 virial%pv_virial = 0.0_dp
907 IF (use_sym_strain_2d) THEN
908 virial%pv_virial = numer_pv_2d
909 ELSE
910 DO i = 1, 3
911 DO j = 1, 3
912 DO k = 1, 3
913 virial%pv_virial(i, j) = virial%pv_virial(i, j) - &
914 0.5_dp*(numer_stress(i, k)*cell_local%hmat(j, k) + &
915 numer_stress(j, k)*cell_local%hmat(i, k))
916 END DO
917 END DO
918 END DO
919 END IF
920 IF (output_unit > 0) THEN
921 IF (globenv%run_type_id == debug_run) THEN
922 CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%FORCE_UNIT", &
923 c_val=unit_string)
924 CALL write_stress_tensor(virial%pv_virial, output_unit, cell, unit_string, virial%pv_numer)
925 END IF
926 WRITE (output_unit, "(/,A,/)") " **************************** "// &
927 "NUMERICAL STRESS END *****************************"
928 END IF
929
930 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
931 "PRINT%STRESS_TENSOR")
932
933 ! Release storage
934 IF (ASSOCIATED(ref_pos_atom)) THEN
935 DEALLOCATE (ref_pos_atom)
936 END IF
937 IF (ASSOCIATED(ref_pos_core)) THEN
938 DEALLOCATE (ref_pos_core)
939 END IF
940 IF (ASSOCIATED(ref_pos_shell)) THEN
941 DEALLOCATE (ref_pos_shell)
942 END IF
943 IF (ASSOCIATED(cell_local)) CALL cell_release(cell_local)
944
945 END SUBROUTINE force_env_calc_num_pressure
946
947! **************************************************************************************************
948!> \brief creates and initializes a force environment
949!> \param force_env the force env to create
950!> \param root_section ...
951!> \param para_env ...
952!> \param globenv ...
953!> \param fist_env , qs_env: exactly one of these should be
954!> associated, the one that is active
955!> \param qs_env ...
956!> \param meta_env ...
957!> \param sub_force_env ...
958!> \param qmmm_env ...
959!> \param qmmmx_env ...
960!> \param eip_env ...
961!> \param pwdft_env ...
962!> \param force_env_section ...
963!> \param mixed_env ...
964!> \param embed_env ...
965!> \param nnp_env ...
966!> \param ipi_env ...
967!> \par History
968!> 04.2003 created [fawzi]
969!> \author fawzi
970! **************************************************************************************************
971 SUBROUTINE force_env_create(force_env, root_section, para_env, globenv, fist_env, &
972 qs_env, meta_env, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, force_env_section, &
973 mixed_env, embed_env, nnp_env, ipi_env)
974
975 TYPE(force_env_type), POINTER :: force_env
976 TYPE(section_vals_type), POINTER :: root_section
977 TYPE(mp_para_env_type), POINTER :: para_env
978 TYPE(global_environment_type), POINTER :: globenv
979 TYPE(fist_environment_type), OPTIONAL, POINTER :: fist_env
980 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
981 TYPE(meta_env_type), OPTIONAL, POINTER :: meta_env
982 TYPE(force_env_p_type), DIMENSION(:), OPTIONAL, &
983 POINTER :: sub_force_env
984 TYPE(qmmm_env_type), OPTIONAL, POINTER :: qmmm_env
985 TYPE(qmmmx_env_type), OPTIONAL, POINTER :: qmmmx_env
986 TYPE(eip_environment_type), OPTIONAL, POINTER :: eip_env
987 TYPE(pwdft_environment_type), OPTIONAL, POINTER :: pwdft_env
988 TYPE(section_vals_type), POINTER :: force_env_section
989 TYPE(mixed_environment_type), OPTIONAL, POINTER :: mixed_env
990 TYPE(embed_env_type), OPTIONAL, POINTER :: embed_env
991 TYPE(nnp_type), OPTIONAL, POINTER :: nnp_env
992 TYPE(ipi_environment_type), OPTIONAL, POINTER :: ipi_env
993
994 ALLOCATE (force_env)
995 NULLIFY (force_env%fist_env, force_env%qs_env, &
996 force_env%para_env, force_env%globenv, &
997 force_env%meta_env, force_env%sub_force_env, &
998 force_env%qmmm_env, force_env%qmmmx_env, force_env%fp_env, &
999 force_env%force_env_section, force_env%eip_env, force_env%mixed_env, &
1000 force_env%embed_env, force_env%pwdft_env, force_env%nnp_env, &
1001 force_env%root_section)
1002 last_force_env_id = last_force_env_id + 1
1003 force_env%ref_count = 1
1004 force_env%in_use = 0
1005 force_env%additional_potential = 0.0_dp
1006
1007 force_env%globenv => globenv
1008 CALL globenv_retain(force_env%globenv)
1009
1010 force_env%root_section => root_section
1011 CALL section_vals_retain(root_section)
1012
1013 force_env%para_env => para_env
1014 CALL force_env%para_env%retain()
1015
1016 CALL section_vals_retain(force_env_section)
1017 force_env%force_env_section => force_env_section
1018
1019 IF (PRESENT(fist_env)) THEN
1020 cpassert(ASSOCIATED(fist_env))
1021 cpassert(force_env%in_use == 0)
1022 force_env%in_use = use_fist_force
1023 force_env%fist_env => fist_env
1024 END IF
1025 IF (PRESENT(eip_env)) THEN
1026 cpassert(ASSOCIATED(eip_env))
1027 cpassert(force_env%in_use == 0)
1028 force_env%in_use = use_eip_force
1029 force_env%eip_env => eip_env
1030 END IF
1031 IF (PRESENT(pwdft_env)) THEN
1032 cpassert(ASSOCIATED(pwdft_env))
1033 cpassert(force_env%in_use == 0)
1034 force_env%in_use = use_pwdft_force
1035 force_env%pwdft_env => pwdft_env
1036 END IF
1037 IF (PRESENT(qs_env)) THEN
1038 cpassert(ASSOCIATED(qs_env))
1039 cpassert(force_env%in_use == 0)
1040 force_env%in_use = use_qs_force
1041 force_env%qs_env => qs_env
1042 END IF
1043 IF (PRESENT(qmmm_env)) THEN
1044 cpassert(ASSOCIATED(qmmm_env))
1045 cpassert(force_env%in_use == 0)
1046 force_env%in_use = use_qmmm
1047 force_env%qmmm_env => qmmm_env
1048 END IF
1049 IF (PRESENT(qmmmx_env)) THEN
1050 cpassert(ASSOCIATED(qmmmx_env))
1051 cpassert(force_env%in_use == 0)
1052 force_env%in_use = use_qmmmx
1053 force_env%qmmmx_env => qmmmx_env
1054 END IF
1055 IF (PRESENT(mixed_env)) THEN
1056 cpassert(ASSOCIATED(mixed_env))
1057 cpassert(force_env%in_use == 0)
1058 force_env%in_use = use_mixed_force
1059 force_env%mixed_env => mixed_env
1060 END IF
1061 IF (PRESENT(embed_env)) THEN
1062 cpassert(ASSOCIATED(embed_env))
1063 cpassert(force_env%in_use == 0)
1064 force_env%in_use = use_embed
1065 force_env%embed_env => embed_env
1066 END IF
1067 IF (PRESENT(nnp_env)) THEN
1068 cpassert(ASSOCIATED(nnp_env))
1069 cpassert(force_env%in_use == 0)
1070 force_env%in_use = use_nnp_force
1071 force_env%nnp_env => nnp_env
1072 END IF
1073 IF (PRESENT(ipi_env)) THEN
1074 cpassert(ASSOCIATED(ipi_env))
1075 cpassert(force_env%in_use == 0)
1076 force_env%in_use = use_ipi
1077 force_env%ipi_env => ipi_env
1078 END IF
1079 cpassert(force_env%in_use /= 0)
1080
1081 IF (PRESENT(sub_force_env)) THEN
1082 force_env%sub_force_env => sub_force_env
1083 END IF
1084
1085 IF (PRESENT(meta_env)) THEN
1086 force_env%meta_env => meta_env
1087 ELSE
1088 NULLIFY (force_env%meta_env)
1089 END IF
1090
1091 END SUBROUTINE force_env_create
1092
1093! **************************************************************************************************
1094!> \brief ****f* force_env_methods/mixed_energy_forces [1.0]
1095!>
1096!> Computes energy and forces for a mixed force_env type
1097!> \param force_env the force_env that holds the mixed_env type
1098!> \param calculate_forces decides if forces should be calculated
1099!> \par History
1100!> 11.06 created [fschiff]
1101!> 04.07 generalization to an illimited number of force_eval [tlaino]
1102!> 04.07 further generalization to force_eval with different geometrical
1103!> structures [tlaino]
1104!> 04.08 reorganizing the genmix structure (collecting common code)
1105!> 01.16 added CDFT [Nico Holmberg]
1106!> 08.17 added DFT embedding [Vladimir Rybkin]
1107!> \author Florian Schiffmann
1108! **************************************************************************************************
1109 SUBROUTINE mixed_energy_forces(force_env, calculate_forces)
1110
1111 TYPE(force_env_type), POINTER :: force_env
1112 LOGICAL, INTENT(IN) :: calculate_forces
1113
1114 CHARACTER(LEN=default_path_length) :: coupling_function
1115 CHARACTER(LEN=default_string_length) :: def_error, description, this_error
1116 INTEGER :: iforce_eval, iparticle, istate(2), &
1117 jparticle, mixing_type, my_group, &
1118 natom, nforce_eval, source, unit_nr
1119 INTEGER, DIMENSION(:), POINTER :: glob_natoms, itmplist, map_index
1120 LOGICAL :: dip_exists
1121 REAL(kind=dp) :: coupling_parameter, dedf, der_1, der_2, &
1122 dx, energy, err, lambda, lerr, &
1123 restraint_strength, restraint_target, &
1124 sd
1125 REAL(kind=dp), DIMENSION(3) :: dip_mix
1126 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1127 TYPE(cell_type), POINTER :: cell_mix
1128 TYPE(cp_logger_type), POINTER :: logger, my_logger
1129 TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1130 TYPE(cp_result_type), POINTER :: loc_results, results_mix
1131 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1132 TYPE(cp_subsys_type), POINTER :: subsys_mix
1133 TYPE(mixed_energy_type), POINTER :: mixed_energy
1134 TYPE(mixed_force_type), DIMENSION(:), POINTER :: global_forces
1135 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1136 TYPE(particle_list_type), POINTER :: particles_mix
1137 TYPE(section_vals_type), POINTER :: force_env_section, gen_section, &
1138 mapping_section, mixed_section, &
1139 root_section
1140 TYPE(virial_p_type), DIMENSION(:), POINTER :: virials
1141 TYPE(virial_type), POINTER :: loc_virial, virial_mix
1142
1143 logger => cp_get_default_logger()
1144 cpassert(ASSOCIATED(force_env))
1145 ! Get infos about the mixed subsys
1146 CALL force_env_get(force_env=force_env, &
1147 subsys=subsys_mix, &
1148 force_env_section=force_env_section, &
1149 root_section=root_section, &
1150 cell=cell_mix)
1151 CALL cp_subsys_get(subsys=subsys_mix, &
1152 particles=particles_mix, &
1153 virial=virial_mix, &
1154 results=results_mix)
1155 NULLIFY (map_index, glob_natoms, global_forces, itmplist)
1156
1157 nforce_eval = SIZE(force_env%sub_force_env)
1158 mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
1159 mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
1160 ! Global Info
1161 ALLOCATE (subsystems(nforce_eval))
1162 ALLOCATE (particles(nforce_eval))
1163 ! Local Info to sync
1164 ALLOCATE (global_forces(nforce_eval))
1165 ALLOCATE (energies(nforce_eval))
1166 ALLOCATE (glob_natoms(nforce_eval))
1167 ALLOCATE (virials(nforce_eval))
1168 ALLOCATE (results(nforce_eval))
1169 energies = 0.0_dp
1170 glob_natoms = 0
1171 ! Check if mixed CDFT calculation is requested and initialize
1172 CALL mixed_cdft_init(force_env, calculate_forces)
1173
1174 !
1175 IF (.NOT. force_env%mixed_env%do_mixed_cdft) THEN
1176 DO iforce_eval = 1, nforce_eval
1177 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1178 NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
1179 ALLOCATE (virials(iforce_eval)%virial)
1180 CALL cp_result_create(results(iforce_eval)%results)
1181 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1182 ! From this point on the error is the sub_error
1183 my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
1184 my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
1185 ! Copy iterations info (they are updated only in the main mixed_env)
1186 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1187 CALL cp_add_default_logger(my_logger)
1188
1189 ! Get all available subsys
1190 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1191 subsys=subsystems(iforce_eval)%subsys)
1192
1193 ! all force_env share the same cell
1194 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
1195
1196 ! Get available particles
1197 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1198 particles=particles(iforce_eval)%list)
1199
1200 ! Get Mapping index array
1201 natom = SIZE(particles(iforce_eval)%list%els)
1202
1203 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1204 map_index)
1205
1206 ! Mapping particles from iforce_eval environment to the mixed env
1207 DO iparticle = 1, natom
1208 jparticle = map_index(iparticle)
1209 particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
1210 END DO
1211
1212 ! Calculate energy and forces for each sub_force_env
1213 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1214 calc_force=calculate_forces, &
1215 skip_external_control=.true.)
1216
1217 ! Only the rank 0 process collect info for each computation
1218 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1219 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1220 potential_energy=energy)
1221 CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1222 virial=loc_virial, results=loc_results)
1223 energies(iforce_eval) = energy
1224 glob_natoms(iforce_eval) = natom
1225 virials(iforce_eval)%virial = loc_virial
1226 CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1227 END IF
1228 ! Deallocate map_index array
1229 IF (ASSOCIATED(map_index)) THEN
1230 DEALLOCATE (map_index)
1231 END IF
1232 CALL cp_rm_default_logger()
1233 END DO
1234 ELSE
1235 CALL mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
1236 glob_natoms, virials, results)
1237 END IF
1238 ! Handling Parallel execution
1239 CALL force_env%para_env%sync()
1240 ! Post CDFT operations
1241 CALL mixed_cdft_post_energy_forces(force_env)
1242 ! Let's transfer energy, natom, forces, virials
1243 CALL force_env%para_env%sum(energies)
1244 CALL force_env%para_env%sum(glob_natoms)
1245 ! Transfer forces
1246 DO iforce_eval = 1, nforce_eval
1247 ALLOCATE (global_forces(iforce_eval)%forces(3, glob_natoms(iforce_eval)))
1248 global_forces(iforce_eval)%forces = 0.0_dp
1249 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN
1250 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1251 ! Forces
1252 DO iparticle = 1, glob_natoms(iforce_eval)
1253 global_forces(iforce_eval)%forces(:, iparticle) = &
1254 particles(iforce_eval)%list%els(iparticle)%f
1255 END DO
1256 END IF
1257 END IF
1258 CALL force_env%para_env%sum(global_forces(iforce_eval)%forces)
1259 !Transfer only the relevant part of the virial..
1260 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_total)
1261 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_kinetic)
1262 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_virial)
1263 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_xc)
1264 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_fock_4c)
1265 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_constraint)
1266 !Transfer results
1267 source = 0
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()) &
1270 source = force_env%para_env%mepos
1271 END IF
1272 CALL force_env%para_env%sum(source)
1273 CALL cp_results_mp_bcast(results(iforce_eval)%results, source, force_env%para_env)
1274 END DO
1275
1276 force_env%mixed_env%energies = energies
1277 ! Start combining the different sub_force_env
1278 CALL get_mixed_env(mixed_env=force_env%mixed_env, &
1279 mixed_energy=mixed_energy)
1280
1281 !NB: do this for all MIXING_TYPE values, since some need it (e.g. linear mixing
1282 !NB if the first system has fewer atoms than the second)
1283 DO iparticle = 1, SIZE(particles_mix%els)
1284 particles_mix%els(iparticle)%f(:) = 0.0_dp
1285 END DO
1286
1287 CALL section_vals_val_get(mixed_section, "MIXING_TYPE", i_val=mixing_type)
1288 SELECT CASE (mixing_type)
1289 CASE (mix_linear_combination)
1290 ! Support offered only 2 force_eval
1291 cpassert(nforce_eval == 2)
1292 CALL section_vals_val_get(mixed_section, "LINEAR%LAMBDA", r_val=lambda)
1293 mixed_energy%pot = lambda*energies(1) + (1.0_dp - lambda)*energies(2)
1294 ! General Mapping of forces...
1295 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1296 lambda, 1, nforce_eval, map_index, mapping_section, .true.)
1297 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1298 (1.0_dp - lambda), 2, nforce_eval, map_index, mapping_section, .false.)
1299 CASE (mix_minimum)
1300 ! Support offered only 2 force_eval
1301 cpassert(nforce_eval == 2)
1302 IF (energies(1) < energies(2)) THEN
1303 mixed_energy%pot = energies(1)
1304 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1305 1.0_dp, 1, nforce_eval, map_index, mapping_section, .true.)
1306 ELSE
1307 mixed_energy%pot = energies(2)
1308 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1309 1.0_dp, 2, nforce_eval, map_index, mapping_section, .true.)
1310 END IF
1311 CASE (mix_coupled)
1312 ! Support offered only 2 force_eval
1313 cpassert(nforce_eval == 2)
1314 CALL section_vals_val_get(mixed_section, "COUPLING%COUPLING_PARAMETER", &
1315 r_val=coupling_parameter)
1316 sd = sqrt((energies(1) - energies(2))**2 + 4.0_dp*coupling_parameter**2)
1317 der_1 = (1.0_dp - (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1318 der_2 = (1.0_dp + (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1319 mixed_energy%pot = (energies(1) + energies(2) - sd)/2.0_dp
1320 ! General Mapping of forces...
1321 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1322 der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1323 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1324 der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1325 CASE (mix_restrained)
1326 ! Support offered only 2 force_eval
1327 cpassert(nforce_eval == 2)
1328 CALL section_vals_val_get(mixed_section, "RESTRAINT%RESTRAINT_TARGET", &
1329 r_val=restraint_target)
1330 CALL section_vals_val_get(mixed_section, "RESTRAINT%RESTRAINT_STRENGTH", &
1331 r_val=restraint_strength)
1332 mixed_energy%pot = energies(1) + restraint_strength*(energies(1) - energies(2) - restraint_target)**2
1333 der_2 = -2.0_dp*restraint_strength*(energies(1) - energies(2) - restraint_target)
1334 der_1 = 1.0_dp - der_2
1335 ! General Mapping of forces...
1336 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1337 der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1338 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1339 der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1340 CASE (mix_generic)
1341 ! Support any number of force_eval sections
1342 gen_section => section_vals_get_subs_vals(mixed_section, "GENERIC")
1343 CALL get_generic_info(gen_section, "MIXING_FUNCTION", coupling_function, force_env%mixed_env%par, &
1344 force_env%mixed_env%val, energies)
1345 CALL initf(1)
1346 CALL parsef(1, trim(coupling_function), force_env%mixed_env%par)
1347 ! Now the hardest part.. map energy with corresponding force_eval
1348 mixed_energy%pot = evalf(1, force_env%mixed_env%val)
1349 cpassert(evalerrtype <= 0)
1350 CALL zero_virial(virial_mix, reset=.false.)
1351 CALL cp_results_erase(results_mix)
1352 DO iforce_eval = 1, nforce_eval
1353 CALL section_vals_val_get(gen_section, "DX", r_val=dx)
1354 CALL section_vals_val_get(gen_section, "ERROR_LIMIT", r_val=lerr)
1355 dedf = evalfd(1, iforce_eval, force_env%mixed_env%val, dx, err)
1356 IF (abs(err) > lerr) THEN
1357 WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
1358 WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
1359 CALL compress(this_error, .true.)
1360 CALL compress(def_error, .true.)
1361 CALL cp_warn(__location__, &
1362 'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
1363 ' Error '//trim(this_error)//' in computing numerical derivatives larger then'// &
1364 trim(def_error)//' .')
1365 END IF
1366 ! General Mapping of forces...
1367 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1368 dedf, iforce_eval, nforce_eval, map_index, mapping_section, .false.)
1369 force_env%mixed_env%val(iforce_eval) = energies(iforce_eval)
1370 END DO
1371 ! Let's store the needed information..
1372 force_env%mixed_env%dx = dx
1373 force_env%mixed_env%lerr = lerr
1374 force_env%mixed_env%coupling_function = trim(coupling_function)
1375 CALL finalizef()
1376 CASE (mix_cdft)
1377 ! Supports any number of force_evals for calculation of CDFT properties, but forces only from two
1378 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LAMBDA", r_val=lambda)
1379 ! Get the states which determine the forces
1380 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%FORCE_STATES", i_vals=itmplist)
1381 IF (SIZE(itmplist) /= 2) &
1382 CALL cp_abort(__location__, &
1383 "Keyword FORCE_STATES takes exactly two input values.")
1384 IF (any(itmplist < 0)) &
1385 cpabort("Invalid force_eval index.")
1386 istate = itmplist
1387 IF (istate(1) > nforce_eval .OR. istate(2) > nforce_eval) &
1388 cpabort("Invalid force_eval index.")
1389 mixed_energy%pot = lambda*energies(istate(1)) + (1.0_dp - lambda)*energies(istate(2))
1390 ! General Mapping of forces...
1391 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1392 lambda, istate(1), nforce_eval, map_index, mapping_section, .true.)
1393 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1394 (1.0_dp - lambda), istate(2), nforce_eval, map_index, mapping_section, .false.)
1395 CASE DEFAULT
1396 cpabort("")
1397 END SELECT
1398 !Simply deallocate and loose the pointer references..
1399 DO iforce_eval = 1, nforce_eval
1400 DEALLOCATE (global_forces(iforce_eval)%forces)
1401 IF (ASSOCIATED(virials(iforce_eval)%virial)) DEALLOCATE (virials(iforce_eval)%virial)
1402 CALL cp_result_release(results(iforce_eval)%results)
1403 END DO
1404 DEALLOCATE (global_forces)
1405 DEALLOCATE (subsystems)
1406 DEALLOCATE (particles)
1407 DEALLOCATE (energies)
1408 DEALLOCATE (glob_natoms)
1409 DEALLOCATE (virials)
1410 DEALLOCATE (results)
1411 ! Print Section
1412 unit_nr = cp_print_key_unit_nr(logger, mixed_section, "PRINT%DIPOLE", &
1413 extension=".data", middle_name="MIXED_DIPOLE", log_filename=.false.)
1414 IF (unit_nr > 0) THEN
1415 description = '[DIPOLE]'
1416 dip_exists = test_for_result(results=results_mix, description=description)
1417 IF (dip_exists) THEN
1418 CALL get_results(results=results_mix, description=description, values=dip_mix)
1419 WRITE (unit_nr, '(/,1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE ( A.U.)|", dip_mix
1420 WRITE (unit_nr, '( 1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE (Debye)|", dip_mix*debye
1421 ELSE
1422 WRITE (unit_nr, *) "NO FORCE_EVAL section calculated the dipole"
1423 END IF
1424 END IF
1425 CALL cp_print_key_finished_output(unit_nr, logger, mixed_section, "PRINT%DIPOLE")
1426 END SUBROUTINE mixed_energy_forces
1427
1428! **************************************************************************************************
1429!> \brief Driver routine for mixed CDFT energy and force calculations
1430!> \param force_env the force_env that holds the mixed_env
1431!> \param calculate_forces if forces should be calculated
1432!> \param particles system particles
1433!> \param energies the energies of the CDFT states
1434!> \param glob_natoms the total number of particles
1435!> \param virials the virials stored in subsys
1436!> \param results results stored in subsys
1437!> \par History
1438!> 01.17 created [Nico Holmberg]
1439!> \author Nico Holmberg
1440! **************************************************************************************************
1441 SUBROUTINE mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
1442 glob_natoms, virials, results)
1443 TYPE(force_env_type), POINTER :: force_env
1444 LOGICAL, INTENT(IN) :: calculate_forces
1445 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1446 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1447 INTEGER, DIMENSION(:), POINTER :: glob_natoms
1448 TYPE(virial_p_type), DIMENSION(:), POINTER :: virials
1449 TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1450
1451 INTEGER :: iforce_eval, iparticle, jparticle, &
1452 my_group, natom, nforce_eval
1453 INTEGER, DIMENSION(:), POINTER :: map_index
1454 REAL(kind=dp) :: energy
1455 TYPE(cell_type), POINTER :: cell_mix
1456 TYPE(cp_logger_type), POINTER :: logger, my_logger
1457 TYPE(cp_result_type), POINTER :: loc_results, results_mix
1458 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1459 TYPE(cp_subsys_type), POINTER :: subsys_mix
1460 TYPE(particle_list_type), POINTER :: particles_mix
1461 TYPE(section_vals_type), POINTER :: force_env_section, mapping_section, &
1462 mixed_section, root_section
1463 TYPE(virial_type), POINTER :: loc_virial, virial_mix
1464
1465 logger => cp_get_default_logger()
1466 cpassert(ASSOCIATED(force_env))
1467 ! Get infos about the mixed subsys
1468 CALL force_env_get(force_env=force_env, &
1469 subsys=subsys_mix, &
1470 force_env_section=force_env_section, &
1471 root_section=root_section, &
1472 cell=cell_mix)
1473 CALL cp_subsys_get(subsys=subsys_mix, &
1474 particles=particles_mix, &
1475 virial=virial_mix, &
1476 results=results_mix)
1477 NULLIFY (map_index)
1478 nforce_eval = SIZE(force_env%sub_force_env)
1479 mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
1480 mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
1481 ALLOCATE (subsystems(nforce_eval))
1482 DO iforce_eval = 1, nforce_eval
1483 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1484 NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
1485 ALLOCATE (virials(iforce_eval)%virial)
1486 CALL cp_result_create(results(iforce_eval)%results)
1487 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1488 ! Get all available subsys
1489 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1490 subsys=subsystems(iforce_eval)%subsys)
1491
1492 ! all force_env share the same cell
1493 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
1494
1495 ! Get available particles
1496 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1497 particles=particles(iforce_eval)%list)
1498
1499 ! Get Mapping index array
1500 natom = SIZE(particles(iforce_eval)%list%els)
1501 ! Serial mode need to deallocate first
1502 IF (ASSOCIATED(map_index)) &
1503 DEALLOCATE (map_index)
1504 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1505 map_index)
1506
1507 ! Mapping particles from iforce_eval environment to the mixed env
1508 DO iparticle = 1, natom
1509 jparticle = map_index(iparticle)
1510 particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
1511 END DO
1512 ! Mixed CDFT + QMMM: Need to translate now
1513 IF (force_env%mixed_env%do_mixed_qmmm_cdft) &
1514 CALL apply_qmmm_translate(force_env%sub_force_env(iforce_eval)%force_env%qmmm_env)
1515 END DO
1516 ! For mixed CDFT calculations parallelized over CDFT states
1517 ! build weight and gradient on all processors before splitting into groups and
1518 ! starting energy calculation
1519 CALL mixed_cdft_build_weight(force_env, calculate_forces)
1520 DO iforce_eval = 1, nforce_eval
1521 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1522 ! From this point on the error is the sub_error
1523 IF (force_env%mixed_env%cdft_control%run_type == mixed_cdft_serial .AND. iforce_eval >= 2) THEN
1524 my_logger => force_env%mixed_env%cdft_control%sub_logger(iforce_eval - 1)%p
1525 ELSE
1526 my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
1527 my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
1528 END IF
1529 ! Copy iterations info (they are updated only in the main mixed_env)
1530 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1531 CALL cp_add_default_logger(my_logger)
1532 ! Serial CDFT calculation: transfer weight/gradient
1533 CALL mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
1534 ! Calculate energy and forces for each sub_force_env
1535 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1536 calc_force=calculate_forces, &
1537 skip_external_control=.true.)
1538 ! Only the rank 0 process collect info for each computation
1539 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1540 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1541 potential_energy=energy)
1542 CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1543 virial=loc_virial, results=loc_results)
1544 energies(iforce_eval) = energy
1545 glob_natoms(iforce_eval) = natom
1546 virials(iforce_eval)%virial = loc_virial
1547 CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1548 END IF
1549 ! Deallocate map_index array
1550 IF (ASSOCIATED(map_index)) THEN
1551 DEALLOCATE (map_index)
1552 END IF
1553 CALL cp_rm_default_logger()
1554 END DO
1555 DEALLOCATE (subsystems)
1556
1557 END SUBROUTINE mixed_cdft_energy_forces
1558
1559! **************************************************************************************************
1560!> \brief Perform additional tasks for mixed CDFT calculations after solving the electronic structure
1561!> of both CDFT states
1562!> \param force_env the force_env that holds the CDFT states
1563!> \par History
1564!> 01.17 created [Nico Holmberg]
1565!> \author Nico Holmberg
1566! **************************************************************************************************
1567 SUBROUTINE mixed_cdft_post_energy_forces(force_env)
1568 TYPE(force_env_type), POINTER :: force_env
1569
1570 INTEGER :: iforce_eval, nforce_eval, nvar
1571 TYPE(dft_control_type), POINTER :: dft_control
1572 TYPE(qs_environment_type), POINTER :: qs_env
1573
1574 cpassert(ASSOCIATED(force_env))
1575 NULLIFY (qs_env, dft_control)
1576 IF (force_env%mixed_env%do_mixed_cdft) THEN
1577 nforce_eval = SIZE(force_env%sub_force_env)
1578 nvar = force_env%mixed_env%cdft_control%nconstraint
1579 ! Transfer cdft strengths for writing restart
1580 IF (.NOT. ASSOCIATED(force_env%mixed_env%strength)) &
1581 ALLOCATE (force_env%mixed_env%strength(nforce_eval, nvar))
1582 force_env%mixed_env%strength = 0.0_dp
1583 DO iforce_eval = 1, nforce_eval
1584 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1585 IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1586 qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
1587 ELSE
1588 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1589 END IF
1590 CALL get_qs_env(qs_env, dft_control=dft_control)
1591 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) &
1592 force_env%mixed_env%strength(iforce_eval, :) = dft_control%qs_control%cdft_control%strength(:)
1593 END DO
1594 CALL force_env%para_env%sum(force_env%mixed_env%strength)
1595 ! Mixed CDFT: calculate ET coupling
1596 IF (force_env%mixed_env%do_mixed_et) THEN
1597 IF (modulo(force_env%mixed_env%cdft_control%sim_step, force_env%mixed_env%et_freq) == 0) &
1598 CALL mixed_cdft_calculate_coupling(force_env)
1599 END IF
1600 END IF
1601
1602 END SUBROUTINE mixed_cdft_post_energy_forces
1603
1604! **************************************************************************************************
1605!> \brief Computes the total energy for an embedded calculation
1606!> \param force_env ...
1607!> \author Vladimir Rybkin
1608! **************************************************************************************************
1609 SUBROUTINE embed_energy(force_env)
1610
1611 TYPE(force_env_type), POINTER :: force_env
1612
1613 INTEGER :: iforce_eval, iparticle, jparticle, &
1614 my_group, natom, nforce_eval
1615 INTEGER, DIMENSION(:), POINTER :: glob_natoms, map_index
1616 LOGICAL :: converged_embed
1617 REAL(kind=dp) :: energy
1618 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1619 TYPE(cell_type), POINTER :: cell_embed
1620 TYPE(cp_logger_type), POINTER :: logger, my_logger
1621 TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1622 TYPE(cp_result_type), POINTER :: loc_results, results_embed
1623 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1624 TYPE(cp_subsys_type), POINTER :: subsys_embed
1625 TYPE(dft_control_type), POINTER :: dft_control
1626 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1627 TYPE(particle_list_type), POINTER :: particles_embed
1628 TYPE(pw_env_type), POINTER :: pw_env
1629 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1630 TYPE(pw_r3d_rs_type), POINTER :: embed_pot, spin_embed_pot
1631 TYPE(section_vals_type), POINTER :: embed_section, force_env_section, &
1632 mapping_section, root_section
1633
1634 logger => cp_get_default_logger()
1635 cpassert(ASSOCIATED(force_env))
1636 ! Get infos about the embedding subsys
1637 CALL force_env_get(force_env=force_env, &
1638 subsys=subsys_embed, &
1639 force_env_section=force_env_section, &
1640 root_section=root_section, &
1641 cell=cell_embed)
1642 CALL cp_subsys_get(subsys=subsys_embed, &
1643 particles=particles_embed, &
1644 results=results_embed)
1645 NULLIFY (map_index, glob_natoms)
1646
1647 nforce_eval = SIZE(force_env%sub_force_env)
1648 embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1649 mapping_section => section_vals_get_subs_vals(embed_section, "MAPPING")
1650 ! Global Info
1651 ALLOCATE (subsystems(nforce_eval))
1652 ALLOCATE (particles(nforce_eval))
1653 ! Local Info to sync
1654 ALLOCATE (energies(nforce_eval))
1655 ALLOCATE (glob_natoms(nforce_eval))
1656 ALLOCATE (results(nforce_eval))
1657 energies = 0.0_dp
1658 glob_natoms = 0
1659
1660 DO iforce_eval = 1, nforce_eval
1661 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1662 NULLIFY (results(iforce_eval)%results)
1663 CALL cp_result_create(results(iforce_eval)%results)
1664 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1665 ! From this point on the error is the sub_error
1666 my_group = force_env%embed_env%group_distribution(force_env%para_env%mepos)
1667 my_logger => force_env%embed_env%sub_logger(my_group + 1)%p
1668 ! Copy iterations info (they are updated only in the main embed_env)
1669 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1670 CALL cp_add_default_logger(my_logger)
1671
1672 ! Get all available subsys
1673 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1674 subsys=subsystems(iforce_eval)%subsys)
1675
1676 ! Check if we import density from previous force calculations
1677 ! Only for QUICKSTEP
1678 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env%qs_env)) THEN
1679 NULLIFY (dft_control)
1680 CALL get_qs_env(force_env%sub_force_env(iforce_eval)%force_env%qs_env, dft_control=dft_control)
1681 IF (dft_control%qs_control%ref_embed_subsys) THEN
1682 IF (iforce_eval == 2) cpabort("Density importing force_eval can't be the first.")
1683 END IF
1684 END IF
1685
1686 ! all force_env share the same cell
1687 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_embed)
1688
1689 ! Get available particles
1690 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1691 particles=particles(iforce_eval)%list)
1692
1693 ! Get Mapping index array
1694 natom = SIZE(particles(iforce_eval)%list%els)
1695
1696 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1697 map_index, .true.)
1698
1699 ! Mapping particles from iforce_eval environment to the embed env
1700 DO iparticle = 1, natom
1701 jparticle = map_index(iparticle)
1702 particles(iforce_eval)%list%els(iparticle)%r = particles_embed%els(jparticle)%r
1703 END DO
1704
1705 ! Calculate energy and forces for each sub_force_env
1706 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1707 calc_force=.false., &
1708 skip_external_control=.true.)
1709
1710 ! Call DFT embedding
1711 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env%qs_env)) THEN
1712 NULLIFY (dft_control)
1713 CALL get_qs_env(force_env%sub_force_env(iforce_eval)%force_env%qs_env, dft_control=dft_control)
1714 IF (dft_control%qs_control%ref_embed_subsys) THEN
1715 ! Now we can optimize the embedding potential
1716 CALL dft_embedding(force_env, iforce_eval, energies, converged_embed)
1717 IF (.NOT. converged_embed) cpabort("Embedding potential optimization not converged.")
1718 END IF
1719 ! Deallocate embedding potential on the high-level subsystem
1720 IF (dft_control%qs_control%high_level_embed_subsys) THEN
1721 CALL get_qs_env(qs_env=force_env%sub_force_env(iforce_eval)%force_env%qs_env, &
1722 embed_pot=embed_pot, spin_embed_pot=spin_embed_pot, pw_env=pw_env)
1723 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1724 CALL auxbas_pw_pool%give_back_pw(embed_pot)
1725 IF (ASSOCIATED(embed_pot)) THEN
1726 CALL embed_pot%release()
1727 DEALLOCATE (embed_pot)
1728 END IF
1729 IF (ASSOCIATED(spin_embed_pot)) THEN
1730 CALL auxbas_pw_pool%give_back_pw(spin_embed_pot)
1731 CALL spin_embed_pot%release()
1732 DEALLOCATE (spin_embed_pot)
1733 END IF
1734 END IF
1735 END IF
1736
1737 ! Only the rank 0 process collect info for each computation
1738 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1739 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1740 potential_energy=energy)
1741 CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1742 results=loc_results)
1743 energies(iforce_eval) = energy
1744 glob_natoms(iforce_eval) = natom
1745 CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1746 END IF
1747 ! Deallocate map_index array
1748 IF (ASSOCIATED(map_index)) THEN
1749 DEALLOCATE (map_index)
1750 END IF
1751 CALL cp_rm_default_logger()
1752 END DO
1753
1754 ! Handling Parallel execution
1755 CALL force_env%para_env%sync()
1756 ! Let's transfer energy, natom
1757 CALL force_env%para_env%sum(energies)
1758 CALL force_env%para_env%sum(glob_natoms)
1759
1760 force_env%embed_env%energies = energies
1761
1762 !NB if the first system has fewer atoms than the second)
1763 DO iparticle = 1, SIZE(particles_embed%els)
1764 particles_embed%els(iparticle)%f(:) = 0.0_dp
1765 END DO
1766
1767 ! ONIOM type of mixing in embedding: E = E_total + E_cluster_high - E_cluster
1768 force_env%embed_env%pot_energy = energies(3) + energies(4) - energies(2)
1769
1770 !Simply deallocate and loose the pointer references..
1771 DO iforce_eval = 1, nforce_eval
1772 CALL cp_result_release(results(iforce_eval)%results)
1773 END DO
1774 DEALLOCATE (subsystems)
1775 DEALLOCATE (particles)
1776 DEALLOCATE (energies)
1777 DEALLOCATE (glob_natoms)
1778 DEALLOCATE (results)
1779
1780 END SUBROUTINE embed_energy
1781
1782! **************************************************************************************************
1783!> \brief ...
1784!> \param force_env ...
1785!> \param ref_subsys_number ...
1786!> \param energies ...
1787!> \param converged_embed ...
1788! **************************************************************************************************
1789 SUBROUTINE dft_embedding(force_env, ref_subsys_number, energies, converged_embed)
1790 TYPE(force_env_type), POINTER :: force_env
1791 INTEGER :: ref_subsys_number
1792 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1793 LOGICAL :: converged_embed
1794
1795 INTEGER :: embed_method
1796 TYPE(section_vals_type), POINTER :: embed_section, force_env_section
1797
1798 ! Find out which embedding scheme is used
1799 CALL force_env_get(force_env=force_env, &
1800 force_env_section=force_env_section)
1801 embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1802
1803 CALL section_vals_val_get(embed_section, "EMBED_METHOD", i_val=embed_method)
1804 SELECT CASE (embed_method)
1805 CASE (dfet)
1806 ! Density functional embedding
1807 CALL dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1808 CASE (dmfet)
1809 ! Density matrix embedding theory
1810 CALL dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1811 END SELECT
1812
1813 END SUBROUTINE dft_embedding
1814! **************************************************************************************************
1815!> \brief ... Main driver for DFT embedding
1816!> \param force_env ...
1817!> \param ref_subsys_number ...
1818!> \param energies ...
1819!> \param converged_embed ...
1820!> \author Vladimir Rybkin
1821! **************************************************************************************************
1822 SUBROUTINE dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1823 TYPE(force_env_type), POINTER :: force_env
1824 INTEGER :: ref_subsys_number
1825 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1826 LOGICAL :: converged_embed
1827
1828 CHARACTER(LEN=*), PARAMETER :: routinen = 'dfet_embedding'
1829
1830 INTEGER :: cluster_subsys_num, handle, &
1831 i_force_eval, i_iter, i_spin, &
1832 nforce_eval, nspins, nspins_subsys, &
1833 output_unit
1834 REAL(kind=dp) :: cluster_energy
1835 REAL(kind=dp), DIMENSION(:), POINTER :: rhs
1836 TYPE(cp_logger_type), POINTER :: logger
1837 TYPE(dft_control_type), POINTER :: dft_control
1838 TYPE(opt_embed_pot_type) :: opt_embed
1839 TYPE(pw_env_type), POINTER :: pw_env
1840 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1841 TYPE(pw_r3d_rs_type) :: diff_rho_r, diff_rho_spin
1842 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_ref, rho_r_subsys
1843 TYPE(pw_r3d_rs_type), POINTER :: embed_pot, embed_pot_subsys, &
1844 spin_embed_pot, spin_embed_pot_subsys
1845 TYPE(qs_energy_type), POINTER :: energy
1846 TYPE(qs_rho_type), POINTER :: rho, subsys_rho
1847 TYPE(section_vals_type), POINTER :: dft_section, embed_section, &
1848 force_env_section, input, &
1849 mapping_section, opt_embed_section
1850
1851 CALL timeset(routinen, handle)
1852
1853 CALL cite_reference(huang2011)
1854 CALL cite_reference(heaton_burgess2007)
1855
1856 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
1857
1858 ! Reveal input file
1859 NULLIFY (logger)
1860 logger => cp_get_default_logger()
1861 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
1862 extension=".Log")
1863
1864 NULLIFY (dft_section, input, opt_embed_section)
1865 NULLIFY (energy, dft_control)
1866
1867 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1868 pw_env=pw_env, dft_control=dft_control, rho=rho, energy=energy, &
1869 input=input)
1870 nspins = dft_control%nspins
1871
1872 dft_section => section_vals_get_subs_vals(input, "DFT")
1873 opt_embed_section => section_vals_get_subs_vals(input, &
1874 "DFT%QS%OPT_EMBED")
1875 ! Rho_r is the reference
1876 CALL qs_rho_get(rho_struct=rho, rho_r=rho_r_ref)
1877
1878 ! We need to understand how to treat spins states
1879 CALL understand_spin_states(force_env, ref_subsys_number, opt_embed%change_spin, opt_embed%open_shell_embed, &
1880 opt_embed%all_nspins)
1881
1882 ! Prepare everything for the potential maximization
1883 CALL prepare_embed_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, opt_embed, &
1884 opt_embed_section)
1885
1886 ! Initialize embedding potential
1887 CALL init_embed_pot(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, &
1888 opt_embed%add_const_pot, opt_embed%Fermi_Amaldi, opt_embed%const_pot, &
1889 opt_embed%open_shell_embed, spin_embed_pot, &
1890 opt_embed%pot_diff, opt_embed%Coulomb_guess, opt_embed%grid_opt)
1891
1892 ! Read embedding potential vector from the file
1893 IF (opt_embed%read_embed_pot .OR. opt_embed%read_embed_pot_cube) CALL read_embed_pot( &
1894 force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, spin_embed_pot, &
1895 opt_embed_section, opt_embed)
1896
1897 ! Prepare the pw object to store density differences
1898 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1899 CALL auxbas_pw_pool%create_pw(diff_rho_r)
1900 CALL pw_zero(diff_rho_r)
1901 IF (opt_embed%open_shell_embed) THEN
1902 CALL auxbas_pw_pool%create_pw(diff_rho_spin)
1903 CALL pw_zero(diff_rho_spin)
1904 END IF
1905
1906 ! Check the preliminary density differences
1907 DO i_spin = 1, nspins
1908 CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
1909 END DO
1910 IF (opt_embed%open_shell_embed) THEN ! Spin part
1911 IF (nspins == 2) THEN ! Reference systems has an open shell, else the reference diff_rho_spin is zero
1912 CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
1913 CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
1914 END IF
1915 END IF
1916
1917 DO i_force_eval = 1, ref_subsys_number - 1
1918 NULLIFY (subsys_rho, rho_r_subsys, dft_control)
1919 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, energy=energy, &
1920 dft_control=dft_control)
1921 nspins_subsys = dft_control%nspins
1922 ! Add subsystem densities
1923 CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
1924 DO i_spin = 1, nspins_subsys
1925 CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
1926 END DO
1927 IF (opt_embed%open_shell_embed) THEN ! Spin part
1928 IF (nspins_subsys == 2) THEN ! The subsystem makes contribution if it is spin-polarized
1929 ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
1930 IF ((i_force_eval == 2) .AND. (opt_embed%change_spin)) THEN
1931 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1932 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1933 ELSE
1934 ! First subsystem (always) and second subsystem (without spin change)
1935 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1936 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1937 END IF
1938 END IF
1939 END IF
1940 END DO
1941
1942 ! Print density difference
1943 CALL print_rho_diff(diff_rho_r, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1944 IF (opt_embed%open_shell_embed) THEN ! Spin part
1945 CALL print_rho_spin_diff(diff_rho_spin, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1946 END IF
1947
1948 ! Construct electrostatic guess if needed
1949 IF (opt_embed%Coulomb_guess) THEN
1950 ! Reveal resp charges for total system
1951 nforce_eval = SIZE(force_env%sub_force_env)
1952 NULLIFY (rhs)
1953 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, rhs=rhs)
1954 ! Get the mapping
1955 CALL force_env_get(force_env=force_env, &
1956 force_env_section=force_env_section)
1957 embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1958 mapping_section => section_vals_get_subs_vals(embed_section, "MAPPING")
1959
1960 DO i_force_eval = 1, ref_subsys_number - 1
1961 IF (i_force_eval == 1) THEN
1962 CALL coulomb_guess(embed_pot, rhs, mapping_section, &
1963 force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1964 ELSE
1965 CALL coulomb_guess(opt_embed%pot_diff, rhs, mapping_section, &
1966 force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1967 END IF
1968 END DO
1969 CALL pw_axpy(opt_embed%pot_diff, embed_pot)
1970 IF (.NOT. opt_embed%grid_opt) CALL pw_copy(embed_pot, opt_embed%const_pot)
1971
1972 END IF
1973
1974 ! Difference guess
1975 IF (opt_embed%diff_guess) THEN
1976 CALL pw_copy(diff_rho_r, embed_pot)
1977 IF (.NOT. opt_embed%grid_opt) CALL pw_copy(embed_pot, opt_embed%const_pot)
1978 ! Open shell
1979 IF (opt_embed%open_shell_embed) CALL pw_copy(diff_rho_spin, spin_embed_pot)
1980 END IF
1981
1982 ! Calculate subsystems with trial embedding potential
1983 DO i_iter = 1, opt_embed%n_iter
1984 opt_embed%i_iter = i_iter
1985
1986 ! Set the density difference as the negative reference one
1987 CALL pw_zero(diff_rho_r)
1988 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, dft_control=dft_control)
1989 nspins = dft_control%nspins
1990 DO i_spin = 1, nspins
1991 CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
1992 END DO
1993 IF (opt_embed%open_shell_embed) THEN ! Spin part
1994 CALL pw_zero(diff_rho_spin)
1995 IF (nspins == 2) THEN ! Reference systems has an open shell, else the reference diff_rho_spin is zero
1996 CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
1997 CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
1998 END IF
1999 END IF
2000
2001 DO i_force_eval = 1, ref_subsys_number - 1
2002 NULLIFY (dft_control)
2003 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
2004 nspins_subsys = dft_control%nspins
2005
2006 IF ((i_force_eval == 2) .AND. (opt_embed%change_spin)) THEN
2007 ! Here we change the sign of the spin embedding potential due to spin change:
2008 ! only in spin_embed_subsys
2009 CALL make_subsys_embed_pot(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2010 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
2011 opt_embed%open_shell_embed, .true.)
2012 ELSE ! Regular case
2013 CALL make_subsys_embed_pot(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2014 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
2015 opt_embed%open_shell_embed, .false.)
2016 END IF
2017
2018 ! Switch on external potential in the subsystems
2019 dft_control%apply_embed_pot = .true.
2020
2021 ! Add the embedding potential
2022 CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, embed_pot=embed_pot_subsys)
2023 IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys == 2)) THEN ! Spin part
2024 CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2025 spin_embed_pot=spin_embed_pot_subsys)
2026 END IF
2027
2028 ! Get the previous subsystem densities
2029 CALL get_prev_density(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
2030
2031 ! Calculate the new density
2032 CALL force_env_calc_energy_force(force_env=force_env%sub_force_env(i_force_eval)%force_env, &
2033 calc_force=.false., &
2034 skip_external_control=.true.)
2035
2036 CALL get_max_subsys_diff(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
2037
2038 ! Extract subsystem density and energy
2039 NULLIFY (rho_r_subsys, energy)
2040
2041 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, &
2042 energy=energy)
2043 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + energy%total
2044
2045 ! Find out which subsystem is the cluster
2046 IF (dft_control%qs_control%cluster_embed_subsys) THEN
2047 cluster_subsys_num = i_force_eval
2048 cluster_energy = energy%total
2049 END IF
2050
2051 ! Add subsystem densities
2052 CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
2053 DO i_spin = 1, nspins_subsys
2054 CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
2055 END DO
2056 IF (opt_embed%open_shell_embed) THEN ! Spin part
2057 IF (nspins_subsys == 2) THEN ! The subsystem makes contribution if it is spin-polarized
2058 ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
2059 IF ((i_force_eval == 2) .AND. (opt_embed%change_spin)) THEN
2060 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
2061 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
2062 ELSE
2063 ! First subsystem (always) and second subsystem (without spin change)
2064 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
2065 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
2066 END IF
2067 END IF
2068 END IF
2069
2070 ! Release embedding potential for subsystem
2071 CALL embed_pot_subsys%release()
2072 DEALLOCATE (embed_pot_subsys)
2073 IF (opt_embed%open_shell_embed) THEN
2074 CALL spin_embed_pot_subsys%release()
2075 DEALLOCATE (spin_embed_pot_subsys)
2076 END IF
2077
2078 END DO ! i_force_eval
2079
2080 ! Print embedding potential for restart
2081 CALL print_embed_restart(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2082 opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
2083 spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .false.)
2084 CALL print_pot_simple_grid(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2085 embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .false., &
2086 force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
2087
2088 ! Integrate the potential over density differences and add to w functional; also add regularization contribution
2089 DO i_spin = 1, nspins ! Sum over nspins for the reference system, not subsystem!
2090 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) - pw_integral_ab(embed_pot, rho_r_ref(i_spin))
2091 END DO
2092 ! Spin part
2093 IF (opt_embed%open_shell_embed) THEN
2094 ! If reference system is not spin-polarized then it does not make a contribution to W functional
2095 IF (nspins == 2) THEN
2096 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) &
2097 - pw_integral_ab(spin_embed_pot, rho_r_ref(1)) &
2098 + pw_integral_ab(spin_embed_pot, rho_r_ref(2))
2099 END IF
2100 END IF
2101 ! Finally, add the regularization term
2102 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + opt_embed%reg_term
2103
2104 ! Print information and check convergence
2105 CALL print_emb_opt_info(output_unit, i_iter, opt_embed)
2106 CALL conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit)
2107 IF (opt_embed%converged) EXIT
2108
2109 ! Update the trust radius and control the step
2110 IF ((i_iter > 1) .AND. (.NOT. opt_embed%steep_desc)) CALL step_control(opt_embed)
2111
2112 ! Print density difference
2113 CALL print_rho_diff(diff_rho_r, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
2114 IF (opt_embed%open_shell_embed) THEN ! Spin part
2115 CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
2116 END IF
2117
2118 ! Calculate potential gradient if the step has been accepted. Otherwise, we reuse the previous one
2119
2120 IF (opt_embed%accept_step .AND. (.NOT. opt_embed%grid_opt)) &
2121 CALL calculate_embed_pot_grad(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2122 diff_rho_r, diff_rho_spin, opt_embed)
2123 ! Take the embedding step
2124 CALL opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, &
2125 force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2126
2127 END DO ! i_iter
2128
2129 ! Print final embedding potential for restart
2130 CALL print_embed_restart(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2131 opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
2132 spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .true.)
2133 CALL print_pot_simple_grid(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2134 embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .true., &
2135 force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
2136
2137 ! Print final density difference
2138 !CALL print_rho_diff(diff_rho_r, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .TRUE.)
2139 IF (opt_embed%open_shell_embed) THEN ! Spin part
2140 CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .true.)
2141 END IF
2142
2143 ! Give away plane waves pools
2144 CALL diff_rho_r%release()
2145 IF (opt_embed%open_shell_embed) THEN
2146 CALL diff_rho_spin%release()
2147 END IF
2148
2149 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
2150 "PRINT%PROGRAM_RUN_INFO")
2151
2152 ! If converged send the embedding potential to the higher-level calculation.
2153 IF (opt_embed%converged) THEN
2154 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, dft_control=dft_control, &
2155 pw_env=pw_env)
2156 nspins_subsys = dft_control%nspins
2157 dft_control%apply_embed_pot = .true.
2158 ! The embedded subsystem corresponds to subsystem #2, where spin change is possible
2159 CALL make_subsys_embed_pot(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, &
2160 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
2161 opt_embed%open_shell_embed, opt_embed%change_spin)
2162
2163 IF (opt_embed%Coulomb_guess) THEN
2164 CALL pw_axpy(opt_embed%pot_diff, embed_pot_subsys, -1.0_dp, allow_noncompatible_grids=.true.)
2165 END IF
2166
2167 CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, embed_pot=embed_pot_subsys)
2168
2169 IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys == 2)) THEN
2170 CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, &
2171 spin_embed_pot=spin_embed_pot_subsys)
2172 END IF
2173
2174 ! Substitute the correct energy in energies: only on rank 0
2175 IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source()) THEN
2176 energies(cluster_subsys_num) = cluster_energy
2177 END IF
2178 END IF
2179
2180 ! Deallocate and release opt_embed content
2181 CALL release_opt_embed(opt_embed)
2182
2183 ! Deallocate embedding potential
2184 CALL embed_pot%release()
2185 DEALLOCATE (embed_pot)
2186 IF (opt_embed%open_shell_embed) THEN
2187 CALL spin_embed_pot%release()
2188 DEALLOCATE (spin_embed_pot)
2189 END IF
2190
2191 converged_embed = opt_embed%converged
2192
2193 CALL timestop(handle)
2194
2195 END SUBROUTINE dfet_embedding
2196
2197! **************************************************************************************************
2198!> \brief Main driver for the DMFET embedding
2199!> \param force_env ...
2200!> \param ref_subsys_number ...
2201!> \param energies ...
2202!> \param converged_embed ...
2203!> \author Vladimir Rybkin
2204! **************************************************************************************************
2205 SUBROUTINE dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
2206 TYPE(force_env_type), POINTER :: force_env
2207 INTEGER :: ref_subsys_number
2208 REAL(kind=dp), DIMENSION(:), POINTER :: energies
2209 LOGICAL :: converged_embed
2210
2211 CHARACTER(LEN=*), PARAMETER :: routinen = 'dmfet_embedding'
2212
2213 INTEGER :: cluster_subsys_num, handle, &
2214 i_force_eval, i_iter, output_unit
2215 LOGICAL :: subsys_open_shell
2216 REAL(kind=dp) :: cluster_energy
2217 TYPE(cp_logger_type), POINTER :: logger
2218 TYPE(dft_control_type), POINTER :: dft_control
2219 TYPE(mp_para_env_type), POINTER :: para_env
2220 TYPE(opt_dmfet_pot_type) :: opt_dmfet
2221 TYPE(qs_energy_type), POINTER :: energy
2222 TYPE(section_vals_type), POINTER :: dft_section, input, opt_dmfet_section
2223
2224 CALL timeset(routinen, handle)
2225
2226 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2227 para_env=para_env)
2228
2229 ! Reveal input file
2230 NULLIFY (logger)
2231 logger => cp_get_default_logger()
2232 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
2233 extension=".Log")
2234
2235 NULLIFY (dft_section, input, opt_dmfet_section)
2236 NULLIFY (energy)
2237
2238 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2239 energy=energy, input=input)
2240
2241 dft_section => section_vals_get_subs_vals(input, "DFT")
2242 opt_dmfet_section => section_vals_get_subs_vals(input, &
2243 "DFT%QS%OPT_DMFET")
2244
2245 ! We need to understand how to treat spins states
2246 CALL understand_spin_states(force_env, ref_subsys_number, opt_dmfet%change_spin, opt_dmfet%open_shell_embed, &
2247 opt_dmfet%all_nspins)
2248
2249 ! Prepare for the potential optimization
2250 CALL prepare_dmfet_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2251 opt_dmfet, opt_dmfet_section)
2252
2253 ! Get the reference density matrix/matrices
2254 subsys_open_shell = subsys_spin(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2255 CALL build_full_dm(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2256 opt_dmfet%dm_total, subsys_open_shell, opt_dmfet%open_shell_embed, opt_dmfet%dm_total_beta)
2257
2258 ! Check the preliminary DM difference
2259 CALL cp_fm_copy_general(opt_dmfet%dm_total, opt_dmfet%dm_diff, para_env)
2260 IF (opt_dmfet%open_shell_embed) CALL cp_fm_copy_general(opt_dmfet%dm_total_beta, &
2261 opt_dmfet%dm_diff_beta, para_env)
2262
2263 DO i_force_eval = 1, ref_subsys_number - 1
2264
2265 ! Get the subsystem density matrix/matrices
2266 subsys_open_shell = subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2267
2268 CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2269 opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2270 opt_dmfet%dm_subsys_beta)
2271
2272 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2273
2274 IF (opt_dmfet%open_shell_embed) CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, &
2275 1.0_dp, opt_dmfet%dm_subsys_beta)
2276
2277 END DO
2278
2279 ! Main loop of iterative matrix potential optimization
2280 DO i_iter = 1, opt_dmfet%n_iter
2281
2282 opt_dmfet%i_iter = i_iter
2283
2284 ! Set the dm difference as the reference one
2285 CALL cp_fm_copy_general(opt_dmfet%dm_total, opt_dmfet%dm_diff, para_env)
2286
2287 IF (opt_dmfet%open_shell_embed) CALL cp_fm_copy_general(opt_dmfet%dm_total_beta, &
2288 opt_dmfet%dm_diff_beta, para_env)
2289
2290 ! Loop over force evaluations
2291 DO i_force_eval = 1, ref_subsys_number - 1
2292
2293 ! Switch on external potential in the subsystems
2294 NULLIFY (dft_control)
2295 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
2296 dft_control%apply_dmfet_pot = .true.
2297
2298 ! Calculate the new density
2299 CALL force_env_calc_energy_force(force_env=force_env%sub_force_env(i_force_eval)%force_env, &
2300 calc_force=.false., &
2301 skip_external_control=.true.)
2302
2303 ! Extract subsystem density matrix and energy
2304 NULLIFY (energy)
2305
2306 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, energy=energy)
2307 opt_dmfet%w_func(i_iter) = opt_dmfet%w_func(i_iter) + energy%total
2308
2309 ! Find out which subsystem is the cluster
2310 IF (dft_control%qs_control%cluster_embed_subsys) THEN
2311 cluster_subsys_num = i_force_eval
2312 cluster_energy = energy%total
2313 END IF
2314
2315 ! Add subsystem density matrices
2316 subsys_open_shell = subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2317
2318 CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2319 opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2320 opt_dmfet%dm_subsys_beta)
2321
2322 IF (opt_dmfet%open_shell_embed) THEN ! Open-shell embedding
2323 ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
2324 IF ((i_force_eval == 2) .AND. (opt_dmfet%change_spin)) THEN
2325 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, 1.0_dp, opt_dmfet%dm_subsys)
2326 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys_beta)
2327 ELSE
2328 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2329 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, 1.0_dp, opt_dmfet%dm_subsys_beta)
2330 END IF
2331 ELSE ! Closed-shell embedding
2332 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2333 END IF
2334
2335 END DO ! i_force_eval
2336
2337 CALL check_dmfet(opt_dmfet, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2338
2339 END DO ! i_iter
2340
2341 ! Substitute the correct energy in energies: only on rank 0
2342 IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source()) THEN
2343 energies(cluster_subsys_num) = cluster_energy
2344 END IF
2345
2346 CALL release_dmfet_opt(opt_dmfet)
2347
2348 converged_embed = .false.
2349
2350 END SUBROUTINE dmfet_embedding
2351
2352END 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:125
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:184
cell_methods::cell_create
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition cell_methods.F:94
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:565
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:535
cell_types::cell_release
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:608
cell_types::cell_clone
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition cell_types.F:118
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:170
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:1518
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:169
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::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)
...
Definition cp_output_handling.F:853
cp_output_handling::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:73
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:1067
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
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_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:434
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:362
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:286
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:1239
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:405
eip_silicon::eip_lenosky
subroutine, public eip_lenosky(eip_env)
Interface routine of Goedecker's Lenosky force field to CP2K.
Definition eip_silicon.F:241
eip_silicon::eip_tersoff
subroutine, public eip_tersoff(eip_env)
Interface routine of the Tersoff force field to CP2K.
Definition eip_silicon.F:567
eip_silicon::eip_bazant
subroutine, public eip_bazant(eip_env)
Interface routine of Goedecker's Bazant EDIP to CP2K.
Definition eip_silicon.F:75
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:701
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:974
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:537
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:342
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:666
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
force_fields_util
Definition force_fields_util.F:19
force_fields_util::get_generic_info
subroutine, public get_generic_info(gen_section, func_name, xfunction, parameters, values, var_values, size_variables, i_rep_sec, input_variables)
Reads from the input structure all information for generic functions.
Definition force_fields_util.F:1199
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
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:122
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:17
input_constants::use_bazant_eip
integer, parameter, public use_bazant_eip
Definition input_constants.F:740
input_constants::driver_run
integer, parameter, public driver_run
Definition input_constants.F:127
input_constants::mixed_cdft_serial
integer, parameter, public mixed_cdft_serial
Definition input_constants.F:804
input_constants::mix_cdft
integer, parameter, public mix_cdft
Definition input_constants.F:813
input_constants::do_method_ofgpw
integer, parameter, public do_method_ofgpw
Definition input_constants.F:259
input_constants::mix_linear_combination
integer, parameter, public mix_linear_combination
Definition input_constants.F:813
input_constants::do_method_rigpw
integer, parameter, public do_method_rigpw
Definition input_constants.F:259
input_constants::do_method_gpw
integer, parameter, public do_method_gpw
Definition input_constants.F:259
input_constants::dmfet
integer, parameter, public dmfet
Definition input_constants.F:1307
input_constants::use_tersoff_eip
integer, parameter, public use_tersoff_eip
Definition input_constants.F:740
input_constants::use_lenosky_eip
integer, parameter, public use_lenosky_eip
Definition input_constants.F:740
input_constants::ehrenfest
integer, parameter, public ehrenfest
Definition input_constants.F:127
input_constants::use_stillinger_weber_eip
integer, parameter, public use_stillinger_weber_eip
Definition input_constants.F:740
input_constants::mix_coupled
integer, parameter, public mix_coupled
Definition input_constants.F:813
input_constants::debug_run
integer, parameter, public debug_run
Definition input_constants.F:127
input_constants::mix_restrained
integer, parameter, public mix_restrained
Definition input_constants.F:813
input_constants::do_method_gapw
integer, parameter, public do_method_gapw
Definition input_constants.F:259
input_constants::mix_generic
integer, parameter, public mix_generic
Definition input_constants.F:813
input_constants::mol_dyn_run
integer, parameter, public mol_dyn_run
Definition input_constants.F:127
input_constants::dfet
integer, parameter, public dfet
Definition input_constants.F:1307
input_constants::cell_opt_run
integer, parameter, public cell_opt_run
Definition input_constants.F:127
input_constants::do_method_lrigpw
integer, parameter, public do_method_lrigpw
Definition input_constants.F:259
input_constants::mix_minimum
integer, parameter, public mix_minimum
Definition input_constants.F:813
input_constants::do_method_gapw_xc
integer, parameter, public do_method_gapw_xc
Definition input_constants.F:259
input_constants::geo_opt_run
integer, parameter, public geo_opt_run
Definition input_constants.F:127
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:638
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:731
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:1047
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:169
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:787
kpoint_methods::kpoint_initialize
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
Definition kpoint_methods.F:121
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:599
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:640
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:371
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)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:525
kpsym
K-points and crystal symmetry routines based on.
Definition kpsym.F:28
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:440
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:158
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
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:1109
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:399
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:139
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:512
nnp_environment_types
Data types for neural network potentials.
Definition nnp_environment_types.F:13
nnp_force
Methods dealing with Neural Network potentials.
Definition nnp_force.F:13
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:62
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:3120
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:1025
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:3065
optimize_embedding_potential::get_max_subsys_diff
subroutine, public get_max_subsys_diff(opt_embed, force_env, subsys_num)
...
Definition optimize_embedding_potential.F:3447
optimize_embedding_potential::print_emb_opt_info
subroutine, public print_emb_opt_info(output_unit, step_num, opt_embed)
...
Definition optimize_embedding_potential.F:3384
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:207
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:664
optimize_embedding_potential::get_prev_density
subroutine, public get_prev_density(opt_embed, force_env, subsys_num)
...
Definition optimize_embedding_potential.F:3417
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:934
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:3014
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:3480
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:2171
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:1468
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:1073
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:3218
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:292
optimize_embedding_potential::release_opt_embed
subroutine, public release_opt_embed(opt_embed)
Deallocate stuff for optimizing embedding potential.
Definition optimize_embedding_potential.F:857
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:374
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)
Construct basis set rotation matrices.
Definition qs_basis_rotation_methods.F:52
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, 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:530
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:1106
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:111
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:53
cp_control_types::dft_control_type
Definition cp_control_types.F:744
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:121
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:171
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
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:151
nnp_environment_types::nnp_type
Main data type collecting all relevant data for neural network potentials.
Definition nnp_environment_types.F:83
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:220
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
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:103
virial_types::virial_p_type
Definition virial_types.F:52
virial_types::virial_type
Definition virial_types.F:26