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