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