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