(git:e7e05ae)
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-2024 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 ! **************************************************************************************************
17 MODULE 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_control_types, ONLY: dft_control_type
34  USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
35  USE cp_fm_types, ONLY: cp_fm_copy_general
36  USE cp_iter_types, ONLY: cp_iteration_info_copy_iter
37  USE cp_log_handling, ONLY: cp_add_default_logger,&
38  cp_get_default_logger,&
39  cp_logger_type,&
40  cp_rm_default_logger,&
41  cp_to_string
42  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
43  cp_print_key_unit_nr,&
44  low_print_level
45  USE cp_result_methods, ONLY: cp_results_erase,&
46  cp_results_mp_bcast,&
47  get_results,&
48  test_for_result
49  USE cp_result_types, ONLY: cp_result_copy,&
50  cp_result_create,&
51  cp_result_p_type,&
52  cp_result_release,&
53  cp_result_type
54  USE cp_subsys_types, ONLY: cp_subsys_get,&
55  cp_subsys_p_type,&
56  cp_subsys_set,&
57  cp_subsys_type
58  USE eip_environment_types, ONLY: eip_environment_type
59  USE eip_silicon, ONLY: eip_bazant,&
60  eip_lenosky
61  USE embed_types, ONLY: embed_env_type,&
62  opt_dmfet_pot_type,&
63  opt_embed_pot_type
64  USE external_potential_methods, ONLY: add_external_potential
65  USE fist_environment_types, ONLY: fist_environment_type
66  USE fist_force, ONLY: fist_calc_energy_force
67  USE force_env_types, ONLY: &
68  force_env_get, force_env_get_natom, force_env_p_type, force_env_set, force_env_type, &
69  use_eip_force, use_embed, use_fist_force, use_mixed_force, use_nnp_force, use_prog_name, &
70  use_pwdft_force, use_qmmm, use_qmmmx, use_qs_force
71  USE force_env_utils, ONLY: rescale_forces,&
72  write_atener,&
73  write_forces
74  USE force_fields_util, ONLY: get_generic_info
75  USE fp_methods, ONLY: fp_eval
76  USE fparser, ONLY: evalerrtype,&
77  evalf,&
78  evalfd,&
79  finalizef,&
80  initf,&
81  parsef
82  USE global_types, ONLY: global_environment_type,&
83  globenv_retain
84  USE grrm_utils, ONLY: write_grrm
85  USE input_constants, ONLY: &
86  debug_run, dfet, dmfet, mix_cdft, mix_coupled, mix_generic, mix_linear_combination, &
87  mix_minimum, mix_restrained, mixed_cdft_serial, use_bazant_eip, use_lenosky_eip
88  USE input_section_types, ONLY: section_vals_get_subs_vals,&
89  section_vals_retain,&
90  section_vals_type,&
91  section_vals_val_get
92  USE kahan_sum, ONLY: accurate_sum
93  USE kinds, ONLY: default_path_length,&
94  default_string_length,&
95  dp
96  USE machine, ONLY: m_memory
97  USE mathlib, ONLY: abnormal_value
98  USE message_passing, ONLY: mp_para_env_type
99  USE metadynamics_types, ONLY: meta_env_type
100  USE mixed_cdft_methods, ONLY: mixed_cdft_build_weight,&
101  mixed_cdft_calculate_coupling,&
102  mixed_cdft_init
103  USE mixed_energy_types, ONLY: mixed_energy_type,&
104  mixed_force_type
105  USE mixed_environment_types, ONLY: get_mixed_env,&
106  mixed_environment_type
107  USE mixed_environment_utils, ONLY: get_subsys_map_index,&
108  mixed_map_forces
109  USE molecule_kind_list_types, ONLY: molecule_kind_list_type
110  USE molecule_kind_types, ONLY: get_molecule_kind,&
111  molecule_kind_type
112  USE nnp_environment_types, ONLY: nnp_type
113  USE nnp_force, ONLY: nnp_calc_energy_force
114  USE optimize_dmfet_potential, ONLY: build_full_dm,&
115  check_dmfet,&
116  prepare_dmfet_opt,&
117  release_dmfet_opt,&
118  subsys_spin
119  USE optimize_embedding_potential, ONLY: &
120  coulomb_guess, calculate_embed_pot_grad, conv_check_embed, get_max_subsys_diff, &
121  get_prev_density, init_embed_pot, make_subsys_embed_pot, opt_embed_step, &
122  prepare_embed_opt, print_emb_opt_info, print_embed_restart, print_pot_simple_grid, &
123  print_rho_diff, print_rho_spin_diff, read_embed_pot, release_opt_embed, step_control, &
124  understand_spin_states
125  USE particle_list_types, ONLY: particle_list_p_type,&
126  particle_list_type
127  USE physcon, ONLY: debye
128  USE pw_env_types, ONLY: pw_env_get,&
129  pw_env_type
130  USE pw_methods, ONLY: pw_axpy,&
131  pw_copy,&
132  pw_integral_ab,&
133  pw_zero
134  USE pw_pool_types, ONLY: pw_pool_type
135  USE pw_types, ONLY: pw_r3d_rs_type
136  USE pwdft_environment, ONLY: pwdft_calc_energy_force
137  USE pwdft_environment_types, ONLY: pwdft_environment_type
138  USE qmmm_force, ONLY: qmmm_calc_energy_force
139  USE qmmm_types, ONLY: qmmm_env_type
140  USE qmmm_util, ONLY: apply_qmmm_translate
141  USE qmmmx_force, ONLY: qmmmx_calc_energy_force
142  USE qmmmx_types, ONLY: qmmmx_env_type
143  USE qs_energy_types, ONLY: qs_energy_type
144  USE qs_environment_types, ONLY: get_qs_env,&
145  qs_environment_type,&
146  set_qs_env
147  USE qs_force, ONLY: qs_calc_energy_force
148  USE qs_rho_types, ONLY: qs_rho_get,&
149  qs_rho_type
150  USE restraint, ONLY: restraint_control
151  USE scine_utils, ONLY: write_scine
152  USE string_utilities, ONLY: compress
153  USE virial_methods, ONLY: write_stress_tensor,&
154  write_stress_tensor_components
155  USE virial_types, ONLY: symmetrize_virial,&
156  virial_p_type,&
157  virial_type,&
158  zero_virial
159 #include "./base/base_uses.f90"
160 
161  IMPLICIT NONE
162 
163  PRIVATE
164 
165  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_env_methods'
166 
167  PUBLIC :: force_env_create, &
168  force_env_calc_energy_force, &
169  force_env_calc_num_pressure
170 
171  INTEGER, SAVE, PRIVATE :: last_force_env_id = 0
172 
173 CONTAINS
174 
175 ! **************************************************************************************************
176 !> \brief Interface routine for force and energy calculations
177 !> \param force_env the force_env of which you want the energy and forces
178 !> \param calc_force if false the forces *might* be left unchanged
179 !> or be invalid, no guarantees can be given. Defaults to true
180 !> \param consistent_energies Performs an additional qs_ks_update_qs_env, so
181 !> that the energies are appropriate to the forces, they are in the
182 !> non-selfconsistent case not consistent to each other! [08.2005, TdK]
183 !> \param skip_external_control ...
184 !> \param eval_energy_forces ...
185 !> \param require_consistent_energy_force ...
186 !> \param linres ...
187 !> \param calc_stress_tensor ...
188 !> \author CJM & fawzi
189 ! **************************************************************************************************
190  RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
191  consistent_energies, skip_external_control, eval_energy_forces, &
192  require_consistent_energy_force, linres, calc_stress_tensor)
193 
194  TYPE(force_env_type), POINTER :: force_env
195  LOGICAL, INTENT(IN), OPTIONAL :: calc_force, consistent_energies, skip_external_control, &
196  eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor
197 
198  REAL(kind=dp), PARAMETER :: ateps = 1.0e-6_dp
199 
200  INTEGER :: i, ikind, j, nat, ndigits, nfixed_atoms, &
201  nfixed_atoms_total, nkind, &
202  output_unit, print_forces, print_grrm, &
203  print_scine
204  LOGICAL :: calculate_forces, calculate_stress_tensor, energy_consistency, eval_ef, &
205  linres_run, my_skip, print_components
206  REAL(kind=dp) :: checksum, e_entropy, e_gap, e_pot, &
207  sum_energy, sum_pv_virial, &
208  sum_stress_tensor
209  REAL(kind=dp), DIMENSION(3) :: grand_total_force, total_force
210  REAL(kind=dp), DIMENSION(3, 3) :: atomic_stress_tensor, diff_stress_tensor
211  TYPE(atprop_type), POINTER :: atprop_env
212  TYPE(cell_type), POINTER :: cell
213  TYPE(cp_logger_type), POINTER :: logger
214  TYPE(cp_subsys_type), POINTER :: subsys
215  TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
216  TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
217  TYPE(molecule_kind_type), POINTER :: molecule_kind
218  TYPE(particle_list_type), POINTER :: core_particles, particles, &
219  shell_particles
220  TYPE(virial_type), POINTER :: virial
221 
222  NULLIFY (logger, virial, subsys, atprop_env, cell)
223  logger => cp_get_default_logger()
224  eval_ef = .true.
225  my_skip = .false.
226  calculate_forces = .true.
227  energy_consistency = .false.
228  linres_run = .false.
229  e_gap = -1.0_dp
230  e_entropy = -1.0_dp
231 
232  IF (PRESENT(eval_energy_forces)) eval_ef = eval_energy_forces
233  IF (PRESENT(skip_external_control)) my_skip = skip_external_control
234  IF (PRESENT(calc_force)) calculate_forces = calc_force
235  IF (PRESENT(calc_stress_tensor)) THEN
236  calculate_stress_tensor = calc_stress_tensor
237  ELSE
238  calculate_stress_tensor = calculate_forces
239  END IF
240  IF (PRESENT(consistent_energies)) energy_consistency = consistent_energies
241  IF (PRESENT(linres)) linres_run = linres
242 
243  cpassert(ASSOCIATED(force_env))
244  cpassert(force_env%ref_count > 0)
245  CALL force_env_get(force_env, subsys=subsys)
246  CALL force_env_set(force_env, additional_potential=0.0_dp)
247  CALL cp_subsys_get(subsys, virial=virial, atprop=atprop_env, cell=cell)
248  IF (virial%pv_availability) CALL zero_virial(virial, reset=.false.)
249 
250  nat = force_env_get_natom(force_env)
251  CALL atprop_init(atprop_env, nat)
252  IF (eval_ef) THEN
253  SELECT CASE (force_env%in_use)
254  CASE (use_fist_force)
255  CALL fist_calc_energy_force(force_env%fist_env)
256  CASE (use_qs_force)
257  CALL qs_calc_energy_force(force_env%qs_env, calculate_forces, energy_consistency, linres_run)
258  CASE (use_pwdft_force)
259  IF (virial%pv_availability .AND. calculate_stress_tensor) THEN
260  CALL pwdft_calc_energy_force(force_env%pwdft_env, calculate_forces,.NOT. virial%pv_numer)
261  ELSE
262  CALL pwdft_calc_energy_force(force_env%pwdft_env, calculate_forces, .false.)
263  END IF
264  e_gap = force_env%pwdft_env%energy%band_gap
265  e_entropy = force_env%pwdft_env%energy%entropy
266  CASE (use_eip_force)
267  IF (force_env%eip_env%eip_model == use_lenosky_eip) THEN
268  CALL eip_lenosky(force_env%eip_env)
269  ELSE IF (force_env%eip_env%eip_model == use_bazant_eip) THEN
270  CALL eip_bazant(force_env%eip_env)
271  END IF
272  CASE (use_qmmm)
273  CALL qmmm_calc_energy_force(force_env%qmmm_env, &
274  calculate_forces, energy_consistency, linres=linres_run)
275  CASE (use_qmmmx)
276  CALL qmmmx_calc_energy_force(force_env%qmmmx_env, &
277  calculate_forces, energy_consistency, linres=linres_run, &
278  require_consistent_energy_force=require_consistent_energy_force)
279  CASE (use_mixed_force)
280  CALL mixed_energy_forces(force_env, calculate_forces)
281  CASE (use_nnp_force)
282  CALL nnp_calc_energy_force(force_env%nnp_env, &
283  calculate_forces)
284  CASE (use_embed)
285  CALL embed_energy(force_env)
286  CASE default
287  cpabort("")
288  END SELECT
289  END IF
290  ! In case it is requested, we evaluate the stress tensor numerically
291  IF (virial%pv_availability) THEN
292  IF (virial%pv_numer .AND. calculate_stress_tensor) THEN
293  ! Compute the numerical stress tensor
294  CALL force_env_calc_num_pressure(force_env)
295  ELSE
296  IF (calculate_forces) THEN
297  ! Symmetrize analytical stress tensor
298  CALL symmetrize_virial(virial)
299  ELSE
300  IF (calculate_stress_tensor) THEN
301  CALL cp_warn(__location__, "The calculation of the stress tensor "// &
302  "requires the calculation of the forces")
303  END IF
304  END IF
305  END IF
306  END IF
307 
308  !sample peak memory
309  CALL m_memory()
310 
311  ! Some additional tasks..
312  IF (.NOT. my_skip) THEN
313  ! Flexible Partitioning
314  IF (ASSOCIATED(force_env%fp_env)) THEN
315  IF (force_env%fp_env%use_fp) THEN
316  CALL fp_eval(force_env%fp_env, subsys, cell)
317  END IF
318  END IF
319  ! Constraints ONLY of Fixed Atom type
320  CALL fix_atom_control(force_env)
321  ! All Restraints
322  CALL restraint_control(force_env)
323  ! Virtual Sites
324  CALL vsite_force_control(force_env)
325  ! External Potential
326  CALL add_external_potential(force_env)
327  ! Rescale forces if requested
328  CALL rescale_forces(force_env)
329  END IF
330 
331  CALL force_env_get(force_env, potential_energy=e_pot)
332 
333  ! Print always Energy in the same format for all methods
334  output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
335  extension=".Log")
336  IF (output_unit > 0) THEN
337  WRITE (output_unit, '(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) energy [a.u.]: ",T55,F26.15)') &
338  adjustr(trim(use_prog_name(force_env%in_use))), e_pot
339  IF (e_gap .GT. -.1_dp) THEN
340  WRITE (output_unit, '(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) gap [a.u.]: ",T55,F26.15)') &
341  adjustr(trim(use_prog_name(force_env%in_use))), e_gap
342  END IF
343  IF (e_entropy .GT. -0.1_dp) THEN
344  WRITE (output_unit, '(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) free energy [a.u.]: ",T55,F26.15)') &
345  adjustr(trim(use_prog_name(force_env%in_use))), e_pot - e_entropy
346  END IF
347  END IF
348  CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
349  "PRINT%PROGRAM_RUN_INFO")
350 
351  ! terminate the run if the value of the potential is abnormal
352  IF (abnormal_value(e_pot)) &
353  cpabort("Potential energy is an abnormal value (NaN/Inf).")
354 
355  ! Print forces, if requested
356  print_forces = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%FORCES", &
357  extension=".xyz")
358  IF ((print_forces > 0) .AND. calculate_forces) THEN
359  CALL force_env_get(force_env, subsys=subsys)
360  CALL cp_subsys_get(subsys, &
361  core_particles=core_particles, &
362  particles=particles, &
363  shell_particles=shell_particles)
364  ! Variable precision output of the forces
365  CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%NDIGITS", &
366  i_val=ndigits)
367  IF (ASSOCIATED(core_particles) .OR. ASSOCIATED(shell_particles)) THEN
368  CALL write_forces(particles, print_forces, "ATOMIC", ndigits, total_force, &
369  zero_force_core_shell_atom=.true.)
370  grand_total_force(1:3) = total_force(1:3)
371  IF (ASSOCIATED(core_particles)) THEN
372  CALL write_forces(core_particles, print_forces, "CORE", ndigits, total_force, &
373  zero_force_core_shell_atom=.false.)
374  grand_total_force(:) = grand_total_force(:) + total_force(:)
375  END IF
376  IF (ASSOCIATED(shell_particles)) THEN
377  CALL write_forces(shell_particles, print_forces, "SHELL", ndigits, total_force, &
378  zero_force_core_shell_atom=.false., &
379  grand_total_force=grand_total_force)
380  END IF
381  ELSE
382  CALL write_forces(particles, print_forces, "ATOMIC", ndigits, total_force)
383  END IF
384  END IF
385  CALL cp_print_key_finished_output(print_forces, logger, force_env%force_env_section, "PRINT%FORCES")
386 
387  ! Write stress tensor
388  IF (virial%pv_availability) THEN
389  ! If the virial is defined but we are not computing forces let's zero the
390  ! virial for consistency
391  IF (calculate_forces .AND. calculate_stress_tensor) THEN
392  output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
393  extension=".stress_tensor")
394  IF (output_unit > 0) THEN
395  CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%COMPONENTS", &
396  l_val=print_components)
397  IF (print_components) THEN
398  IF ((.NOT. virial%pv_numer) .AND. (force_env%in_use == use_qs_force)) THEN
399  CALL write_stress_tensor_components(virial, output_unit, cell)
400  END IF
401  END IF
402  CALL write_stress_tensor(virial%pv_virial, output_unit, cell, virial%pv_numer)
403  END IF
404  CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
405  "PRINT%STRESS_TENSOR")
406  ELSE
407  CALL zero_virial(virial, reset=.false.)
408  END IF
409  ELSE
410  output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
411  extension=".stress_tensor")
412  IF (output_unit > 0) THEN
413  CALL cp_warn(__location__, "To print the stress tensor switch on the "// &
414  "virial evaluation with the keyword: STRESS_TENSOR")
415  END IF
416  CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
417  "PRINT%STRESS_TENSOR")
418  END IF
419 
420  ! Atomic energy
421  output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
422  extension=".Log")
423  IF (atprop_env%energy) THEN
424  CALL force_env%para_env%sum(atprop_env%atener)
425  CALL force_env_get(force_env, potential_energy=e_pot)
426  IF (output_unit > 0) THEN
427  IF (logger%iter_info%print_level >= low_print_level) THEN
428  CALL cp_subsys_get(subsys=subsys, particles=particles)
429  CALL write_atener(output_unit, particles, atprop_env%atener, "Mulliken Atomic Energies")
430  END IF
431  sum_energy = accurate_sum(atprop_env%atener(:))
432  checksum = abs(e_pot - sum_energy)
433  WRITE (unit=output_unit, fmt="(/,(T2,A,T56,F25.13))") &
434  "Potential energy (Atomic):", sum_energy, &
435  "Potential energy (Total) :", e_pot, &
436  "Difference :", checksum
437  cpassert((checksum < ateps*abs(e_pot)))
438  END IF
439  CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
440  "PRINT%PROGRAM_RUN_INFO")
441  END IF
442 
443  ! Print GRMM interface file
444  print_grrm = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%GRRM", &
445  file_position="REWIND", extension=".rrm")
446  IF (print_grrm > 0) THEN
447  CALL force_env_get(force_env, subsys=subsys)
448  CALL cp_subsys_get(subsys=subsys, particles=particles, &
449  molecule_kinds=molecule_kinds)
450  ! Count the number of fixed atoms
451  nfixed_atoms_total = 0
452  nkind = molecule_kinds%n_els
453  molecule_kind_set => molecule_kinds%els
454  DO ikind = 1, nkind
455  molecule_kind => molecule_kind_set(ikind)
456  CALL get_molecule_kind(molecule_kind, nfixd=nfixed_atoms)
457  nfixed_atoms_total = nfixed_atoms_total + nfixed_atoms
458  END DO
459  !
460  CALL write_grrm(print_grrm, force_env, particles%els, e_pot, fixed_atoms=nfixed_atoms_total)
461  END IF
462  CALL cp_print_key_finished_output(print_grrm, logger, force_env%force_env_section, "PRINT%GRRM")
463 
464  print_scine = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%SCINE", &
465  file_position="REWIND", extension=".scine")
466  IF (print_scine > 0) THEN
467  CALL force_env_get(force_env, subsys=subsys)
468  CALL cp_subsys_get(subsys=subsys, particles=particles)
469  !
470  CALL write_scine(print_scine, force_env, particles%els, e_pot)
471  END IF
472  CALL cp_print_key_finished_output(print_scine, logger, force_env%force_env_section, "PRINT%SCINE")
473 
474  ! Atomic stress
475  output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
476  extension=".Log")
477 
478  IF (atprop_env%stress) THEN
479  CALL force_env%para_env%sum(atprop_env%atstress)
480  ! symmetrize (same as pv_virial)
481  DO i = 1, SIZE(atprop_env%atstress, 3)
482  atprop_env%atstress(:, :, i) = 0.5_dp*(atprop_env%atstress(:, :, i) &
483  + transpose(atprop_env%atstress(:, :, i)))
484  END DO
485  IF (output_unit > 0) THEN
486  IF (logger%iter_info%print_level > low_print_level) THEN
487  DO i = 1, SIZE(atprop_env%atstress, 3)
488  WRITE (unit=output_unit, fmt="(/,T2,I0,T16,A1,2(19X,A1))") i, "X", "Y", "Z"
489  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "X", (atprop_env%atstress(1, j, i), j=1, 3)
490  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Y", (atprop_env%atstress(2, j, i), j=1, 3)
491  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Z", (atprop_env%atstress(3, j, i), j=1, 3)
492  WRITE (unit=output_unit, fmt="(T2,A,F20.13)") "1/3 Trace(Atomic stress tensor):", &
493  (atprop_env%atstress(1, 1, i) + atprop_env%atstress(2, 2, i) + atprop_env%atstress(3, 3, i))/3.0_dp
494  END DO
495  END IF
496  atomic_stress_tensor(:, :) = 0.0_dp
497  DO i = 1, 3
498  atomic_stress_tensor(i, i) = accurate_sum(atprop_env%atstress(i, i, :))
499  DO j = i + 1, 3
500  atomic_stress_tensor(i, j) = accurate_sum(atprop_env%atstress(i, j, :))
501  atomic_stress_tensor(j, i) = atomic_stress_tensor(i, j)
502  END DO
503  END DO
504  WRITE (unit=output_unit, fmt="(/,T2,A,T15,A1,2(19X,A1))") "Atomic", "X", "Y", "Z"
505  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "X", (atomic_stress_tensor(1, i), i=1, 3)
506  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Y", (atomic_stress_tensor(2, i), i=1, 3)
507  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Z", (atomic_stress_tensor(3, i), i=1, 3)
508  WRITE (unit=output_unit, fmt="(T2,A,10X,F20.13)") "1/3 Trace(Atomic stress tensor):", &
509  (atomic_stress_tensor(1, 1) + atomic_stress_tensor(2, 2) + atomic_stress_tensor(3, 3))/3.0_dp
510  sum_stress_tensor = accurate_sum(atomic_stress_tensor(:, :))
511  IF (virial%pv_availability .AND. calculate_forces) THEN
512  WRITE (unit=output_unit, fmt="(/,T2,A,T16,A1,2(19X,A1))") "Total", "X", "Y", "Z"
513  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "X", (virial%pv_virial(1, i), i=1, 3)
514  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Y", (virial%pv_virial(2, i), i=1, 3)
515  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Z", (virial%pv_virial(3, i), i=1, 3)
516  WRITE (unit=output_unit, fmt="(T2,A,10X,F20.13)") "1/3 Trace(Total stress tensor): ", &
517  (virial%pv_virial(1, 1) + virial%pv_virial(2, 2) + virial%pv_virial(3, 3))/3.0_dp
518  sum_pv_virial = sum(virial%pv_virial(:, :))
519  diff_stress_tensor(:, :) = abs(virial%pv_virial(:, :) - atomic_stress_tensor(:, :))
520  WRITE (unit=output_unit, fmt="(/,T2,A,T16,A1,2(19X,A1))") "Diff", "X", "Y", "Z"
521  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "X", (diff_stress_tensor(1, i), i=1, 3)
522  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Y", (diff_stress_tensor(2, i), i=1, 3)
523  WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Z", (diff_stress_tensor(3, i), i=1, 3)
524  WRITE (unit=output_unit, fmt="(T2,A,10X,F20.13)") "1/3 Trace(Diff) : ", &
525  (diff_stress_tensor(1, 1) + diff_stress_tensor(2, 2) + diff_stress_tensor(3, 3))/3.0_dp
526  checksum = accurate_sum(diff_stress_tensor(:, :))
527  WRITE (unit=output_unit, fmt="(/,(T2,A,11X,F25.13))") &
528  "Checksum stress (Atomic) :", sum_stress_tensor, &
529  "Checksum stress (Total) :", sum_pv_virial, &
530  "Difference :", checksum
531  cpassert(checksum < ateps)
532  END IF
533  END IF
534  CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
535  "PRINT%PROGRAM_RUN_INFO")
536  END IF
537 
538  END SUBROUTINE force_env_calc_energy_force
539 
540 ! **************************************************************************************************
541 !> \brief Evaluates the stress tensor and pressure numerically
542 !> \param force_env ...
543 !> \param dx ...
544 !> \par History
545 !> 10.2005 created [JCS]
546 !> 05.2009 Teodoro Laino [tlaino] - rewriting for general force_env
547 !>
548 !> \author JCS
549 ! **************************************************************************************************
550  SUBROUTINE force_env_calc_num_pressure(force_env, dx)
551 
552  TYPE(force_env_type), POINTER :: force_env
553  REAL(kind=dp), INTENT(IN), OPTIONAL :: dx
554 
555  REAL(kind=dp), PARAMETER :: default_dx = 0.001_dp
556 
557  INTEGER :: i, ip, iq, j, k, natom, ncore, nshell, &
558  output_unit, symmetry_id
559  REAL(kind=dp) :: dx_w
560  REAL(kind=dp), DIMENSION(2) :: numer_energy
561  REAL(kind=dp), DIMENSION(3) :: s
562  REAL(kind=dp), DIMENSION(3, 3) :: numer_stress
563  REAL(kind=dp), DIMENSION(:, :), POINTER :: ref_pos_atom, ref_pos_core, ref_pos_shell
564  TYPE(cell_type), POINTER :: cell, cell_local
565  TYPE(cp_logger_type), POINTER :: logger
566  TYPE(cp_subsys_type), POINTER :: subsys
567  TYPE(global_environment_type), POINTER :: globenv
568  TYPE(particle_list_type), POINTER :: core_particles, particles, &
569  shell_particles
570  TYPE(virial_type), POINTER :: virial
571 
572  NULLIFY (cell_local)
573  NULLIFY (core_particles)
574  NULLIFY (particles)
575  NULLIFY (shell_particles)
576  NULLIFY (ref_pos_atom)
577  NULLIFY (ref_pos_core)
578  NULLIFY (ref_pos_shell)
579  natom = 0
580  ncore = 0
581  nshell = 0
582  numer_stress = 0.0_dp
583 
584  logger => cp_get_default_logger()
585 
586  dx_w = default_dx
587  IF (PRESENT(dx)) dx_w = dx
588  CALL force_env_get(force_env, subsys=subsys, globenv=globenv)
589  CALL cp_subsys_get(subsys, &
590  core_particles=core_particles, &
591  particles=particles, &
592  shell_particles=shell_particles, &
593  virial=virial)
594  output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
595  extension=".stress_tensor")
596  IF (output_unit > 0) THEN
597  WRITE (output_unit, '(/A,A/)') ' **************************** ', &
598  'NUMERICAL STRESS ********************************'
599  END IF
600 
601  ! Save all original particle positions
602  natom = particles%n_els
603  ALLOCATE (ref_pos_atom(natom, 3))
604  DO i = 1, natom
605  ref_pos_atom(i, :) = particles%els(i)%r
606  END DO
607  IF (ASSOCIATED(core_particles)) THEN
608  ncore = core_particles%n_els
609  ALLOCATE (ref_pos_core(ncore, 3))
610  DO i = 1, ncore
611  ref_pos_core(i, :) = core_particles%els(i)%r
612  END DO
613  END IF
614  IF (ASSOCIATED(shell_particles)) THEN
615  nshell = shell_particles%n_els
616  ALLOCATE (ref_pos_shell(nshell, 3))
617  DO i = 1, nshell
618  ref_pos_shell(i, :) = shell_particles%els(i)%r
619  END DO
620  END IF
621  CALL force_env_get(force_env, cell=cell)
622  ! Save cell symmetry (distorted cell has no symmetry)
623  symmetry_id = cell%symmetry_id
624  cell%symmetry_id = cell_sym_triclinic
625  !
626  CALL cell_create(cell_local)
627  CALL cell_clone(cell, cell_local)
628  ! First change box
629  DO ip = 1, 3
630  DO iq = 1, 3
631  IF (virial%pv_diagonal .AND. (ip /= iq)) cycle
632  DO k = 1, 2
633  cell%hmat(ip, iq) = cell_local%hmat(ip, iq) - (-1.0_dp)**k*dx_w
634  CALL init_cell(cell)
635  ! Scale positions
636  DO i = 1, natom
637  CALL real_to_scaled(s, ref_pos_atom(i, 1:3), cell_local)
638  CALL scaled_to_real(particles%els(i)%r, s, cell)
639  END DO
640  DO i = 1, ncore
641  CALL real_to_scaled(s, ref_pos_core(i, 1:3), cell_local)
642  CALL scaled_to_real(core_particles%els(i)%r, s, cell)
643  END DO
644  DO i = 1, nshell
645  CALL real_to_scaled(s, ref_pos_shell(i, 1:3), cell_local)
646  CALL scaled_to_real(shell_particles%els(i)%r, s, cell)
647  END DO
648  ! Compute energies
649  CALL force_env_calc_energy_force(force_env, &
650  calc_force=.false., &
651  consistent_energies=.true., &
652  calc_stress_tensor=.false.)
653  CALL force_env_get(force_env, potential_energy=numer_energy(k))
654  ! Reset cell
655  cell%hmat(ip, iq) = cell_local%hmat(ip, iq)
656  END DO
657  CALL init_cell(cell)
658  numer_stress(ip, iq) = 0.5_dp*(numer_energy(1) - numer_energy(2))/dx_w
659  IF (output_unit > 0) THEN
660  IF (globenv%run_type_id == debug_run) THEN
661  WRITE (unit=output_unit, fmt="(/,T2,A,T19,A,F7.4,A,T44,A,F7.4,A,T69,A)") &
662  "DEBUG|", "E("//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
663  "E("//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
664  "f(numerical)"
665  WRITE (unit=output_unit, fmt="(T2,A,2(1X,F24.8),1X,F22.8)") &
666  "DEBUG|", numer_energy(1:2), numer_stress(ip, iq)
667  ELSE
668  WRITE (unit=output_unit, fmt="(/,T7,A,F7.4,A,T27,A,F7.4,A,T49,A)") &
669  "E("//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
670  "E("//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
671  "f(numerical)"
672  WRITE (unit=output_unit, fmt="(3(1X,F19.8))") &
673  numer_energy(1:2), numer_stress(ip, iq)
674  END IF
675  END IF
676  END DO
677  END DO
678 
679  ! Reset positions and rebuild original environment
680  cell%symmetry_id = symmetry_id
681  CALL init_cell(cell)
682  DO i = 1, natom
683  particles%els(i)%r = ref_pos_atom(i, :)
684  END DO
685  DO i = 1, ncore
686  core_particles%els(i)%r = ref_pos_core(i, :)
687  END DO
688  DO i = 1, nshell
689  shell_particles%els(i)%r = ref_pos_shell(i, :)
690  END DO
691  CALL force_env_calc_energy_force(force_env, &
692  calc_force=.false., &
693  consistent_energies=.true., &
694  calc_stress_tensor=.false.)
695 
696  ! Computing pv_test
697  virial%pv_virial = 0.0_dp
698  DO i = 1, 3
699  DO j = 1, 3
700  DO k = 1, 3
701  virial%pv_virial(i, j) = virial%pv_virial(i, j) - &
702  0.5_dp*(numer_stress(i, k)*cell_local%hmat(j, k) + &
703  numer_stress(j, k)*cell_local%hmat(i, k))
704  END DO
705  END DO
706  END DO
707 
708  IF (output_unit > 0) THEN
709  IF (globenv%run_type_id == debug_run) THEN
710  CALL write_stress_tensor(virial%pv_virial, output_unit, cell, virial%pv_numer)
711  END IF
712  WRITE (output_unit, '(/,A,/)') ' **************************** '// &
713  'NUMERICAL STRESS END *****************************'
714  END IF
715 
716  CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
717  "PRINT%STRESS_TENSOR")
718 
719  ! Release storage
720  IF (ASSOCIATED(ref_pos_atom)) THEN
721  DEALLOCATE (ref_pos_atom)
722  END IF
723  IF (ASSOCIATED(ref_pos_core)) THEN
724  DEALLOCATE (ref_pos_core)
725  END IF
726  IF (ASSOCIATED(ref_pos_shell)) THEN
727  DEALLOCATE (ref_pos_shell)
728  END IF
729  IF (ASSOCIATED(cell_local)) CALL cell_release(cell_local)
730 
731  END SUBROUTINE force_env_calc_num_pressure
732 
733 ! **************************************************************************************************
734 !> \brief creates and initializes a force environment
735 !> \param force_env the force env to create
736 !> \param root_section ...
737 !> \param para_env ...
738 !> \param globenv ...
739 !> \param fist_env , qs_env: exactly one of these should be
740 !> associated, the one that is active
741 !> \param qs_env ...
742 !> \param meta_env ...
743 !> \param sub_force_env ...
744 !> \param qmmm_env ...
745 !> \param qmmmx_env ...
746 !> \param eip_env ...
747 !> \param pwdft_env ...
748 !> \param force_env_section ...
749 !> \param mixed_env ...
750 !> \param embed_env ...
751 !> \param nnp_env ...
752 !> \par History
753 !> 04.2003 created [fawzi]
754 !> \author fawzi
755 ! **************************************************************************************************
756  SUBROUTINE force_env_create(force_env, root_section, para_env, globenv, fist_env, &
757  qs_env, meta_env, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, force_env_section, &
758  mixed_env, embed_env, nnp_env)
759 
760  TYPE(force_env_type), POINTER :: force_env
761  TYPE(section_vals_type), POINTER :: root_section
762  TYPE(mp_para_env_type), POINTER :: para_env
763  TYPE(global_environment_type), POINTER :: globenv
764  TYPE(fist_environment_type), OPTIONAL, POINTER :: fist_env
765  TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
766  TYPE(meta_env_type), OPTIONAL, POINTER :: meta_env
767  TYPE(force_env_p_type), DIMENSION(:), OPTIONAL, &
768  POINTER :: sub_force_env
769  TYPE(qmmm_env_type), OPTIONAL, POINTER :: qmmm_env
770  TYPE(qmmmx_env_type), OPTIONAL, POINTER :: qmmmx_env
771  TYPE(eip_environment_type), OPTIONAL, POINTER :: eip_env
772  TYPE(pwdft_environment_type), OPTIONAL, POINTER :: pwdft_env
773  TYPE(section_vals_type), POINTER :: force_env_section
774  TYPE(mixed_environment_type), OPTIONAL, POINTER :: mixed_env
775  TYPE(embed_env_type), OPTIONAL, POINTER :: embed_env
776  TYPE(nnp_type), OPTIONAL, POINTER :: nnp_env
777 
778  ALLOCATE (force_env)
779  NULLIFY (force_env%fist_env, force_env%qs_env, &
780  force_env%para_env, force_env%globenv, &
781  force_env%meta_env, force_env%sub_force_env, &
782  force_env%qmmm_env, force_env%qmmmx_env, force_env%fp_env, &
783  force_env%force_env_section, force_env%eip_env, force_env%mixed_env, &
784  force_env%embed_env, force_env%pwdft_env, force_env%nnp_env, &
785  force_env%root_section)
786  last_force_env_id = last_force_env_id + 1
787  force_env%ref_count = 1
788  force_env%in_use = 0
789  force_env%additional_potential = 0.0_dp
790 
791  force_env%globenv => globenv
792  CALL globenv_retain(force_env%globenv)
793 
794  force_env%root_section => root_section
795  CALL section_vals_retain(root_section)
796 
797  force_env%para_env => para_env
798  CALL force_env%para_env%retain()
799 
800  CALL section_vals_retain(force_env_section)
801  force_env%force_env_section => force_env_section
802 
803  IF (PRESENT(fist_env)) THEN
804  cpassert(ASSOCIATED(fist_env))
805  cpassert(force_env%in_use == 0)
806  force_env%in_use = use_fist_force
807  force_env%fist_env => fist_env
808  END IF
809  IF (PRESENT(eip_env)) THEN
810  cpassert(ASSOCIATED(eip_env))
811  cpassert(force_env%in_use == 0)
812  force_env%in_use = use_eip_force
813  force_env%eip_env => eip_env
814  END IF
815  IF (PRESENT(pwdft_env)) THEN
816  cpassert(ASSOCIATED(pwdft_env))
817  cpassert(force_env%in_use == 0)
818  force_env%in_use = use_pwdft_force
819  force_env%pwdft_env => pwdft_env
820  END IF
821  IF (PRESENT(qs_env)) THEN
822  cpassert(ASSOCIATED(qs_env))
823  cpassert(force_env%in_use == 0)
824  force_env%in_use = use_qs_force
825  force_env%qs_env => qs_env
826  END IF
827  IF (PRESENT(qmmm_env)) THEN
828  cpassert(ASSOCIATED(qmmm_env))
829  cpassert(force_env%in_use == 0)
830  force_env%in_use = use_qmmm
831  force_env%qmmm_env => qmmm_env
832  END IF
833  IF (PRESENT(qmmmx_env)) THEN
834  cpassert(ASSOCIATED(qmmmx_env))
835  cpassert(force_env%in_use == 0)
836  force_env%in_use = use_qmmmx
837  force_env%qmmmx_env => qmmmx_env
838  END IF
839  IF (PRESENT(mixed_env)) THEN
840  cpassert(ASSOCIATED(mixed_env))
841  cpassert(force_env%in_use == 0)
842  force_env%in_use = use_mixed_force
843  force_env%mixed_env => mixed_env
844  END IF
845  IF (PRESENT(embed_env)) THEN
846  cpassert(ASSOCIATED(embed_env))
847  cpassert(force_env%in_use == 0)
848  force_env%in_use = use_embed
849  force_env%embed_env => embed_env
850  END IF
851  IF (PRESENT(nnp_env)) THEN
852  cpassert(ASSOCIATED(nnp_env))
853  cpassert(force_env%in_use == 0)
854  force_env%in_use = use_nnp_force
855  force_env%nnp_env => nnp_env
856  END IF
857  cpassert(force_env%in_use /= 0)
858 
859  IF (PRESENT(sub_force_env)) THEN
860  force_env%sub_force_env => sub_force_env
861  END IF
862 
863  IF (PRESENT(meta_env)) THEN
864  force_env%meta_env => meta_env
865  ELSE
866  NULLIFY (force_env%meta_env)
867  END IF
868 
869  END SUBROUTINE force_env_create
870 
871 ! **************************************************************************************************
872 !> \brief ****f* force_env_methods/mixed_energy_forces [1.0]
873 !>
874 !> Computes energy and forces for a mixed force_env type
875 !> \param force_env the force_env that holds the mixed_env type
876 !> \param calculate_forces decides if forces should be calculated
877 !> \par History
878 !> 11.06 created [fschiff]
879 !> 04.07 generalization to an illimited number of force_eval [tlaino]
880 !> 04.07 further generalization to force_eval with different geometrical
881 !> structures [tlaino]
882 !> 04.08 reorganizing the genmix structure (collecting common code)
883 !> 01.16 added CDFT [Nico Holmberg]
884 !> 08.17 added DFT embedding [Vladimir Rybkin]
885 !> \author Florian Schiffmann
886 ! **************************************************************************************************
887  SUBROUTINE mixed_energy_forces(force_env, calculate_forces)
888 
889  TYPE(force_env_type), POINTER :: force_env
890  LOGICAL, INTENT(IN) :: calculate_forces
891 
892  CHARACTER(LEN=default_path_length) :: coupling_function
893  CHARACTER(LEN=default_string_length) :: def_error, description, this_error
894  INTEGER :: iforce_eval, iparticle, istate(2), &
895  jparticle, mixing_type, my_group, &
896  natom, nforce_eval, source, unit_nr
897  INTEGER, DIMENSION(:), POINTER :: glob_natoms, itmplist, map_index
898  LOGICAL :: dip_exists
899  REAL(kind=dp) :: coupling_parameter, dedf, der_1, der_2, &
900  dx, energy, err, lambda, lerr, &
901  restraint_strength, restraint_target, &
902  sd
903  REAL(kind=dp), DIMENSION(3) :: dip_mix
904  REAL(kind=dp), DIMENSION(:), POINTER :: energies
905  TYPE(cell_type), POINTER :: cell_mix
906  TYPE(cp_logger_type), POINTER :: logger, my_logger
907  TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
908  TYPE(cp_result_type), POINTER :: loc_results, results_mix
909  TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
910  TYPE(cp_subsys_type), POINTER :: subsys_mix
911  TYPE(mixed_energy_type), POINTER :: mixed_energy
912  TYPE(mixed_force_type), DIMENSION(:), POINTER :: global_forces
913  TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
914  TYPE(particle_list_type), POINTER :: particles_mix
915  TYPE(section_vals_type), POINTER :: force_env_section, gen_section, &
916  mapping_section, mixed_section, &
917  root_section
918  TYPE(virial_p_type), DIMENSION(:), POINTER :: virials
919  TYPE(virial_type), POINTER :: loc_virial, virial_mix
920 
921  logger => cp_get_default_logger()
922  cpassert(ASSOCIATED(force_env))
923  ! Get infos about the mixed subsys
924  CALL force_env_get(force_env=force_env, &
925  subsys=subsys_mix, &
926  force_env_section=force_env_section, &
927  root_section=root_section, &
928  cell=cell_mix)
929  CALL cp_subsys_get(subsys=subsys_mix, &
930  particles=particles_mix, &
931  virial=virial_mix, &
932  results=results_mix)
933  NULLIFY (map_index, glob_natoms, global_forces, itmplist)
934 
935  nforce_eval = SIZE(force_env%sub_force_env)
936  mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
937  mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
938  ! Global Info
939  ALLOCATE (subsystems(nforce_eval))
940  ALLOCATE (particles(nforce_eval))
941  ! Local Info to sync
942  ALLOCATE (global_forces(nforce_eval))
943  ALLOCATE (energies(nforce_eval))
944  ALLOCATE (glob_natoms(nforce_eval))
945  ALLOCATE (virials(nforce_eval))
946  ALLOCATE (results(nforce_eval))
947  energies = 0.0_dp
948  glob_natoms = 0
949  ! Check if mixed CDFT calculation is requested and initialize
950  CALL mixed_cdft_init(force_env, calculate_forces)
951 
952  !
953  IF (.NOT. force_env%mixed_env%do_mixed_cdft) THEN
954  DO iforce_eval = 1, nforce_eval
955  NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
956  NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
957  ALLOCATE (virials(iforce_eval)%virial)
958  CALL cp_result_create(results(iforce_eval)%results)
959  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
960  ! From this point on the error is the sub_error
961  my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
962  my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
963  ! Copy iterations info (they are updated only in the main mixed_env)
964  CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
965  CALL cp_add_default_logger(my_logger)
966 
967  ! Get all available subsys
968  CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
969  subsys=subsystems(iforce_eval)%subsys)
970 
971  ! all force_env share the same cell
972  CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
973 
974  ! Get available particles
975  CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
976  particles=particles(iforce_eval)%list)
977 
978  ! Get Mapping index array
979  natom = SIZE(particles(iforce_eval)%list%els)
980 
981  CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
982  map_index)
983 
984  ! Mapping particles from iforce_eval environment to the mixed env
985  DO iparticle = 1, natom
986  jparticle = map_index(iparticle)
987  particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
988  END DO
989 
990  ! Calculate energy and forces for each sub_force_env
991  CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
992  calc_force=calculate_forces, &
993  skip_external_control=.true.)
994 
995  ! Only the rank 0 process collect info for each computation
996  IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
997  CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
998  potential_energy=energy)
999  CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1000  virial=loc_virial, results=loc_results)
1001  energies(iforce_eval) = energy
1002  glob_natoms(iforce_eval) = natom
1003  virials(iforce_eval)%virial = loc_virial
1004  CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1005  END IF
1006  ! Deallocate map_index array
1007  IF (ASSOCIATED(map_index)) THEN
1008  DEALLOCATE (map_index)
1009  END IF
1010  CALL cp_rm_default_logger()
1011  END DO
1012  ELSE
1013  CALL mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
1014  glob_natoms, virials, results)
1015  END IF
1016  ! Handling Parallel execution
1017  CALL force_env%para_env%sync()
1018  ! Post CDFT operations
1019  CALL mixed_cdft_post_energy_forces(force_env)
1020  ! Let's transfer energy, natom, forces, virials
1021  CALL force_env%para_env%sum(energies)
1022  CALL force_env%para_env%sum(glob_natoms)
1023  ! Transfer forces
1024  DO iforce_eval = 1, nforce_eval
1025  ALLOCATE (global_forces(iforce_eval)%forces(3, glob_natoms(iforce_eval)))
1026  global_forces(iforce_eval)%forces = 0.0_dp
1027  IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN
1028  IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1029  ! Forces
1030  DO iparticle = 1, glob_natoms(iforce_eval)
1031  global_forces(iforce_eval)%forces(:, iparticle) = &
1032  particles(iforce_eval)%list%els(iparticle)%f
1033  END DO
1034  END IF
1035  END IF
1036  CALL force_env%para_env%sum(global_forces(iforce_eval)%forces)
1037  !Transfer only the relevant part of the virial..
1038  CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_total)
1039  CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_kinetic)
1040  CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_virial)
1041  CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_xc)
1042  CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_fock_4c)
1043  CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_constraint)
1044  !Transfer results
1045  source = 0
1046  IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN
1047  IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) &
1048  source = force_env%para_env%mepos
1049  END IF
1050  CALL force_env%para_env%sum(source)
1051  CALL cp_results_mp_bcast(results(iforce_eval)%results, source, force_env%para_env)
1052  END DO
1053 
1054  force_env%mixed_env%energies = energies
1055  ! Start combining the different sub_force_env
1056  CALL get_mixed_env(mixed_env=force_env%mixed_env, &
1057  mixed_energy=mixed_energy)
1058 
1059  !NB: do this for all MIXING_TYPE values, since some need it (e.g. linear mixing
1060  !NB if the first system has fewer atoms than the second)
1061  DO iparticle = 1, SIZE(particles_mix%els)
1062  particles_mix%els(iparticle)%f(:) = 0.0_dp
1063  END DO
1064 
1065  CALL section_vals_val_get(mixed_section, "MIXING_TYPE", i_val=mixing_type)
1066  SELECT CASE (mixing_type)
1067  CASE (mix_linear_combination)
1068  ! Support offered only 2 force_eval
1069  cpassert(nforce_eval == 2)
1070  CALL section_vals_val_get(mixed_section, "LINEAR%LAMBDA", r_val=lambda)
1071  mixed_energy%pot = lambda*energies(1) + (1.0_dp - lambda)*energies(2)
1072  ! General Mapping of forces...
1073  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1074  lambda, 1, nforce_eval, map_index, mapping_section, .true.)
1075  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1076  (1.0_dp - lambda), 2, nforce_eval, map_index, mapping_section, .false.)
1077  CASE (mix_minimum)
1078  ! Support offered only 2 force_eval
1079  cpassert(nforce_eval == 2)
1080  IF (energies(1) < energies(2)) THEN
1081  mixed_energy%pot = energies(1)
1082  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1083  1.0_dp, 1, nforce_eval, map_index, mapping_section, .true.)
1084  ELSE
1085  mixed_energy%pot = energies(2)
1086  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1087  1.0_dp, 2, nforce_eval, map_index, mapping_section, .true.)
1088  END IF
1089  CASE (mix_coupled)
1090  ! Support offered only 2 force_eval
1091  cpassert(nforce_eval == 2)
1092  CALL section_vals_val_get(mixed_section, "COUPLING%COUPLING_PARAMETER", &
1093  r_val=coupling_parameter)
1094  sd = sqrt((energies(1) - energies(2))**2 + 4.0_dp*coupling_parameter**2)
1095  der_1 = (1.0_dp - (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1096  der_2 = (1.0_dp + (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1097  mixed_energy%pot = (energies(1) + energies(2) - sd)/2.0_dp
1098  ! General Mapping of forces...
1099  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1100  der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1101  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1102  der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1103  CASE (mix_restrained)
1104  ! Support offered only 2 force_eval
1105  cpassert(nforce_eval == 2)
1106  CALL section_vals_val_get(mixed_section, "RESTRAINT%RESTRAINT_TARGET", &
1107  r_val=restraint_target)
1108  CALL section_vals_val_get(mixed_section, "RESTRAINT%RESTRAINT_STRENGTH", &
1109  r_val=restraint_strength)
1110  mixed_energy%pot = energies(1) + restraint_strength*(energies(1) - energies(2) - restraint_target)**2
1111  der_2 = -2.0_dp*restraint_strength*(energies(1) - energies(2) - restraint_target)
1112  der_1 = 1.0_dp - der_2
1113  ! General Mapping of forces...
1114  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1115  der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1116  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1117  der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1118  CASE (mix_generic)
1119  ! Support any number of force_eval sections
1120  gen_section => section_vals_get_subs_vals(mixed_section, "GENERIC")
1121  CALL get_generic_info(gen_section, "MIXING_FUNCTION", coupling_function, force_env%mixed_env%par, &
1122  force_env%mixed_env%val, energies)
1123  CALL initf(1)
1124  CALL parsef(1, trim(coupling_function), force_env%mixed_env%par)
1125  ! Now the hardest part.. map energy with corresponding force_eval
1126  mixed_energy%pot = evalf(1, force_env%mixed_env%val)
1127  cpassert(evalerrtype <= 0)
1128  CALL zero_virial(virial_mix, reset=.false.)
1129  CALL cp_results_erase(results_mix)
1130  DO iforce_eval = 1, nforce_eval
1131  CALL section_vals_val_get(gen_section, "DX", r_val=dx)
1132  CALL section_vals_val_get(gen_section, "ERROR_LIMIT", r_val=lerr)
1133  dedf = evalfd(1, iforce_eval, force_env%mixed_env%val, dx, err)
1134  IF (abs(err) > lerr) THEN
1135  WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
1136  WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
1137  CALL compress(this_error, .true.)
1138  CALL compress(def_error, .true.)
1139  CALL cp_warn(__location__, &
1140  'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
1141  ' Error '//trim(this_error)//' in computing numerical derivatives larger then'// &
1142  trim(def_error)//' .')
1143  END IF
1144  ! General Mapping of forces...
1145  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1146  dedf, iforce_eval, nforce_eval, map_index, mapping_section, .false.)
1147  force_env%mixed_env%val(iforce_eval) = energies(iforce_eval)
1148  END DO
1149  ! Let's store the needed information..
1150  force_env%mixed_env%dx = dx
1151  force_env%mixed_env%lerr = lerr
1152  force_env%mixed_env%coupling_function = trim(coupling_function)
1153  CALL finalizef()
1154  CASE (mix_cdft)
1155  ! Supports any number of force_evals for calculation of CDFT properties, but forces only from two
1156  CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LAMBDA", r_val=lambda)
1157  ! Get the states which determine the forces
1158  CALL section_vals_val_get(mixed_section, "MIXED_CDFT%FORCE_STATES", i_vals=itmplist)
1159  IF (SIZE(itmplist) /= 2) &
1160  CALL cp_abort(__location__, &
1161  "Keyword FORCE_STATES takes exactly two input values.")
1162  IF (any(itmplist .LT. 0)) &
1163  cpabort("Invalid force_eval index.")
1164  istate = itmplist
1165  IF (istate(1) > nforce_eval .OR. istate(2) > nforce_eval) &
1166  cpabort("Invalid force_eval index.")
1167  mixed_energy%pot = lambda*energies(istate(1)) + (1.0_dp - lambda)*energies(istate(2))
1168  ! General Mapping of forces...
1169  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1170  lambda, istate(1), nforce_eval, map_index, mapping_section, .true.)
1171  CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1172  (1.0_dp - lambda), istate(2), nforce_eval, map_index, mapping_section, .false.)
1173  CASE DEFAULT
1174  cpabort("")
1175  END SELECT
1176  !Simply deallocate and loose the pointer references..
1177  DO iforce_eval = 1, nforce_eval
1178  DEALLOCATE (global_forces(iforce_eval)%forces)
1179  IF (ASSOCIATED(virials(iforce_eval)%virial)) DEALLOCATE (virials(iforce_eval)%virial)
1180  CALL cp_result_release(results(iforce_eval)%results)
1181  END DO
1182  DEALLOCATE (global_forces)
1183  DEALLOCATE (subsystems)
1184  DEALLOCATE (particles)
1185  DEALLOCATE (energies)
1186  DEALLOCATE (glob_natoms)
1187  DEALLOCATE (virials)
1188  DEALLOCATE (results)
1189  ! Print Section
1190  unit_nr = cp_print_key_unit_nr(logger, mixed_section, "PRINT%DIPOLE", &
1191  extension=".data", middle_name="MIXED_DIPOLE", log_filename=.false.)
1192  IF (unit_nr > 0) THEN
1193  description = '[DIPOLE]'
1194  dip_exists = test_for_result(results=results_mix, description=description)
1195  IF (dip_exists) THEN
1196  CALL get_results(results=results_mix, description=description, values=dip_mix)
1197  WRITE (unit_nr, '(/,1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE ( A.U.)|", dip_mix
1198  WRITE (unit_nr, '( 1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE (Debye)|", dip_mix*debye
1199  ELSE
1200  WRITE (unit_nr, *) "NO FORCE_EVAL section calculated the dipole"
1201  END IF
1202  END IF
1203  CALL cp_print_key_finished_output(unit_nr, logger, mixed_section, "PRINT%DIPOLE")
1204  END SUBROUTINE mixed_energy_forces
1205 
1206 ! **************************************************************************************************
1207 !> \brief Driver routine for mixed CDFT energy and force calculations
1208 !> \param force_env the force_env that holds the mixed_env
1209 !> \param calculate_forces if forces should be calculated
1210 !> \param particles system particles
1211 !> \param energies the energies of the CDFT states
1212 !> \param glob_natoms the total number of particles
1213 !> \param virials the virials stored in subsys
1214 !> \param results results stored in subsys
1215 !> \par History
1216 !> 01.17 created [Nico Holmberg]
1217 !> \author Nico Holmberg
1218 ! **************************************************************************************************
1219  SUBROUTINE mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
1220  glob_natoms, virials, results)
1221  TYPE(force_env_type), POINTER :: force_env
1222  LOGICAL, INTENT(IN) :: calculate_forces
1223  TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1224  REAL(kind=dp), DIMENSION(:), POINTER :: energies
1225  INTEGER, DIMENSION(:), POINTER :: glob_natoms
1226  TYPE(virial_p_type), DIMENSION(:), POINTER :: virials
1227  TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1228 
1229  INTEGER :: iforce_eval, iparticle, jparticle, &
1230  my_group, natom, nforce_eval
1231  INTEGER, DIMENSION(:), POINTER :: map_index
1232  REAL(kind=dp) :: energy
1233  TYPE(cell_type), POINTER :: cell_mix
1234  TYPE(cp_logger_type), POINTER :: logger, my_logger
1235  TYPE(cp_result_type), POINTER :: loc_results, results_mix
1236  TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1237  TYPE(cp_subsys_type), POINTER :: subsys_mix
1238  TYPE(particle_list_type), POINTER :: particles_mix
1239  TYPE(section_vals_type), POINTER :: force_env_section, mapping_section, &
1240  mixed_section, root_section
1241  TYPE(virial_type), POINTER :: loc_virial, virial_mix
1242 
1243  logger => cp_get_default_logger()
1244  cpassert(ASSOCIATED(force_env))
1245  ! Get infos about the mixed subsys
1246  CALL force_env_get(force_env=force_env, &
1247  subsys=subsys_mix, &
1248  force_env_section=force_env_section, &
1249  root_section=root_section, &
1250  cell=cell_mix)
1251  CALL cp_subsys_get(subsys=subsys_mix, &
1252  particles=particles_mix, &
1253  virial=virial_mix, &
1254  results=results_mix)
1255  NULLIFY (map_index)
1256  nforce_eval = SIZE(force_env%sub_force_env)
1257  mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
1258  mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
1259  ALLOCATE (subsystems(nforce_eval))
1260  DO iforce_eval = 1, nforce_eval
1261  NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1262  NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
1263  ALLOCATE (virials(iforce_eval)%virial)
1264  CALL cp_result_create(results(iforce_eval)%results)
1265  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1266  ! Get all available subsys
1267  CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1268  subsys=subsystems(iforce_eval)%subsys)
1269 
1270  ! all force_env share the same cell
1271  CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
1272 
1273  ! Get available particles
1274  CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1275  particles=particles(iforce_eval)%list)
1276 
1277  ! Get Mapping index array
1278  natom = SIZE(particles(iforce_eval)%list%els)
1279  ! Serial mode need to deallocate first
1280  IF (ASSOCIATED(map_index)) &
1281  DEALLOCATE (map_index)
1282  CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1283  map_index)
1284 
1285  ! Mapping particles from iforce_eval environment to the mixed env
1286  DO iparticle = 1, natom
1287  jparticle = map_index(iparticle)
1288  particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
1289  END DO
1290  ! Mixed CDFT + QMMM: Need to translate now
1291  IF (force_env%mixed_env%do_mixed_qmmm_cdft) &
1292  CALL apply_qmmm_translate(force_env%sub_force_env(iforce_eval)%force_env%qmmm_env)
1293  END DO
1294  ! For mixed CDFT calculations parallelized over CDFT states
1295  ! build weight and gradient on all processors before splitting into groups and
1296  ! starting energy calculation
1297  CALL mixed_cdft_build_weight(force_env, calculate_forces)
1298  DO iforce_eval = 1, nforce_eval
1299  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1300  ! From this point on the error is the sub_error
1301  IF (force_env%mixed_env%cdft_control%run_type == mixed_cdft_serial .AND. iforce_eval .GE. 2) THEN
1302  my_logger => force_env%mixed_env%cdft_control%sub_logger(iforce_eval - 1)%p
1303  ELSE
1304  my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
1305  my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
1306  END IF
1307  ! Copy iterations info (they are updated only in the main mixed_env)
1308  CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1309  CALL cp_add_default_logger(my_logger)
1310  ! Serial CDFT calculation: transfer weight/gradient
1311  CALL mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
1312  ! Calculate energy and forces for each sub_force_env
1313  CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1314  calc_force=calculate_forces, &
1315  skip_external_control=.true.)
1316  ! Only the rank 0 process collect info for each computation
1317  IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1318  CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1319  potential_energy=energy)
1320  CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1321  virial=loc_virial, results=loc_results)
1322  energies(iforce_eval) = energy
1323  glob_natoms(iforce_eval) = natom
1324  virials(iforce_eval)%virial = loc_virial
1325  CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1326  END IF
1327  ! Deallocate map_index array
1328  IF (ASSOCIATED(map_index)) THEN
1329  DEALLOCATE (map_index)
1330  END IF
1331  CALL cp_rm_default_logger()
1332  END DO
1333  DEALLOCATE (subsystems)
1334 
1335  END SUBROUTINE mixed_cdft_energy_forces
1336 
1337 ! **************************************************************************************************
1338 !> \brief Perform additional tasks for mixed CDFT calculations after solving the electronic structure
1339 !> of both CDFT states
1340 !> \param force_env the force_env that holds the CDFT states
1341 !> \par History
1342 !> 01.17 created [Nico Holmberg]
1343 !> \author Nico Holmberg
1344 ! **************************************************************************************************
1345  SUBROUTINE mixed_cdft_post_energy_forces(force_env)
1346  TYPE(force_env_type), POINTER :: force_env
1347 
1348  INTEGER :: iforce_eval, nforce_eval, nvar
1349  TYPE(dft_control_type), POINTER :: dft_control
1350  TYPE(qs_environment_type), POINTER :: qs_env
1351 
1352  cpassert(ASSOCIATED(force_env))
1353  NULLIFY (qs_env, dft_control)
1354  IF (force_env%mixed_env%do_mixed_cdft) THEN
1355  nforce_eval = SIZE(force_env%sub_force_env)
1356  nvar = force_env%mixed_env%cdft_control%nconstraint
1357  ! Transfer cdft strengths for writing restart
1358  IF (.NOT. ASSOCIATED(force_env%mixed_env%strength)) &
1359  ALLOCATE (force_env%mixed_env%strength(nforce_eval, nvar))
1360  force_env%mixed_env%strength = 0.0_dp
1361  DO iforce_eval = 1, nforce_eval
1362  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1363  IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1364  qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
1365  ELSE
1366  CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1367  END IF
1368  CALL get_qs_env(qs_env, dft_control=dft_control)
1369  IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) &
1370  force_env%mixed_env%strength(iforce_eval, :) = dft_control%qs_control%cdft_control%strength(:)
1371  END DO
1372  CALL force_env%para_env%sum(force_env%mixed_env%strength)
1373  ! Mixed CDFT: calculate ET coupling
1374  IF (force_env%mixed_env%do_mixed_et) THEN
1375  IF (modulo(force_env%mixed_env%cdft_control%sim_step, force_env%mixed_env%et_freq) == 0) &
1376  CALL mixed_cdft_calculate_coupling(force_env)
1377  END IF
1378  END IF
1379 
1380  END SUBROUTINE mixed_cdft_post_energy_forces
1381 
1382 ! **************************************************************************************************
1383 !> \brief Computes the total energy for an embedded calculation
1384 !> \param force_env ...
1385 !> \author Vladimir Rybkin
1386 ! **************************************************************************************************
1387  SUBROUTINE embed_energy(force_env)
1388 
1389  TYPE(force_env_type), POINTER :: force_env
1390 
1391  INTEGER :: iforce_eval, iparticle, jparticle, &
1392  my_group, natom, nforce_eval
1393  INTEGER, DIMENSION(:), POINTER :: glob_natoms, map_index
1394  LOGICAL :: converged_embed
1395  REAL(kind=dp) :: energy
1396  REAL(kind=dp), DIMENSION(:), POINTER :: energies
1397  TYPE(cell_type), POINTER :: cell_embed
1398  TYPE(cp_logger_type), POINTER :: logger, my_logger
1399  TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1400  TYPE(cp_result_type), POINTER :: loc_results, results_embed
1401  TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1402  TYPE(cp_subsys_type), POINTER :: subsys_embed
1403  TYPE(dft_control_type), POINTER :: dft_control
1404  TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1405  TYPE(particle_list_type), POINTER :: particles_embed
1406  TYPE(pw_env_type), POINTER :: pw_env
1407  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1408  TYPE(pw_r3d_rs_type), POINTER :: embed_pot, spin_embed_pot
1409  TYPE(section_vals_type), POINTER :: embed_section, force_env_section, &
1410  mapping_section, root_section
1411 
1412  logger => cp_get_default_logger()
1413  cpassert(ASSOCIATED(force_env))
1414  ! Get infos about the embedding subsys
1415  CALL force_env_get(force_env=force_env, &
1416  subsys=subsys_embed, &
1417  force_env_section=force_env_section, &
1418  root_section=root_section, &
1419  cell=cell_embed)
1420  CALL cp_subsys_get(subsys=subsys_embed, &
1421  particles=particles_embed, &
1422  results=results_embed)
1423  NULLIFY (map_index, glob_natoms)
1424 
1425  nforce_eval = SIZE(force_env%sub_force_env)
1426  embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1427  mapping_section => section_vals_get_subs_vals(embed_section, "MAPPING")
1428  ! Global Info
1429  ALLOCATE (subsystems(nforce_eval))
1430  ALLOCATE (particles(nforce_eval))
1431  ! Local Info to sync
1432  ALLOCATE (energies(nforce_eval))
1433  ALLOCATE (glob_natoms(nforce_eval))
1434  ALLOCATE (results(nforce_eval))
1435  energies = 0.0_dp
1436  glob_natoms = 0
1437 
1438  DO iforce_eval = 1, nforce_eval
1439  NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1440  NULLIFY (results(iforce_eval)%results)
1441  CALL cp_result_create(results(iforce_eval)%results)
1442  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1443  ! From this point on the error is the sub_error
1444  my_group = force_env%embed_env%group_distribution(force_env%para_env%mepos)
1445  my_logger => force_env%embed_env%sub_logger(my_group + 1)%p
1446  ! Copy iterations info (they are updated only in the main embed_env)
1447  CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1448  CALL cp_add_default_logger(my_logger)
1449 
1450  ! Get all available subsys
1451  CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1452  subsys=subsystems(iforce_eval)%subsys)
1453 
1454  ! Check if we import density from previous force calculations
1455  ! Only for QUICKSTEP
1456  IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env%qs_env)) THEN
1457  NULLIFY (dft_control)
1458  CALL get_qs_env(force_env%sub_force_env(iforce_eval)%force_env%qs_env, dft_control=dft_control)
1459  IF (dft_control%qs_control%ref_embed_subsys) THEN
1460  IF (iforce_eval .EQ. 2) cpabort("Density importing force_eval can't be the first.")
1461  END IF
1462  END IF
1463 
1464  ! all force_env share the same cell
1465  CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_embed)
1466 
1467  ! Get available particles
1468  CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1469  particles=particles(iforce_eval)%list)
1470 
1471  ! Get Mapping index array
1472  natom = SIZE(particles(iforce_eval)%list%els)
1473 
1474  CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1475  map_index, .true.)
1476 
1477  ! Mapping particles from iforce_eval environment to the embed env
1478  DO iparticle = 1, natom
1479  jparticle = map_index(iparticle)
1480  particles(iforce_eval)%list%els(iparticle)%r = particles_embed%els(jparticle)%r
1481  END DO
1482 
1483  ! Calculate energy and forces for each sub_force_env
1484  CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1485  calc_force=.false., &
1486  skip_external_control=.true.)
1487 
1488  ! Call DFT embedding
1489  IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env%qs_env)) THEN
1490  NULLIFY (dft_control)
1491  CALL get_qs_env(force_env%sub_force_env(iforce_eval)%force_env%qs_env, dft_control=dft_control)
1492  IF (dft_control%qs_control%ref_embed_subsys) THEN
1493  ! Now we can optimize the embedding potential
1494  CALL dft_embedding(force_env, iforce_eval, energies, converged_embed)
1495  IF (.NOT. converged_embed) cpabort("Embedding potential optimization not converged.")
1496  END IF
1497  ! Deallocate embedding potential on the high-level subsystem
1498  IF (dft_control%qs_control%high_level_embed_subsys) THEN
1499  CALL get_qs_env(qs_env=force_env%sub_force_env(iforce_eval)%force_env%qs_env, &
1500  embed_pot=embed_pot, spin_embed_pot=spin_embed_pot, pw_env=pw_env)
1501  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1502  CALL auxbas_pw_pool%give_back_pw(embed_pot)
1503  IF (ASSOCIATED(embed_pot)) THEN
1504  CALL embed_pot%release()
1505  DEALLOCATE (embed_pot)
1506  END IF
1507  IF (ASSOCIATED(spin_embed_pot)) THEN
1508  CALL auxbas_pw_pool%give_back_pw(spin_embed_pot)
1509  CALL spin_embed_pot%release()
1510  DEALLOCATE (spin_embed_pot)
1511  END IF
1512  END IF
1513  END IF
1514 
1515  ! Only the rank 0 process collect info for each computation
1516  IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1517  CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1518  potential_energy=energy)
1519  CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1520  results=loc_results)
1521  energies(iforce_eval) = energy
1522  glob_natoms(iforce_eval) = natom
1523  CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1524  END IF
1525  ! Deallocate map_index array
1526  IF (ASSOCIATED(map_index)) THEN
1527  DEALLOCATE (map_index)
1528  END IF
1529  CALL cp_rm_default_logger()
1530  END DO
1531 
1532  ! Handling Parallel execution
1533  CALL force_env%para_env%sync()
1534  ! Let's transfer energy, natom
1535  CALL force_env%para_env%sum(energies)
1536  CALL force_env%para_env%sum(glob_natoms)
1537 
1538  force_env%embed_env%energies = energies
1539 
1540  !NB if the first system has fewer atoms than the second)
1541  DO iparticle = 1, SIZE(particles_embed%els)
1542  particles_embed%els(iparticle)%f(:) = 0.0_dp
1543  END DO
1544 
1545  ! ONIOM type of mixing in embedding: E = E_total + E_cluster_high - E_cluster
1546  force_env%embed_env%pot_energy = energies(3) + energies(4) - energies(2)
1547 
1548  !Simply deallocate and loose the pointer references..
1549  DO iforce_eval = 1, nforce_eval
1550  CALL cp_result_release(results(iforce_eval)%results)
1551  END DO
1552  DEALLOCATE (subsystems)
1553  DEALLOCATE (particles)
1554  DEALLOCATE (energies)
1555  DEALLOCATE (glob_natoms)
1556  DEALLOCATE (results)
1557 
1558  END SUBROUTINE embed_energy
1559 
1560 ! **************************************************************************************************
1561 !> \brief ...
1562 !> \param force_env ...
1563 !> \param ref_subsys_number ...
1564 !> \param energies ...
1565 !> \param converged_embed ...
1566 ! **************************************************************************************************
1567  SUBROUTINE dft_embedding(force_env, ref_subsys_number, energies, converged_embed)
1568  TYPE(force_env_type), POINTER :: force_env
1569  INTEGER :: ref_subsys_number
1570  REAL(kind=dp), DIMENSION(:), POINTER :: energies
1571  LOGICAL :: converged_embed
1572 
1573  INTEGER :: embed_method
1574  TYPE(section_vals_type), POINTER :: embed_section, force_env_section
1575 
1576  ! Find out which embedding scheme is used
1577  CALL force_env_get(force_env=force_env, &
1578  force_env_section=force_env_section)
1579  embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1580 
1581  CALL section_vals_val_get(embed_section, "EMBED_METHOD", i_val=embed_method)
1582  SELECT CASE (embed_method)
1583  CASE (dfet)
1584  ! Density functional embedding
1585  CALL dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1586  CASE (dmfet)
1587  ! Density matrix embedding theory
1588  CALL dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1589  END SELECT
1590 
1591  END SUBROUTINE dft_embedding
1592 ! **************************************************************************************************
1593 !> \brief ... Main driver for DFT embedding
1594 !> \param force_env ...
1595 !> \param ref_subsys_number ...
1596 !> \param energies ...
1597 !> \param converged_embed ...
1598 !> \author Vladimir Rybkin
1599 ! **************************************************************************************************
1600  SUBROUTINE dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1601  TYPE(force_env_type), POINTER :: force_env
1602  INTEGER :: ref_subsys_number
1603  REAL(kind=dp), DIMENSION(:), POINTER :: energies
1604  LOGICAL :: converged_embed
1605 
1606  CHARACTER(LEN=*), PARAMETER :: routinen = 'dfet_embedding'
1607 
1608  INTEGER :: cluster_subsys_num, handle, &
1609  i_force_eval, i_iter, i_spin, &
1610  nforce_eval, nspins, nspins_subsys, &
1611  output_unit
1612  REAL(kind=dp) :: cluster_energy
1613  REAL(kind=dp), DIMENSION(:), POINTER :: rhs
1614  TYPE(cp_logger_type), POINTER :: logger
1615  TYPE(dft_control_type), POINTER :: dft_control
1616  TYPE(opt_embed_pot_type) :: opt_embed
1617  TYPE(pw_env_type), POINTER :: pw_env
1618  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1619  TYPE(pw_r3d_rs_type) :: diff_rho_r, diff_rho_spin
1620  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_ref, rho_r_subsys
1621  TYPE(pw_r3d_rs_type), POINTER :: embed_pot, embed_pot_subsys, &
1622  spin_embed_pot, spin_embed_pot_subsys
1623  TYPE(qs_energy_type), POINTER :: energy
1624  TYPE(qs_rho_type), POINTER :: rho, subsys_rho
1625  TYPE(section_vals_type), POINTER :: dft_section, embed_section, &
1626  force_env_section, input, &
1627  mapping_section, opt_embed_section
1628 
1629  CALL timeset(routinen, handle)
1630 
1631  CALL cite_reference(huang2011)
1632  CALL cite_reference(heaton_burgess2007)
1633 
1634  CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
1635 
1636  ! Reveal input file
1637  NULLIFY (logger)
1638  logger => cp_get_default_logger()
1639  output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
1640  extension=".Log")
1641 
1642  NULLIFY (dft_section, input, opt_embed_section)
1643  NULLIFY (energy, dft_control)
1644 
1645  CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1646  pw_env=pw_env, dft_control=dft_control, rho=rho, energy=energy, &
1647  input=input)
1648  nspins = dft_control%nspins
1649 
1650  dft_section => section_vals_get_subs_vals(input, "DFT")
1651  opt_embed_section => section_vals_get_subs_vals(input, &
1652  "DFT%QS%OPT_EMBED")
1653  ! Rho_r is the reference
1654  CALL qs_rho_get(rho_struct=rho, rho_r=rho_r_ref)
1655 
1656  ! We need to understand how to treat spins states
1657  CALL understand_spin_states(force_env, ref_subsys_number, opt_embed%change_spin, opt_embed%open_shell_embed, &
1658  opt_embed%all_nspins)
1659 
1660  ! Prepare everything for the potential maximization
1661  CALL prepare_embed_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, opt_embed, &
1662  opt_embed_section)
1663 
1664  ! Initialize embedding potential
1665  CALL init_embed_pot(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, &
1666  opt_embed%add_const_pot, opt_embed%Fermi_Amaldi, opt_embed%const_pot, &
1667  opt_embed%open_shell_embed, spin_embed_pot, &
1668  opt_embed%pot_diff, opt_embed%Coulomb_guess, opt_embed%grid_opt)
1669 
1670  ! Read embedding potential vector from the file
1671  IF (opt_embed%read_embed_pot .OR. opt_embed%read_embed_pot_cube) CALL read_embed_pot( &
1672  force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, spin_embed_pot, &
1673  opt_embed_section, opt_embed)
1674 
1675  ! Prepare the pw object to store density differences
1676  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1677  CALL auxbas_pw_pool%create_pw(diff_rho_r)
1678  CALL pw_zero(diff_rho_r)
1679  IF (opt_embed%open_shell_embed) THEN
1680  CALL auxbas_pw_pool%create_pw(diff_rho_spin)
1681  CALL pw_zero(diff_rho_spin)
1682  END IF
1683 
1684  ! Check the preliminary density differences
1685  DO i_spin = 1, nspins
1686  CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
1687  END DO
1688  IF (opt_embed%open_shell_embed) THEN ! Spin part
1689  IF (nspins .EQ. 2) THEN ! Reference systems has an open shell, else the reference diff_rho_spin is zero
1690  CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
1691  CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
1692  END IF
1693  END IF
1694 
1695  DO i_force_eval = 1, ref_subsys_number - 1
1696  NULLIFY (subsys_rho, rho_r_subsys, dft_control)
1697  CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, energy=energy, &
1698  dft_control=dft_control)
1699  nspins_subsys = dft_control%nspins
1700  ! Add subsystem densities
1701  CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
1702  DO i_spin = 1, nspins_subsys
1703  CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
1704  END DO
1705  IF (opt_embed%open_shell_embed) THEN ! Spin part
1706  IF (nspins_subsys .EQ. 2) THEN ! The subsystem makes contribution if it is spin-polarized
1707  ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
1708  IF ((i_force_eval .EQ. 2) .AND. (opt_embed%change_spin)) THEN
1709  CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1710  CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1711  ELSE
1712  ! First subsystem (always) and second subsystem (without spin change)
1713  CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1714  CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1715  END IF
1716  END IF
1717  END IF
1718  END DO
1719 
1720  ! Print density difference
1721  CALL print_rho_diff(diff_rho_r, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1722  IF (opt_embed%open_shell_embed) THEN ! Spin part
1723  CALL print_rho_spin_diff(diff_rho_spin, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1724  END IF
1725 
1726  ! Construct electrostatic guess if needed
1727  IF (opt_embed%Coulomb_guess) THEN
1728  ! Reveal resp charges for total system
1729  nforce_eval = SIZE(force_env%sub_force_env)
1730  NULLIFY (rhs)
1731  CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, rhs=rhs)
1732  ! Get the mapping
1733  CALL force_env_get(force_env=force_env, &
1734  force_env_section=force_env_section)
1735  embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1736  mapping_section => section_vals_get_subs_vals(embed_section, "MAPPING")
1737 
1738  DO i_force_eval = 1, ref_subsys_number - 1
1739  IF (i_force_eval .EQ. 1) THEN
1740  CALL coulomb_guess(embed_pot, rhs, mapping_section, &
1741  force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1742  ELSE
1743  CALL coulomb_guess(opt_embed%pot_diff, rhs, mapping_section, &
1744  force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1745  END IF
1746  END DO
1747  CALL pw_axpy(opt_embed%pot_diff, embed_pot)
1748  IF (.NOT. opt_embed%grid_opt) CALL pw_copy(embed_pot, opt_embed%const_pot)
1749 
1750  END IF
1751 
1752  ! Difference guess
1753  IF (opt_embed%diff_guess) THEN
1754  CALL pw_copy(diff_rho_r, embed_pot)
1755  IF (.NOT. opt_embed%grid_opt) CALL pw_copy(embed_pot, opt_embed%const_pot)
1756  ! Open shell
1757  IF (opt_embed%open_shell_embed) CALL pw_copy(diff_rho_spin, spin_embed_pot)
1758  END IF
1759 
1760  ! Calculate subsystems with trial embedding potential
1761  DO i_iter = 1, opt_embed%n_iter
1762  opt_embed%i_iter = i_iter
1763 
1764  ! Set the density difference as the negative reference one
1765  CALL pw_zero(diff_rho_r)
1766  CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, dft_control=dft_control)
1767  nspins = dft_control%nspins
1768  DO i_spin = 1, nspins
1769  CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
1770  END DO
1771  IF (opt_embed%open_shell_embed) THEN ! Spin part
1772  CALL pw_zero(diff_rho_spin)
1773  IF (nspins .EQ. 2) THEN ! Reference systems has an open shell, else the reference diff_rho_spin is zero
1774  CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
1775  CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
1776  END IF
1777  END IF
1778 
1779  DO i_force_eval = 1, ref_subsys_number - 1
1780  NULLIFY (dft_control)
1781  CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
1782  nspins_subsys = dft_control%nspins
1783 
1784  IF ((i_force_eval .EQ. 2) .AND. (opt_embed%change_spin)) THEN
1785  ! Here we change the sign of the spin embedding potential due to spin change:
1786  ! only in spin_embed_subsys
1787  CALL make_subsys_embed_pot(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
1788  embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
1789  opt_embed%open_shell_embed, .true.)
1790  ELSE ! Regular case
1791  CALL make_subsys_embed_pot(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
1792  embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
1793  opt_embed%open_shell_embed, .false.)
1794  END IF
1795 
1796  ! Switch on external potential in the subsystems
1797  dft_control%apply_embed_pot = .true.
1798 
1799  ! Add the embedding potential
1800  CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, embed_pot=embed_pot_subsys)
1801  IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys .EQ. 2)) THEN ! Spin part
1802  CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
1803  spin_embed_pot=spin_embed_pot_subsys)
1804  END IF
1805 
1806  ! Get the previous subsystem densities
1807  CALL get_prev_density(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
1808 
1809  ! Calculate the new density
1810  CALL force_env_calc_energy_force(force_env=force_env%sub_force_env(i_force_eval)%force_env, &
1811  calc_force=.false., &
1812  skip_external_control=.true.)
1813 
1814  CALL get_max_subsys_diff(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
1815 
1816  ! Extract subsystem density and energy
1817  NULLIFY (rho_r_subsys, energy)
1818 
1819  CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, &
1820  energy=energy)
1821  opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + energy%total
1822 
1823  ! Find out which subsystem is the cluster
1824  IF (dft_control%qs_control%cluster_embed_subsys) THEN
1825  cluster_subsys_num = i_force_eval
1826  cluster_energy = energy%total
1827  END IF
1828 
1829  ! Add subsystem densities
1830  CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
1831  DO i_spin = 1, nspins_subsys
1832  CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
1833  END DO
1834  IF (opt_embed%open_shell_embed) THEN ! Spin part
1835  IF (nspins_subsys .EQ. 2) THEN ! The subsystem makes contribution if it is spin-polarized
1836  ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
1837  IF ((i_force_eval .EQ. 2) .AND. (opt_embed%change_spin)) THEN
1838  CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1839  CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1840  ELSE
1841  ! First subsystem (always) and second subsystem (without spin change)
1842  CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1843  CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1844  END IF
1845  END IF
1846  END IF
1847 
1848  ! Release embedding potential for subsystem
1849  CALL embed_pot_subsys%release()
1850  DEALLOCATE (embed_pot_subsys)
1851  IF (opt_embed%open_shell_embed) THEN
1852  CALL spin_embed_pot_subsys%release()
1853  DEALLOCATE (spin_embed_pot_subsys)
1854  END IF
1855 
1856  END DO ! i_force_eval
1857 
1858  ! Print embedding potential for restart
1859  CALL print_embed_restart(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1860  opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
1861  spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .false.)
1862  CALL print_pot_simple_grid(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1863  embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .false., &
1864  force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
1865 
1866  ! Integrate the potential over density differences and add to w functional; also add regularization contribution
1867  DO i_spin = 1, nspins ! Sum over nspins for the reference system, not subsystem!
1868  opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) - pw_integral_ab(embed_pot, rho_r_ref(i_spin))
1869  END DO
1870  ! Spin part
1871  IF (opt_embed%open_shell_embed) THEN
1872  ! If reference system is not spin-polarized then it does not make a contribution to W functional
1873  IF (nspins .EQ. 2) THEN
1874  opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) &
1875  - pw_integral_ab(spin_embed_pot, rho_r_ref(1)) &
1876  + pw_integral_ab(spin_embed_pot, rho_r_ref(2))
1877  END IF
1878  END IF
1879  ! Finally, add the regularization term
1880  opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + opt_embed%reg_term
1881 
1882  ! Print information and check convergence
1883  CALL print_emb_opt_info(output_unit, i_iter, opt_embed)
1884  CALL conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit)
1885  IF (opt_embed%converged) EXIT
1886 
1887  ! Update the trust radius and control the step
1888  IF ((i_iter .GT. 1) .AND. (.NOT. opt_embed%steep_desc)) CALL step_control(opt_embed)
1889 
1890  ! Print density difference
1891  CALL print_rho_diff(diff_rho_r, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1892  IF (opt_embed%open_shell_embed) THEN ! Spin part
1893  CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1894  END IF
1895 
1896  ! Calculate potential gradient if the step has been accepted. Otherwise, we reuse the previous one
1897 
1898  IF (opt_embed%accept_step .AND. (.NOT. opt_embed%grid_opt)) &
1899  CALL calculate_embed_pot_grad(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1900  diff_rho_r, diff_rho_spin, opt_embed)
1901  ! Take the embedding step
1902  CALL opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, &
1903  force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
1904 
1905  END DO ! i_iter
1906 
1907  ! Print final embedding potential for restart
1908  CALL print_embed_restart(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1909  opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
1910  spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .true.)
1911  CALL print_pot_simple_grid(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1912  embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .true., &
1913  force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
1914 
1915  ! Print final density difference
1916  !CALL print_rho_diff(diff_rho_r, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .TRUE.)
1917  IF (opt_embed%open_shell_embed) THEN ! Spin part
1918  CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .true.)
1919  END IF
1920 
1921  ! Give away plane waves pools
1922  CALL diff_rho_r%release()
1923  IF (opt_embed%open_shell_embed) THEN
1924  CALL diff_rho_spin%release()
1925  END IF
1926 
1927  CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
1928  "PRINT%PROGRAM_RUN_INFO")
1929 
1930  ! If converged send the embedding potential to the higher-level calculation.
1931  IF (opt_embed%converged) THEN
1932  CALL get_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, dft_control=dft_control, &
1933  pw_env=pw_env)
1934  nspins_subsys = dft_control%nspins
1935  dft_control%apply_embed_pot = .true.
1936  ! The embedded subsystem corresponds to subsystem #2, where spin change is possible
1937  CALL make_subsys_embed_pot(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, &
1938  embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
1939  opt_embed%open_shell_embed, opt_embed%change_spin)
1940 
1941  IF (opt_embed%Coulomb_guess) THEN
1942  CALL pw_axpy(opt_embed%pot_diff, embed_pot_subsys, -1.0_dp, allow_noncompatible_grids=.true.)
1943  END IF
1944 
1945  CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, embed_pot=embed_pot_subsys)
1946 
1947  IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys .EQ. 2)) THEN
1948  CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, &
1949  spin_embed_pot=spin_embed_pot_subsys)
1950  END IF
1951 
1952  ! Substitute the correct energy in energies: only on rank 0
1953  IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source()) THEN
1954  energies(cluster_subsys_num) = cluster_energy
1955  END IF
1956  END IF
1957 
1958  ! Deallocate and release opt_embed content
1959  CALL release_opt_embed(opt_embed)
1960 
1961  ! Deallocate embedding potential
1962  CALL embed_pot%release()
1963  DEALLOCATE (embed_pot)
1964  IF (opt_embed%open_shell_embed) THEN
1965  CALL spin_embed_pot%release()
1966  DEALLOCATE (spin_embed_pot)
1967  END IF
1968 
1969  converged_embed = opt_embed%converged
1970 
1971  CALL timestop(handle)
1972 
1973  END SUBROUTINE dfet_embedding
1974 
1975 ! **************************************************************************************************
1976 !> \brief Main driver for the DMFET embedding
1977 !> \param force_env ...
1978 !> \param ref_subsys_number ...
1979 !> \param energies ...
1980 !> \param converged_embed ...
1981 !> \author Vladimir Rybkin
1982 ! **************************************************************************************************
1983  SUBROUTINE dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1984  TYPE(force_env_type), POINTER :: force_env
1985  INTEGER :: ref_subsys_number
1986  REAL(kind=dp), DIMENSION(:), POINTER :: energies
1987  LOGICAL :: converged_embed
1988 
1989  CHARACTER(LEN=*), PARAMETER :: routinen = 'dmfet_embedding'
1990 
1991  INTEGER :: cluster_subsys_num, handle, &
1992  i_force_eval, i_iter, output_unit
1993  LOGICAL :: subsys_open_shell
1994  REAL(kind=dp) :: cluster_energy
1995  TYPE(cp_logger_type), POINTER :: logger
1996  TYPE(dft_control_type), POINTER :: dft_control
1997  TYPE(mp_para_env_type), POINTER :: para_env
1998  TYPE(opt_dmfet_pot_type) :: opt_dmfet
1999  TYPE(qs_energy_type), POINTER :: energy
2000  TYPE(section_vals_type), POINTER :: dft_section, input, opt_dmfet_section
2001 
2002  CALL timeset(routinen, handle)
2003 
2004  CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2005  para_env=para_env)
2006 
2007  ! Reveal input file
2008  NULLIFY (logger)
2009  logger => cp_get_default_logger()
2010  output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
2011  extension=".Log")
2012 
2013  NULLIFY (dft_section, input, opt_dmfet_section)
2014  NULLIFY (energy)
2015 
2016  CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2017  energy=energy, input=input)
2018 
2019  dft_section => section_vals_get_subs_vals(input, "DFT")
2020  opt_dmfet_section => section_vals_get_subs_vals(input, &
2021  "DFT%QS%OPT_DMFET")
2022 
2023  ! We need to understand how to treat spins states
2024  CALL understand_spin_states(force_env, ref_subsys_number, opt_dmfet%change_spin, opt_dmfet%open_shell_embed, &
2025  opt_dmfet%all_nspins)
2026 
2027  ! Prepare for the potential optimization
2028  CALL prepare_dmfet_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2029  opt_dmfet, opt_dmfet_section)
2030 
2031  ! Get the reference density matrix/matrices
2032  subsys_open_shell = subsys_spin(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2033  CALL build_full_dm(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2034  opt_dmfet%dm_total, subsys_open_shell, opt_dmfet%open_shell_embed, opt_dmfet%dm_total_beta)
2035 
2036  ! Check the preliminary DM difference
2037  CALL cp_fm_copy_general(opt_dmfet%dm_total, opt_dmfet%dm_diff, para_env)
2038  IF (opt_dmfet%open_shell_embed) CALL cp_fm_copy_general(opt_dmfet%dm_total_beta, &
2039  opt_dmfet%dm_diff_beta, para_env)
2040 
2041  DO i_force_eval = 1, ref_subsys_number - 1
2042 
2043  ! Get the subsystem density matrix/matrices
2044  subsys_open_shell = subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2045 
2046  CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2047  opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2048  opt_dmfet%dm_subsys_beta)
2049 
2050  CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2051 
2052  IF (opt_dmfet%open_shell_embed) CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, &
2053  1.0_dp, opt_dmfet%dm_subsys_beta)
2054 
2055  END DO
2056 
2057  ! Main loop of iterative matrix potential optimization
2058  DO i_iter = 1, opt_dmfet%n_iter
2059 
2060  opt_dmfet%i_iter = i_iter
2061 
2062  ! Set the dm difference as the reference one
2063  CALL cp_fm_copy_general(opt_dmfet%dm_total, opt_dmfet%dm_diff, para_env)
2064 
2065  IF (opt_dmfet%open_shell_embed) CALL cp_fm_copy_general(opt_dmfet%dm_total_beta, &
2066  opt_dmfet%dm_diff_beta, para_env)
2067 
2068  ! Loop over force evaluations
2069  DO i_force_eval = 1, ref_subsys_number - 1
2070 
2071  ! Switch on external potential in the subsystems
2072  NULLIFY (dft_control)
2073  CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
2074  dft_control%apply_dmfet_pot = .true.
2075 
2076  ! Calculate the new density
2077  CALL force_env_calc_energy_force(force_env=force_env%sub_force_env(i_force_eval)%force_env, &
2078  calc_force=.false., &
2079  skip_external_control=.true.)
2080 
2081  ! Extract subsystem density matrix and energy
2082  NULLIFY (energy)
2083 
2084  CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, energy=energy)
2085  opt_dmfet%w_func(i_iter) = opt_dmfet%w_func(i_iter) + energy%total
2086 
2087  ! Find out which subsystem is the cluster
2088  IF (dft_control%qs_control%cluster_embed_subsys) THEN
2089  cluster_subsys_num = i_force_eval
2090  cluster_energy = energy%total
2091  END IF
2092 
2093  ! Add subsystem density matrices
2094  subsys_open_shell = subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2095 
2096  CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2097  opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2098  opt_dmfet%dm_subsys_beta)
2099 
2100  IF (opt_dmfet%open_shell_embed) THEN ! Open-shell embedding
2101  ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
2102  IF ((i_force_eval .EQ. 2) .AND. (opt_dmfet%change_spin)) THEN
2103  CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, 1.0_dp, opt_dmfet%dm_subsys)
2104  CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys_beta)
2105  ELSE
2106  CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2107  CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, 1.0_dp, opt_dmfet%dm_subsys_beta)
2108  END IF
2109  ELSE ! Closed-shell embedding
2110  CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2111  END IF
2112 
2113  END DO ! i_force_eval
2114 
2115  CALL check_dmfet(opt_dmfet, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2116 
2117  END DO ! i_iter
2118 
2119  ! Substitute the correct energy in energies: only on rank 0
2120  IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source()) THEN
2121  energies(cluster_subsys_num) = cluster_energy
2122  END IF
2123 
2124  CALL release_dmfet_opt(opt_dmfet)
2125 
2126  converged_embed = .false.
2127 
2128  END SUBROUTINE dmfet_embedding
2129 
2130 END 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:117
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:66
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::huang2011
integer, save, public huang2011
Definition: bibliography.F:43
bibliography::heaton_burgess2007
integer, save, public heaton_burgess2007
Definition: bibliography.F:43
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:117
cell_methods::cell_create
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition: cell_methods.F:85
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:516
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:486
cell_types::cell_release
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition: cell_types.F:559
cell_types::cell_clone
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition: cell_types.F:107
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_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:101
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:1538
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:167
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:803
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:1013
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)
sets various propreties of the subsys
Definition: cp_subsys_types.F:205
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)
returns information about various attributes of the given subsys
Definition: cp_subsys_types.F:345
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_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_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:193
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:551
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)
creates and initializes a force environment
Definition: force_env_methods.F:759
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:511
force_env_types::use_qmmm
integer, parameter, public use_qmmm
Definition: force_env_types.F:79
force_env_types::use_prog_name
character(len=10), dimension(501:509), parameter, public use_prog_name
Definition: force_env_types.F:89
force_env_types::use_mixed_force
integer, parameter, public use_mixed_force
Definition: force_env_types.F:79
force_env_types::use_eip_force
integer, parameter, public use_eip_force
Definition: force_env_types.F:79
force_env_types::use_embed
integer, parameter, public use_embed
Definition: force_env_types.F:79
force_env_types::use_qmmmx
integer, parameter, public use_qmmmx
Definition: force_env_types.F:79
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)
returns various attributes about the force environment
Definition: force_env_types.F:329
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:640
force_env_types::use_qs_force
integer, parameter, public use_qs_force
Definition: force_env_types.F:79
force_env_types::use_pwdft_force
integer, parameter, public use_pwdft_force
Definition: force_env_types.F:79
force_env_types::use_nnp_force
integer, parameter, public use_nnp_force
Definition: force_env_types.F:79
force_env_types::use_fist_force
integer, parameter, public use_fist_force
Definition: force_env_types.F:79
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, total_force, grand_total_force, zero_force_core_shell_atom)
Write forces.
Definition: force_env_utils.F:415
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:491
force_env_utils::rescale_forces
subroutine, public rescale_forces(force_env)
Rescale forces if requested.
Definition: force_env_utils.F:369
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:1236
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)
computest the forces and the energy due to the flexible potential & bias, and writes the weights file
Definition: fp_methods.F:54
fparser
This public domain function parser module is intended for applications where a set of mathematical ex...
Definition: fparser.F:17
fparser::evalf
real(rn) function, public evalf(i, Val)
...
Definition: fparser.F:180
fparser::evalerrtype
integer, public evalerrtype
Definition: fparser.F:32
fparser::evalfd
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
Definition: fparser.F:976
fparser::finalizef
subroutine, public finalizef()
...
Definition: fparser.F:101
fparser::initf
subroutine, public initf(n)
...
Definition: fparser.F:130
fparser::parsef
subroutine, public parsef(i, FuncStr, Var)
Parse ith function string FuncStr and compile it into bytecode.
Definition: fparser.F:148
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:119
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:710
input_constants::mixed_cdft_serial
integer, parameter, public mixed_cdft_serial
Definition: input_constants.F:772
input_constants::mix_cdft
integer, parameter, public mix_cdft
Definition: input_constants.F:781
input_constants::mix_linear_combination
integer, parameter, public mix_linear_combination
Definition: input_constants.F:781
input_constants::dmfet
integer, parameter, public dmfet
Definition: input_constants.F:1246
input_constants::use_lenosky_eip
integer, parameter, public use_lenosky_eip
Definition: input_constants.F:710
input_constants::mix_coupled
integer, parameter, public mix_coupled
Definition: input_constants.F:781
input_constants::debug_run
integer, parameter, public debug_run
Definition: input_constants.F:121
input_constants::mix_restrained
integer, parameter, public mix_restrained
Definition: input_constants.F:781
input_constants::mix_generic
integer, parameter, public mix_generic
Definition: input_constants.F:781
input_constants::dfet
integer, parameter, public dfet
Definition: input_constants.F:1246
input_constants::mix_minimum
integer, parameter, public mix_minimum
Definition: input_constants.F:781
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:631
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:724
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:1040
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
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:332
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:150
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:510
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:3119
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:1024
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:3064
optimize_embedding_potential::get_max_subsys_diff
subroutine, public get_max_subsys_diff(opt_embed, force_env, subsys_num)
...
Definition: optimize_embedding_potential.F:3446
optimize_embedding_potential::print_emb_opt_info
subroutine, public print_emb_opt_info(output_unit, step_num, opt_embed)
...
Definition: optimize_embedding_potential.F:3383
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:128
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:663
optimize_embedding_potential::get_prev_density
subroutine, public get_prev_density(opt_embed, force_env, subsys_num)
...
Definition: optimize_embedding_potential.F:3416
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:933
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:3013
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:3479
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:2170
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:1467
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:1072
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:3217
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:291
optimize_embedding_potential::release_opt_embed
subroutine, public release_opt_embed(opt_embed)
Deallocate stuff for optimizing embedding potential.
Definition: optimize_embedding_potential.F:856
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:206
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
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:373
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_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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, 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, 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, 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, rhs)
Set the QUICKSTEP environment.
Definition: qs_environment_types.F:1035
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:107
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
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:3190
virial_methods
Definition: virial_methods.F:12
virial_methods::write_stress_tensor_components
subroutine, public write_stress_tensor_components(virial, iw, cell)
...
Definition: virial_methods.F:163
virial_methods::write_stress_tensor
subroutine, public write_stress_tensor(pv_virial, iw, cell, numerical)
Print stress tensor to output file.
Definition: virial_methods.F:254
virial_types
Definition: virial_types.F:13
virial_types::zero_virial
subroutine, public zero_virial(virial, reset)
...
Definition: virial_types.F:123
virial_types::symmetrize_virial
subroutine, public symmetrize_virial(virial)
Symmetrize the virial components.
Definition: virial_types.F:67