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