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