(git:374b731)
Loading...
Searching...
No Matches
force_env_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
17MODULE 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_mixed_force, use_nnp_force, use_prog_name, &
70 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 kahan_sum, ONLY: accurate_sum
93 USE kinds, ONLY: default_path_length,&
94 default_string_length,&
95 dp
96 USE machine, ONLY: m_memory
97 USE mathlib, ONLY: abnormal_value
98 USE message_passing, ONLY: mp_para_env_type
99 USE metadynamics_types, ONLY: meta_env_type
100 USE mixed_cdft_methods, ONLY: mixed_cdft_build_weight,&
101 mixed_cdft_calculate_coupling,&
102 mixed_cdft_init
103 USE mixed_energy_types, ONLY: mixed_energy_type,&
104 mixed_force_type
105 USE mixed_environment_types, ONLY: get_mixed_env,&
106 mixed_environment_type
107 USE mixed_environment_utils, ONLY: get_subsys_map_index,&
108 mixed_map_forces
109 USE molecule_kind_list_types, ONLY: molecule_kind_list_type
110 USE molecule_kind_types, ONLY: get_molecule_kind,&
111 molecule_kind_type
112 USE nnp_environment_types, ONLY: nnp_type
113 USE nnp_force, ONLY: nnp_calc_energy_force
114 USE optimize_dmfet_potential, ONLY: build_full_dm,&
115 check_dmfet,&
116 prepare_dmfet_opt,&
117 release_dmfet_opt,&
118 subsys_spin
119 USE optimize_embedding_potential, ONLY: &
120 coulomb_guess, calculate_embed_pot_grad, conv_check_embed, get_max_subsys_diff, &
121 get_prev_density, init_embed_pot, make_subsys_embed_pot, opt_embed_step, &
122 prepare_embed_opt, print_emb_opt_info, print_embed_restart, print_pot_simple_grid, &
123 print_rho_diff, print_rho_spin_diff, read_embed_pot, release_opt_embed, step_control, &
124 understand_spin_states
125 USE particle_list_types, ONLY: particle_list_p_type,&
126 particle_list_type
127 USE physcon, ONLY: debye
128 USE pw_env_types, ONLY: pw_env_get,&
129 pw_env_type
130 USE pw_methods, ONLY: pw_axpy,&
131 pw_copy,&
132 pw_integral_ab,&
133 pw_zero
134 USE pw_pool_types, ONLY: pw_pool_type
135 USE pw_types, ONLY: pw_r3d_rs_type
136 USE pwdft_environment, ONLY: pwdft_calc_energy_force
137 USE pwdft_environment_types, ONLY: pwdft_environment_type
138 USE qmmm_force, ONLY: qmmm_calc_energy_force
139 USE qmmm_types, ONLY: qmmm_env_type
140 USE qmmm_util, ONLY: apply_qmmm_translate
141 USE qmmmx_force, ONLY: qmmmx_calc_energy_force
142 USE qmmmx_types, ONLY: qmmmx_env_type
143 USE qs_energy_types, ONLY: qs_energy_type
144 USE qs_environment_types, ONLY: get_qs_env,&
145 qs_environment_type,&
146 set_qs_env
147 USE qs_force, ONLY: qs_calc_energy_force
148 USE qs_rho_types, ONLY: qs_rho_get,&
149 qs_rho_type
150 USE restraint, ONLY: restraint_control
151 USE scine_utils, ONLY: write_scine
152 USE string_utilities, ONLY: compress
153 USE virial_methods, ONLY: write_stress_tensor,&
154 write_stress_tensor_components
155 USE virial_types, ONLY: symmetrize_virial,&
156 virial_p_type,&
157 virial_type,&
158 zero_virial
159#include "./base/base_uses.f90"
160
161 IMPLICIT NONE
162
163 PRIVATE
164
165 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_env_methods'
166
167 PUBLIC :: force_env_create, &
168 force_env_calc_energy_force, &
169 force_env_calc_num_pressure
170
171 INTEGER, SAVE, PRIVATE :: last_force_env_id = 0
172
173CONTAINS
174
175! **************************************************************************************************
176!> \brief Interface routine for force and energy calculations
177!> \param force_env the force_env of which you want the energy and forces
178!> \param calc_force if false the forces *might* be left unchanged
179!> or be invalid, no guarantees can be given. Defaults to true
180!> \param consistent_energies Performs an additional qs_ks_update_qs_env, so
181!> that the energies are appropriate to the forces, they are in the
182!> non-selfconsistent case not consistent to each other! [08.2005, TdK]
183!> \param skip_external_control ...
184!> \param eval_energy_forces ...
185!> \param require_consistent_energy_force ...
186!> \param linres ...
187!> \param calc_stress_tensor ...
188!> \author CJM & fawzi
189! **************************************************************************************************
190 RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
191 consistent_energies, skip_external_control, eval_energy_forces, &
192 require_consistent_energy_force, linres, calc_stress_tensor)
193
194 TYPE(force_env_type), POINTER :: force_env
195 LOGICAL, INTENT(IN), OPTIONAL :: calc_force, consistent_energies, skip_external_control, &
196 eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor
197
198 REAL(kind=dp), PARAMETER :: ateps = 1.0e-6_dp
199
200 INTEGER :: i, ikind, j, nat, ndigits, nfixed_atoms, &
201 nfixed_atoms_total, nkind, &
202 output_unit, print_forces, print_grrm, &
203 print_scine
204 LOGICAL :: calculate_forces, calculate_stress_tensor, energy_consistency, eval_ef, &
205 linres_run, my_skip, print_components
206 REAL(kind=dp) :: checksum, e_entropy, e_gap, e_pot, &
207 sum_energy, sum_pv_virial, &
208 sum_stress_tensor
209 REAL(kind=dp), DIMENSION(3) :: grand_total_force, total_force
210 REAL(kind=dp), DIMENSION(3, 3) :: atomic_stress_tensor, diff_stress_tensor
211 TYPE(atprop_type), POINTER :: atprop_env
212 TYPE(cell_type), POINTER :: cell
213 TYPE(cp_logger_type), POINTER :: logger
214 TYPE(cp_subsys_type), POINTER :: subsys
215 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
216 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
217 TYPE(molecule_kind_type), POINTER :: molecule_kind
218 TYPE(particle_list_type), POINTER :: core_particles, particles, &
219 shell_particles
220 TYPE(virial_type), POINTER :: virial
221
222 NULLIFY (logger, virial, subsys, atprop_env, cell)
223 logger => cp_get_default_logger()
224 eval_ef = .true.
225 my_skip = .false.
226 calculate_forces = .true.
227 energy_consistency = .false.
228 linres_run = .false.
229 e_gap = -1.0_dp
230 e_entropy = -1.0_dp
231
232 IF (PRESENT(eval_energy_forces)) eval_ef = eval_energy_forces
233 IF (PRESENT(skip_external_control)) my_skip = skip_external_control
234 IF (PRESENT(calc_force)) calculate_forces = calc_force
235 IF (PRESENT(calc_stress_tensor)) THEN
236 calculate_stress_tensor = calc_stress_tensor
237 ELSE
238 calculate_stress_tensor = calculate_forces
239 END IF
240 IF (PRESENT(consistent_energies)) energy_consistency = consistent_energies
241 IF (PRESENT(linres)) linres_run = linres
242
243 cpassert(ASSOCIATED(force_env))
244 cpassert(force_env%ref_count > 0)
245 CALL force_env_get(force_env, subsys=subsys)
246 CALL force_env_set(force_env, additional_potential=0.0_dp)
247 CALL cp_subsys_get(subsys, virial=virial, atprop=atprop_env, cell=cell)
248 IF (virial%pv_availability) CALL zero_virial(virial, reset=.false.)
249
250 nat = force_env_get_natom(force_env)
251 CALL atprop_init(atprop_env, nat)
252 IF (eval_ef) THEN
253 SELECT CASE (force_env%in_use)
254 CASE (use_fist_force)
255 CALL fist_calc_energy_force(force_env%fist_env)
256 CASE (use_qs_force)
257 CALL qs_calc_energy_force(force_env%qs_env, calculate_forces, energy_consistency, linres_run)
258 CASE (use_pwdft_force)
259 IF (virial%pv_availability .AND. calculate_stress_tensor) THEN
260 CALL pwdft_calc_energy_force(force_env%pwdft_env, calculate_forces,.NOT. virial%pv_numer)
261 ELSE
262 CALL pwdft_calc_energy_force(force_env%pwdft_env, calculate_forces, .false.)
263 END IF
264 e_gap = force_env%pwdft_env%energy%band_gap
265 e_entropy = force_env%pwdft_env%energy%entropy
266 CASE (use_eip_force)
267 IF (force_env%eip_env%eip_model == use_lenosky_eip) THEN
268 CALL eip_lenosky(force_env%eip_env)
269 ELSE IF (force_env%eip_env%eip_model == use_bazant_eip) THEN
270 CALL eip_bazant(force_env%eip_env)
271 END IF
272 CASE (use_qmmm)
273 CALL qmmm_calc_energy_force(force_env%qmmm_env, &
274 calculate_forces, energy_consistency, linres=linres_run)
275 CASE (use_qmmmx)
276 CALL qmmmx_calc_energy_force(force_env%qmmmx_env, &
277 calculate_forces, energy_consistency, linres=linres_run, &
278 require_consistent_energy_force=require_consistent_energy_force)
279 CASE (use_mixed_force)
280 CALL mixed_energy_forces(force_env, calculate_forces)
281 CASE (use_nnp_force)
282 CALL nnp_calc_energy_force(force_env%nnp_env, &
283 calculate_forces)
284 CASE (use_embed)
285 CALL embed_energy(force_env)
286 CASE default
287 cpabort("")
288 END SELECT
289 END IF
290 ! In case it is requested, we evaluate the stress tensor numerically
291 IF (virial%pv_availability) THEN
292 IF (virial%pv_numer .AND. calculate_stress_tensor) THEN
293 ! Compute the numerical stress tensor
294 CALL force_env_calc_num_pressure(force_env)
295 ELSE
296 IF (calculate_forces) THEN
297 ! Symmetrize analytical stress tensor
298 CALL symmetrize_virial(virial)
299 ELSE
300 IF (calculate_stress_tensor) THEN
301 CALL cp_warn(__location__, "The calculation of the stress tensor "// &
302 "requires the calculation of the forces")
303 END IF
304 END IF
305 END IF
306 END IF
307
308 !sample peak memory
309 CALL m_memory()
310
311 ! Some additional tasks..
312 IF (.NOT. my_skip) THEN
313 ! Flexible Partitioning
314 IF (ASSOCIATED(force_env%fp_env)) THEN
315 IF (force_env%fp_env%use_fp) THEN
316 CALL fp_eval(force_env%fp_env, subsys, cell)
317 END IF
318 END IF
319 ! Constraints ONLY of Fixed Atom type
320 CALL fix_atom_control(force_env)
321 ! All Restraints
322 CALL restraint_control(force_env)
323 ! Virtual Sites
324 CALL vsite_force_control(force_env)
325 ! External Potential
326 CALL add_external_potential(force_env)
327 ! Rescale forces if requested
328 CALL rescale_forces(force_env)
329 END IF
330
331 CALL force_env_get(force_env, potential_energy=e_pot)
332
333 ! Print always Energy in the same format for all methods
334 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
335 extension=".Log")
336 IF (output_unit > 0) THEN
337 WRITE (output_unit, '(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) energy [a.u.]: ",T55,F26.15)') &
338 adjustr(trim(use_prog_name(force_env%in_use))), e_pot
339 IF (e_gap .GT. -.1_dp) THEN
340 WRITE (output_unit, '(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) gap [a.u.]: ",T55,F26.15)') &
341 adjustr(trim(use_prog_name(force_env%in_use))), e_gap
342 END IF
343 IF (e_entropy .GT. -0.1_dp) THEN
344 WRITE (output_unit, '(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) free energy [a.u.]: ",T55,F26.15)') &
345 adjustr(trim(use_prog_name(force_env%in_use))), e_pot - e_entropy
346 END IF
347 END IF
348 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
349 "PRINT%PROGRAM_RUN_INFO")
350
351 ! terminate the run if the value of the potential is abnormal
352 IF (abnormal_value(e_pot)) &
353 cpabort("Potential energy is an abnormal value (NaN/Inf).")
354
355 ! Print forces, if requested
356 print_forces = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%FORCES", &
357 extension=".xyz")
358 IF ((print_forces > 0) .AND. calculate_forces) THEN
359 CALL force_env_get(force_env, subsys=subsys)
360 CALL cp_subsys_get(subsys, &
361 core_particles=core_particles, &
362 particles=particles, &
363 shell_particles=shell_particles)
364 ! Variable precision output of the forces
365 CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%NDIGITS", &
366 i_val=ndigits)
367 IF (ASSOCIATED(core_particles) .OR. ASSOCIATED(shell_particles)) THEN
368 CALL write_forces(particles, print_forces, "ATOMIC", ndigits, total_force, &
369 zero_force_core_shell_atom=.true.)
370 grand_total_force(1:3) = total_force(1:3)
371 IF (ASSOCIATED(core_particles)) THEN
372 CALL write_forces(core_particles, print_forces, "CORE", ndigits, total_force, &
373 zero_force_core_shell_atom=.false.)
374 grand_total_force(:) = grand_total_force(:) + total_force(:)
375 END IF
376 IF (ASSOCIATED(shell_particles)) THEN
377 CALL write_forces(shell_particles, print_forces, "SHELL", ndigits, total_force, &
378 zero_force_core_shell_atom=.false., &
379 grand_total_force=grand_total_force)
380 END IF
381 ELSE
382 CALL write_forces(particles, print_forces, "ATOMIC", ndigits, total_force)
383 END IF
384 END IF
385 CALL cp_print_key_finished_output(print_forces, logger, force_env%force_env_section, "PRINT%FORCES")
386
387 ! Write stress tensor
388 IF (virial%pv_availability) THEN
389 ! If the virial is defined but we are not computing forces let's zero the
390 ! virial for consistency
391 IF (calculate_forces .AND. calculate_stress_tensor) THEN
392 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
393 extension=".stress_tensor")
394 IF (output_unit > 0) THEN
395 CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%COMPONENTS", &
396 l_val=print_components)
397 IF (print_components) THEN
398 IF ((.NOT. virial%pv_numer) .AND. (force_env%in_use == use_qs_force)) THEN
399 CALL write_stress_tensor_components(virial, output_unit, cell)
400 END IF
401 END IF
402 CALL write_stress_tensor(virial%pv_virial, output_unit, cell, virial%pv_numer)
403 END IF
404 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
405 "PRINT%STRESS_TENSOR")
406 ELSE
407 CALL zero_virial(virial, reset=.false.)
408 END IF
409 ELSE
410 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
411 extension=".stress_tensor")
412 IF (output_unit > 0) THEN
413 CALL cp_warn(__location__, "To print the stress tensor switch on the "// &
414 "virial evaluation with the keyword: STRESS_TENSOR")
415 END IF
416 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
417 "PRINT%STRESS_TENSOR")
418 END IF
419
420 ! Atomic energy
421 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
422 extension=".Log")
423 IF (atprop_env%energy) THEN
424 CALL force_env%para_env%sum(atprop_env%atener)
425 CALL force_env_get(force_env, potential_energy=e_pot)
426 IF (output_unit > 0) THEN
427 IF (logger%iter_info%print_level >= low_print_level) THEN
428 CALL cp_subsys_get(subsys=subsys, particles=particles)
429 CALL write_atener(output_unit, particles, atprop_env%atener, "Mulliken Atomic Energies")
430 END IF
431 sum_energy = accurate_sum(atprop_env%atener(:))
432 checksum = abs(e_pot - sum_energy)
433 WRITE (unit=output_unit, fmt="(/,(T2,A,T56,F25.13))") &
434 "Potential energy (Atomic):", sum_energy, &
435 "Potential energy (Total) :", e_pot, &
436 "Difference :", checksum
437 cpassert((checksum < ateps*abs(e_pot)))
438 END IF
439 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
440 "PRINT%PROGRAM_RUN_INFO")
441 END IF
442
443 ! Print GRMM interface file
444 print_grrm = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%GRRM", &
445 file_position="REWIND", extension=".rrm")
446 IF (print_grrm > 0) THEN
447 CALL force_env_get(force_env, subsys=subsys)
448 CALL cp_subsys_get(subsys=subsys, particles=particles, &
449 molecule_kinds=molecule_kinds)
450 ! Count the number of fixed atoms
451 nfixed_atoms_total = 0
452 nkind = molecule_kinds%n_els
453 molecule_kind_set => molecule_kinds%els
454 DO ikind = 1, nkind
455 molecule_kind => molecule_kind_set(ikind)
456 CALL get_molecule_kind(molecule_kind, nfixd=nfixed_atoms)
457 nfixed_atoms_total = nfixed_atoms_total + nfixed_atoms
458 END DO
459 !
460 CALL write_grrm(print_grrm, force_env, particles%els, e_pot, fixed_atoms=nfixed_atoms_total)
461 END IF
462 CALL cp_print_key_finished_output(print_grrm, logger, force_env%force_env_section, "PRINT%GRRM")
463
464 print_scine = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%SCINE", &
465 file_position="REWIND", extension=".scine")
466 IF (print_scine > 0) THEN
467 CALL force_env_get(force_env, subsys=subsys)
468 CALL cp_subsys_get(subsys=subsys, particles=particles)
469 !
470 CALL write_scine(print_scine, force_env, particles%els, e_pot)
471 END IF
472 CALL cp_print_key_finished_output(print_scine, logger, force_env%force_env_section, "PRINT%SCINE")
473
474 ! Atomic stress
475 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
476 extension=".Log")
477
478 IF (atprop_env%stress) THEN
479 CALL force_env%para_env%sum(atprop_env%atstress)
480 ! symmetrize (same as pv_virial)
481 DO i = 1, SIZE(atprop_env%atstress, 3)
482 atprop_env%atstress(:, :, i) = 0.5_dp*(atprop_env%atstress(:, :, i) &
483 + transpose(atprop_env%atstress(:, :, i)))
484 END DO
485 IF (output_unit > 0) THEN
486 IF (logger%iter_info%print_level > low_print_level) THEN
487 DO i = 1, SIZE(atprop_env%atstress, 3)
488 WRITE (unit=output_unit, fmt="(/,T2,I0,T16,A1,2(19X,A1))") i, "X", "Y", "Z"
489 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "X", (atprop_env%atstress(1, j, i), j=1, 3)
490 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Y", (atprop_env%atstress(2, j, i), j=1, 3)
491 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Z", (atprop_env%atstress(3, j, i), j=1, 3)
492 WRITE (unit=output_unit, fmt="(T2,A,F20.13)") "1/3 Trace(Atomic stress tensor):", &
493 (atprop_env%atstress(1, 1, i) + atprop_env%atstress(2, 2, i) + atprop_env%atstress(3, 3, i))/3.0_dp
494 END DO
495 END IF
496 atomic_stress_tensor(:, :) = 0.0_dp
497 DO i = 1, 3
498 atomic_stress_tensor(i, i) = accurate_sum(atprop_env%atstress(i, i, :))
499 DO j = i + 1, 3
500 atomic_stress_tensor(i, j) = accurate_sum(atprop_env%atstress(i, j, :))
501 atomic_stress_tensor(j, i) = atomic_stress_tensor(i, j)
502 END DO
503 END DO
504 WRITE (unit=output_unit, fmt="(/,T2,A,T15,A1,2(19X,A1))") "Atomic", "X", "Y", "Z"
505 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "X", (atomic_stress_tensor(1, i), i=1, 3)
506 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Y", (atomic_stress_tensor(2, i), i=1, 3)
507 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Z", (atomic_stress_tensor(3, i), i=1, 3)
508 WRITE (unit=output_unit, fmt="(T2,A,10X,F20.13)") "1/3 Trace(Atomic stress tensor):", &
509 (atomic_stress_tensor(1, 1) + atomic_stress_tensor(2, 2) + atomic_stress_tensor(3, 3))/3.0_dp
510 sum_stress_tensor = accurate_sum(atomic_stress_tensor(:, :))
511 IF (virial%pv_availability .AND. calculate_forces) THEN
512 WRITE (unit=output_unit, fmt="(/,T2,A,T16,A1,2(19X,A1))") "Total", "X", "Y", "Z"
513 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "X", (virial%pv_virial(1, i), i=1, 3)
514 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Y", (virial%pv_virial(2, i), i=1, 3)
515 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Z", (virial%pv_virial(3, i), i=1, 3)
516 WRITE (unit=output_unit, fmt="(T2,A,10X,F20.13)") "1/3 Trace(Total stress tensor): ", &
517 (virial%pv_virial(1, 1) + virial%pv_virial(2, 2) + virial%pv_virial(3, 3))/3.0_dp
518 sum_pv_virial = sum(virial%pv_virial(:, :))
519 diff_stress_tensor(:, :) = abs(virial%pv_virial(:, :) - atomic_stress_tensor(:, :))
520 WRITE (unit=output_unit, fmt="(/,T2,A,T16,A1,2(19X,A1))") "Diff", "X", "Y", "Z"
521 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "X", (diff_stress_tensor(1, i), i=1, 3)
522 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Y", (diff_stress_tensor(2, i), i=1, 3)
523 WRITE (unit=output_unit, fmt="(A3,3F20.13)") "Z", (diff_stress_tensor(3, i), i=1, 3)
524 WRITE (unit=output_unit, fmt="(T2,A,10X,F20.13)") "1/3 Trace(Diff) : ", &
525 (diff_stress_tensor(1, 1) + diff_stress_tensor(2, 2) + diff_stress_tensor(3, 3))/3.0_dp
526 checksum = accurate_sum(diff_stress_tensor(:, :))
527 WRITE (unit=output_unit, fmt="(/,(T2,A,11X,F25.13))") &
528 "Checksum stress (Atomic) :", sum_stress_tensor, &
529 "Checksum stress (Total) :", sum_pv_virial, &
530 "Difference :", checksum
531 cpassert(checksum < ateps)
532 END IF
533 END IF
534 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
535 "PRINT%PROGRAM_RUN_INFO")
536 END IF
537
538 END SUBROUTINE force_env_calc_energy_force
539
540! **************************************************************************************************
541!> \brief Evaluates the stress tensor and pressure numerically
542!> \param force_env ...
543!> \param dx ...
544!> \par History
545!> 10.2005 created [JCS]
546!> 05.2009 Teodoro Laino [tlaino] - rewriting for general force_env
547!>
548!> \author JCS
549! **************************************************************************************************
550 SUBROUTINE force_env_calc_num_pressure(force_env, dx)
551
552 TYPE(force_env_type), POINTER :: force_env
553 REAL(kind=dp), INTENT(IN), OPTIONAL :: dx
554
555 REAL(kind=dp), PARAMETER :: default_dx = 0.001_dp
556
557 INTEGER :: i, ip, iq, j, k, natom, ncore, nshell, &
558 output_unit, symmetry_id
559 REAL(kind=dp) :: dx_w
560 REAL(kind=dp), DIMENSION(2) :: numer_energy
561 REAL(kind=dp), DIMENSION(3) :: s
562 REAL(kind=dp), DIMENSION(3, 3) :: numer_stress
563 REAL(kind=dp), DIMENSION(:, :), POINTER :: ref_pos_atom, ref_pos_core, ref_pos_shell
564 TYPE(cell_type), POINTER :: cell, cell_local
565 TYPE(cp_logger_type), POINTER :: logger
566 TYPE(cp_subsys_type), POINTER :: subsys
567 TYPE(global_environment_type), POINTER :: globenv
568 TYPE(particle_list_type), POINTER :: core_particles, particles, &
569 shell_particles
570 TYPE(virial_type), POINTER :: virial
571
572 NULLIFY (cell_local)
573 NULLIFY (core_particles)
574 NULLIFY (particles)
575 NULLIFY (shell_particles)
576 NULLIFY (ref_pos_atom)
577 NULLIFY (ref_pos_core)
578 NULLIFY (ref_pos_shell)
579 natom = 0
580 ncore = 0
581 nshell = 0
582 numer_stress = 0.0_dp
583
584 logger => cp_get_default_logger()
585
586 dx_w = default_dx
587 IF (PRESENT(dx)) dx_w = dx
588 CALL force_env_get(force_env, subsys=subsys, globenv=globenv)
589 CALL cp_subsys_get(subsys, &
590 core_particles=core_particles, &
591 particles=particles, &
592 shell_particles=shell_particles, &
593 virial=virial)
594 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
595 extension=".stress_tensor")
596 IF (output_unit > 0) THEN
597 WRITE (output_unit, '(/A,A/)') ' **************************** ', &
598 'NUMERICAL STRESS ********************************'
599 END IF
600
601 ! Save all original particle positions
602 natom = particles%n_els
603 ALLOCATE (ref_pos_atom(natom, 3))
604 DO i = 1, natom
605 ref_pos_atom(i, :) = particles%els(i)%r
606 END DO
607 IF (ASSOCIATED(core_particles)) THEN
608 ncore = core_particles%n_els
609 ALLOCATE (ref_pos_core(ncore, 3))
610 DO i = 1, ncore
611 ref_pos_core(i, :) = core_particles%els(i)%r
612 END DO
613 END IF
614 IF (ASSOCIATED(shell_particles)) THEN
615 nshell = shell_particles%n_els
616 ALLOCATE (ref_pos_shell(nshell, 3))
617 DO i = 1, nshell
618 ref_pos_shell(i, :) = shell_particles%els(i)%r
619 END DO
620 END IF
621 CALL force_env_get(force_env, cell=cell)
622 ! Save cell symmetry (distorted cell has no symmetry)
623 symmetry_id = cell%symmetry_id
624 cell%symmetry_id = cell_sym_triclinic
625 !
626 CALL cell_create(cell_local)
627 CALL cell_clone(cell, cell_local)
628 ! First change box
629 DO ip = 1, 3
630 DO iq = 1, 3
631 IF (virial%pv_diagonal .AND. (ip /= iq)) cycle
632 DO k = 1, 2
633 cell%hmat(ip, iq) = cell_local%hmat(ip, iq) - (-1.0_dp)**k*dx_w
634 CALL init_cell(cell)
635 ! Scale positions
636 DO i = 1, natom
637 CALL real_to_scaled(s, ref_pos_atom(i, 1:3), cell_local)
638 CALL scaled_to_real(particles%els(i)%r, s, cell)
639 END DO
640 DO i = 1, ncore
641 CALL real_to_scaled(s, ref_pos_core(i, 1:3), cell_local)
642 CALL scaled_to_real(core_particles%els(i)%r, s, cell)
643 END DO
644 DO i = 1, nshell
645 CALL real_to_scaled(s, ref_pos_shell(i, 1:3), cell_local)
646 CALL scaled_to_real(shell_particles%els(i)%r, s, cell)
647 END DO
648 ! Compute energies
649 CALL force_env_calc_energy_force(force_env, &
650 calc_force=.false., &
651 consistent_energies=.true., &
652 calc_stress_tensor=.false.)
653 CALL force_env_get(force_env, potential_energy=numer_energy(k))
654 ! Reset cell
655 cell%hmat(ip, iq) = cell_local%hmat(ip, iq)
656 END DO
657 CALL init_cell(cell)
658 numer_stress(ip, iq) = 0.5_dp*(numer_energy(1) - numer_energy(2))/dx_w
659 IF (output_unit > 0) THEN
660 IF (globenv%run_type_id == debug_run) THEN
661 WRITE (unit=output_unit, fmt="(/,T2,A,T19,A,F7.4,A,T44,A,F7.4,A,T69,A)") &
662 "DEBUG|", "E("//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
663 "E("//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
664 "f(numerical)"
665 WRITE (unit=output_unit, fmt="(T2,A,2(1X,F24.8),1X,F22.8)") &
666 "DEBUG|", numer_energy(1:2), numer_stress(ip, iq)
667 ELSE
668 WRITE (unit=output_unit, fmt="(/,T7,A,F7.4,A,T27,A,F7.4,A,T49,A)") &
669 "E("//achar(119 + ip)//achar(119 + iq)//" +", dx_w, ")", &
670 "E("//achar(119 + ip)//achar(119 + iq)//" -", dx_w, ")", &
671 "f(numerical)"
672 WRITE (unit=output_unit, fmt="(3(1X,F19.8))") &
673 numer_energy(1:2), numer_stress(ip, iq)
674 END IF
675 END IF
676 END DO
677 END DO
678
679 ! Reset positions and rebuild original environment
680 cell%symmetry_id = symmetry_id
681 CALL init_cell(cell)
682 DO i = 1, natom
683 particles%els(i)%r = ref_pos_atom(i, :)
684 END DO
685 DO i = 1, ncore
686 core_particles%els(i)%r = ref_pos_core(i, :)
687 END DO
688 DO i = 1, nshell
689 shell_particles%els(i)%r = ref_pos_shell(i, :)
690 END DO
691 CALL force_env_calc_energy_force(force_env, &
692 calc_force=.false., &
693 consistent_energies=.true., &
694 calc_stress_tensor=.false.)
695
696 ! Computing pv_test
697 virial%pv_virial = 0.0_dp
698 DO i = 1, 3
699 DO j = 1, 3
700 DO k = 1, 3
701 virial%pv_virial(i, j) = virial%pv_virial(i, j) - &
702 0.5_dp*(numer_stress(i, k)*cell_local%hmat(j, k) + &
703 numer_stress(j, k)*cell_local%hmat(i, k))
704 END DO
705 END DO
706 END DO
707
708 IF (output_unit > 0) THEN
709 IF (globenv%run_type_id == debug_run) THEN
710 CALL write_stress_tensor(virial%pv_virial, output_unit, cell, virial%pv_numer)
711 END IF
712 WRITE (output_unit, '(/,A,/)') ' **************************** '// &
713 'NUMERICAL STRESS END *****************************'
714 END IF
715
716 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
717 "PRINT%STRESS_TENSOR")
718
719 ! Release storage
720 IF (ASSOCIATED(ref_pos_atom)) THEN
721 DEALLOCATE (ref_pos_atom)
722 END IF
723 IF (ASSOCIATED(ref_pos_core)) THEN
724 DEALLOCATE (ref_pos_core)
725 END IF
726 IF (ASSOCIATED(ref_pos_shell)) THEN
727 DEALLOCATE (ref_pos_shell)
728 END IF
729 IF (ASSOCIATED(cell_local)) CALL cell_release(cell_local)
730
731 END SUBROUTINE force_env_calc_num_pressure
732
733! **************************************************************************************************
734!> \brief creates and initializes a force environment
735!> \param force_env the force env to create
736!> \param root_section ...
737!> \param para_env ...
738!> \param globenv ...
739!> \param fist_env , qs_env: exactly one of these should be
740!> associated, the one that is active
741!> \param qs_env ...
742!> \param meta_env ...
743!> \param sub_force_env ...
744!> \param qmmm_env ...
745!> \param qmmmx_env ...
746!> \param eip_env ...
747!> \param pwdft_env ...
748!> \param force_env_section ...
749!> \param mixed_env ...
750!> \param embed_env ...
751!> \param nnp_env ...
752!> \par History
753!> 04.2003 created [fawzi]
754!> \author fawzi
755! **************************************************************************************************
756 SUBROUTINE force_env_create(force_env, root_section, para_env, globenv, fist_env, &
757 qs_env, meta_env, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, force_env_section, &
758 mixed_env, embed_env, nnp_env)
759
760 TYPE(force_env_type), POINTER :: force_env
761 TYPE(section_vals_type), POINTER :: root_section
762 TYPE(mp_para_env_type), POINTER :: para_env
763 TYPE(global_environment_type), POINTER :: globenv
764 TYPE(fist_environment_type), OPTIONAL, POINTER :: fist_env
765 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
766 TYPE(meta_env_type), OPTIONAL, POINTER :: meta_env
767 TYPE(force_env_p_type), DIMENSION(:), OPTIONAL, &
768 POINTER :: sub_force_env
769 TYPE(qmmm_env_type), OPTIONAL, POINTER :: qmmm_env
770 TYPE(qmmmx_env_type), OPTIONAL, POINTER :: qmmmx_env
771 TYPE(eip_environment_type), OPTIONAL, POINTER :: eip_env
772 TYPE(pwdft_environment_type), OPTIONAL, POINTER :: pwdft_env
773 TYPE(section_vals_type), POINTER :: force_env_section
774 TYPE(mixed_environment_type), OPTIONAL, POINTER :: mixed_env
775 TYPE(embed_env_type), OPTIONAL, POINTER :: embed_env
776 TYPE(nnp_type), OPTIONAL, POINTER :: nnp_env
777
778 ALLOCATE (force_env)
779 NULLIFY (force_env%fist_env, force_env%qs_env, &
780 force_env%para_env, force_env%globenv, &
781 force_env%meta_env, force_env%sub_force_env, &
782 force_env%qmmm_env, force_env%qmmmx_env, force_env%fp_env, &
783 force_env%force_env_section, force_env%eip_env, force_env%mixed_env, &
784 force_env%embed_env, force_env%pwdft_env, force_env%nnp_env, &
785 force_env%root_section)
786 last_force_env_id = last_force_env_id + 1
787 force_env%ref_count = 1
788 force_env%in_use = 0
789 force_env%additional_potential = 0.0_dp
790
791 force_env%globenv => globenv
792 CALL globenv_retain(force_env%globenv)
793
794 force_env%root_section => root_section
795 CALL section_vals_retain(root_section)
796
797 force_env%para_env => para_env
798 CALL force_env%para_env%retain()
799
800 CALL section_vals_retain(force_env_section)
801 force_env%force_env_section => force_env_section
802
803 IF (PRESENT(fist_env)) THEN
804 cpassert(ASSOCIATED(fist_env))
805 cpassert(force_env%in_use == 0)
806 force_env%in_use = use_fist_force
807 force_env%fist_env => fist_env
808 END IF
809 IF (PRESENT(eip_env)) THEN
810 cpassert(ASSOCIATED(eip_env))
811 cpassert(force_env%in_use == 0)
812 force_env%in_use = use_eip_force
813 force_env%eip_env => eip_env
814 END IF
815 IF (PRESENT(pwdft_env)) THEN
816 cpassert(ASSOCIATED(pwdft_env))
817 cpassert(force_env%in_use == 0)
818 force_env%in_use = use_pwdft_force
819 force_env%pwdft_env => pwdft_env
820 END IF
821 IF (PRESENT(qs_env)) THEN
822 cpassert(ASSOCIATED(qs_env))
823 cpassert(force_env%in_use == 0)
824 force_env%in_use = use_qs_force
825 force_env%qs_env => qs_env
826 END IF
827 IF (PRESENT(qmmm_env)) THEN
828 cpassert(ASSOCIATED(qmmm_env))
829 cpassert(force_env%in_use == 0)
830 force_env%in_use = use_qmmm
831 force_env%qmmm_env => qmmm_env
832 END IF
833 IF (PRESENT(qmmmx_env)) THEN
834 cpassert(ASSOCIATED(qmmmx_env))
835 cpassert(force_env%in_use == 0)
836 force_env%in_use = use_qmmmx
837 force_env%qmmmx_env => qmmmx_env
838 END IF
839 IF (PRESENT(mixed_env)) THEN
840 cpassert(ASSOCIATED(mixed_env))
841 cpassert(force_env%in_use == 0)
842 force_env%in_use = use_mixed_force
843 force_env%mixed_env => mixed_env
844 END IF
845 IF (PRESENT(embed_env)) THEN
846 cpassert(ASSOCIATED(embed_env))
847 cpassert(force_env%in_use == 0)
848 force_env%in_use = use_embed
849 force_env%embed_env => embed_env
850 END IF
851 IF (PRESENT(nnp_env)) THEN
852 cpassert(ASSOCIATED(nnp_env))
853 cpassert(force_env%in_use == 0)
854 force_env%in_use = use_nnp_force
855 force_env%nnp_env => nnp_env
856 END IF
857 cpassert(force_env%in_use /= 0)
858
859 IF (PRESENT(sub_force_env)) THEN
860 force_env%sub_force_env => sub_force_env
861 END IF
862
863 IF (PRESENT(meta_env)) THEN
864 force_env%meta_env => meta_env
865 ELSE
866 NULLIFY (force_env%meta_env)
867 END IF
868
869 END SUBROUTINE force_env_create
870
871! **************************************************************************************************
872!> \brief ****f* force_env_methods/mixed_energy_forces [1.0]
873!>
874!> Computes energy and forces for a mixed force_env type
875!> \param force_env the force_env that holds the mixed_env type
876!> \param calculate_forces decides if forces should be calculated
877!> \par History
878!> 11.06 created [fschiff]
879!> 04.07 generalization to an illimited number of force_eval [tlaino]
880!> 04.07 further generalization to force_eval with different geometrical
881!> structures [tlaino]
882!> 04.08 reorganizing the genmix structure (collecting common code)
883!> 01.16 added CDFT [Nico Holmberg]
884!> 08.17 added DFT embedding [Vladimir Rybkin]
885!> \author Florian Schiffmann
886! **************************************************************************************************
887 SUBROUTINE mixed_energy_forces(force_env, calculate_forces)
888
889 TYPE(force_env_type), POINTER :: force_env
890 LOGICAL, INTENT(IN) :: calculate_forces
891
892 CHARACTER(LEN=default_path_length) :: coupling_function
893 CHARACTER(LEN=default_string_length) :: def_error, description, this_error
894 INTEGER :: iforce_eval, iparticle, istate(2), &
895 jparticle, mixing_type, my_group, &
896 natom, nforce_eval, source, unit_nr
897 INTEGER, DIMENSION(:), POINTER :: glob_natoms, itmplist, map_index
898 LOGICAL :: dip_exists
899 REAL(kind=dp) :: coupling_parameter, dedf, der_1, der_2, &
900 dx, energy, err, lambda, lerr, &
901 restraint_strength, restraint_target, &
902 sd
903 REAL(kind=dp), DIMENSION(3) :: dip_mix
904 REAL(kind=dp), DIMENSION(:), POINTER :: energies
905 TYPE(cell_type), POINTER :: cell_mix
906 TYPE(cp_logger_type), POINTER :: logger, my_logger
907 TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
908 TYPE(cp_result_type), POINTER :: loc_results, results_mix
909 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
910 TYPE(cp_subsys_type), POINTER :: subsys_mix
911 TYPE(mixed_energy_type), POINTER :: mixed_energy
912 TYPE(mixed_force_type), DIMENSION(:), POINTER :: global_forces
913 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
914 TYPE(particle_list_type), POINTER :: particles_mix
915 TYPE(section_vals_type), POINTER :: force_env_section, gen_section, &
916 mapping_section, mixed_section, &
917 root_section
918 TYPE(virial_p_type), DIMENSION(:), POINTER :: virials
919 TYPE(virial_type), POINTER :: loc_virial, virial_mix
920
921 logger => cp_get_default_logger()
922 cpassert(ASSOCIATED(force_env))
923 ! Get infos about the mixed subsys
924 CALL force_env_get(force_env=force_env, &
925 subsys=subsys_mix, &
926 force_env_section=force_env_section, &
927 root_section=root_section, &
928 cell=cell_mix)
929 CALL cp_subsys_get(subsys=subsys_mix, &
930 particles=particles_mix, &
931 virial=virial_mix, &
932 results=results_mix)
933 NULLIFY (map_index, glob_natoms, global_forces, itmplist)
934
935 nforce_eval = SIZE(force_env%sub_force_env)
936 mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
937 mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
938 ! Global Info
939 ALLOCATE (subsystems(nforce_eval))
940 ALLOCATE (particles(nforce_eval))
941 ! Local Info to sync
942 ALLOCATE (global_forces(nforce_eval))
943 ALLOCATE (energies(nforce_eval))
944 ALLOCATE (glob_natoms(nforce_eval))
945 ALLOCATE (virials(nforce_eval))
946 ALLOCATE (results(nforce_eval))
947 energies = 0.0_dp
948 glob_natoms = 0
949 ! Check if mixed CDFT calculation is requested and initialize
950 CALL mixed_cdft_init(force_env, calculate_forces)
951
952 !
953 IF (.NOT. force_env%mixed_env%do_mixed_cdft) THEN
954 DO iforce_eval = 1, nforce_eval
955 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
956 NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
957 ALLOCATE (virials(iforce_eval)%virial)
958 CALL cp_result_create(results(iforce_eval)%results)
959 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
960 ! From this point on the error is the sub_error
961 my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
962 my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
963 ! Copy iterations info (they are updated only in the main mixed_env)
964 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
965 CALL cp_add_default_logger(my_logger)
966
967 ! Get all available subsys
968 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
969 subsys=subsystems(iforce_eval)%subsys)
970
971 ! all force_env share the same cell
972 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
973
974 ! Get available particles
975 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
976 particles=particles(iforce_eval)%list)
977
978 ! Get Mapping index array
979 natom = SIZE(particles(iforce_eval)%list%els)
980
981 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
982 map_index)
983
984 ! Mapping particles from iforce_eval environment to the mixed env
985 DO iparticle = 1, natom
986 jparticle = map_index(iparticle)
987 particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
988 END DO
989
990 ! Calculate energy and forces for each sub_force_env
991 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
992 calc_force=calculate_forces, &
993 skip_external_control=.true.)
994
995 ! Only the rank 0 process collect info for each computation
996 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
997 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
998 potential_energy=energy)
999 CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1000 virial=loc_virial, results=loc_results)
1001 energies(iforce_eval) = energy
1002 glob_natoms(iforce_eval) = natom
1003 virials(iforce_eval)%virial = loc_virial
1004 CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1005 END IF
1006 ! Deallocate map_index array
1007 IF (ASSOCIATED(map_index)) THEN
1008 DEALLOCATE (map_index)
1009 END IF
1010 CALL cp_rm_default_logger()
1011 END DO
1012 ELSE
1013 CALL mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
1014 glob_natoms, virials, results)
1015 END IF
1016 ! Handling Parallel execution
1017 CALL force_env%para_env%sync()
1018 ! Post CDFT operations
1019 CALL mixed_cdft_post_energy_forces(force_env)
1020 ! Let's transfer energy, natom, forces, virials
1021 CALL force_env%para_env%sum(energies)
1022 CALL force_env%para_env%sum(glob_natoms)
1023 ! Transfer forces
1024 DO iforce_eval = 1, nforce_eval
1025 ALLOCATE (global_forces(iforce_eval)%forces(3, glob_natoms(iforce_eval)))
1026 global_forces(iforce_eval)%forces = 0.0_dp
1027 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN
1028 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1029 ! Forces
1030 DO iparticle = 1, glob_natoms(iforce_eval)
1031 global_forces(iforce_eval)%forces(:, iparticle) = &
1032 particles(iforce_eval)%list%els(iparticle)%f
1033 END DO
1034 END IF
1035 END IF
1036 CALL force_env%para_env%sum(global_forces(iforce_eval)%forces)
1037 !Transfer only the relevant part of the virial..
1038 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_total)
1039 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_kinetic)
1040 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_virial)
1041 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_xc)
1042 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_fock_4c)
1043 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_constraint)
1044 !Transfer results
1045 source = 0
1046 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN
1047 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) &
1048 source = force_env%para_env%mepos
1049 END IF
1050 CALL force_env%para_env%sum(source)
1051 CALL cp_results_mp_bcast(results(iforce_eval)%results, source, force_env%para_env)
1052 END DO
1053
1054 force_env%mixed_env%energies = energies
1055 ! Start combining the different sub_force_env
1056 CALL get_mixed_env(mixed_env=force_env%mixed_env, &
1057 mixed_energy=mixed_energy)
1058
1059 !NB: do this for all MIXING_TYPE values, since some need it (e.g. linear mixing
1060 !NB if the first system has fewer atoms than the second)
1061 DO iparticle = 1, SIZE(particles_mix%els)
1062 particles_mix%els(iparticle)%f(:) = 0.0_dp
1063 END DO
1064
1065 CALL section_vals_val_get(mixed_section, "MIXING_TYPE", i_val=mixing_type)
1066 SELECT CASE (mixing_type)
1067 CASE (mix_linear_combination)
1068 ! Support offered only 2 force_eval
1069 cpassert(nforce_eval == 2)
1070 CALL section_vals_val_get(mixed_section, "LINEAR%LAMBDA", r_val=lambda)
1071 mixed_energy%pot = lambda*energies(1) + (1.0_dp - lambda)*energies(2)
1072 ! General Mapping of forces...
1073 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1074 lambda, 1, nforce_eval, map_index, mapping_section, .true.)
1075 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1076 (1.0_dp - lambda), 2, nforce_eval, map_index, mapping_section, .false.)
1077 CASE (mix_minimum)
1078 ! Support offered only 2 force_eval
1079 cpassert(nforce_eval == 2)
1080 IF (energies(1) < energies(2)) THEN
1081 mixed_energy%pot = energies(1)
1082 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1083 1.0_dp, 1, nforce_eval, map_index, mapping_section, .true.)
1084 ELSE
1085 mixed_energy%pot = energies(2)
1086 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1087 1.0_dp, 2, nforce_eval, map_index, mapping_section, .true.)
1088 END IF
1089 CASE (mix_coupled)
1090 ! Support offered only 2 force_eval
1091 cpassert(nforce_eval == 2)
1092 CALL section_vals_val_get(mixed_section, "COUPLING%COUPLING_PARAMETER", &
1093 r_val=coupling_parameter)
1094 sd = sqrt((energies(1) - energies(2))**2 + 4.0_dp*coupling_parameter**2)
1095 der_1 = (1.0_dp - (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1096 der_2 = (1.0_dp + (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1097 mixed_energy%pot = (energies(1) + energies(2) - sd)/2.0_dp
1098 ! General Mapping of forces...
1099 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1100 der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1101 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1102 der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1103 CASE (mix_restrained)
1104 ! Support offered only 2 force_eval
1105 cpassert(nforce_eval == 2)
1106 CALL section_vals_val_get(mixed_section, "RESTRAINT%RESTRAINT_TARGET", &
1107 r_val=restraint_target)
1108 CALL section_vals_val_get(mixed_section, "RESTRAINT%RESTRAINT_STRENGTH", &
1109 r_val=restraint_strength)
1110 mixed_energy%pot = energies(1) + restraint_strength*(energies(1) - energies(2) - restraint_target)**2
1111 der_2 = -2.0_dp*restraint_strength*(energies(1) - energies(2) - restraint_target)
1112 der_1 = 1.0_dp - der_2
1113 ! General Mapping of forces...
1114 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1115 der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1116 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1117 der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1118 CASE (mix_generic)
1119 ! Support any number of force_eval sections
1120 gen_section => section_vals_get_subs_vals(mixed_section, "GENERIC")
1121 CALL get_generic_info(gen_section, "MIXING_FUNCTION", coupling_function, force_env%mixed_env%par, &
1122 force_env%mixed_env%val, energies)
1123 CALL initf(1)
1124 CALL parsef(1, trim(coupling_function), force_env%mixed_env%par)
1125 ! Now the hardest part.. map energy with corresponding force_eval
1126 mixed_energy%pot = evalf(1, force_env%mixed_env%val)
1127 cpassert(evalerrtype <= 0)
1128 CALL zero_virial(virial_mix, reset=.false.)
1129 CALL cp_results_erase(results_mix)
1130 DO iforce_eval = 1, nforce_eval
1131 CALL section_vals_val_get(gen_section, "DX", r_val=dx)
1132 CALL section_vals_val_get(gen_section, "ERROR_LIMIT", r_val=lerr)
1133 dedf = evalfd(1, iforce_eval, force_env%mixed_env%val, dx, err)
1134 IF (abs(err) > lerr) THEN
1135 WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
1136 WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
1137 CALL compress(this_error, .true.)
1138 CALL compress(def_error, .true.)
1139 CALL cp_warn(__location__, &
1140 'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
1141 ' Error '//trim(this_error)//' in computing numerical derivatives larger then'// &
1142 trim(def_error)//' .')
1143 END IF
1144 ! General Mapping of forces...
1145 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1146 dedf, iforce_eval, nforce_eval, map_index, mapping_section, .false.)
1147 force_env%mixed_env%val(iforce_eval) = energies(iforce_eval)
1148 END DO
1149 ! Let's store the needed information..
1150 force_env%mixed_env%dx = dx
1151 force_env%mixed_env%lerr = lerr
1152 force_env%mixed_env%coupling_function = trim(coupling_function)
1153 CALL finalizef()
1154 CASE (mix_cdft)
1155 ! Supports any number of force_evals for calculation of CDFT properties, but forces only from two
1156 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LAMBDA", r_val=lambda)
1157 ! Get the states which determine the forces
1158 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%FORCE_STATES", i_vals=itmplist)
1159 IF (SIZE(itmplist) /= 2) &
1160 CALL cp_abort(__location__, &
1161 "Keyword FORCE_STATES takes exactly two input values.")
1162 IF (any(itmplist .LT. 0)) &
1163 cpabort("Invalid force_eval index.")
1164 istate = itmplist
1165 IF (istate(1) > nforce_eval .OR. istate(2) > nforce_eval) &
1166 cpabort("Invalid force_eval index.")
1167 mixed_energy%pot = lambda*energies(istate(1)) + (1.0_dp - lambda)*energies(istate(2))
1168 ! General Mapping of forces...
1169 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1170 lambda, istate(1), nforce_eval, map_index, mapping_section, .true.)
1171 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1172 (1.0_dp - lambda), istate(2), nforce_eval, map_index, mapping_section, .false.)
1173 CASE DEFAULT
1174 cpabort("")
1175 END SELECT
1176 !Simply deallocate and loose the pointer references..
1177 DO iforce_eval = 1, nforce_eval
1178 DEALLOCATE (global_forces(iforce_eval)%forces)
1179 IF (ASSOCIATED(virials(iforce_eval)%virial)) DEALLOCATE (virials(iforce_eval)%virial)
1180 CALL cp_result_release(results(iforce_eval)%results)
1181 END DO
1182 DEALLOCATE (global_forces)
1183 DEALLOCATE (subsystems)
1184 DEALLOCATE (particles)
1185 DEALLOCATE (energies)
1186 DEALLOCATE (glob_natoms)
1187 DEALLOCATE (virials)
1188 DEALLOCATE (results)
1189 ! Print Section
1190 unit_nr = cp_print_key_unit_nr(logger, mixed_section, "PRINT%DIPOLE", &
1191 extension=".data", middle_name="MIXED_DIPOLE", log_filename=.false.)
1192 IF (unit_nr > 0) THEN
1193 description = '[DIPOLE]'
1194 dip_exists = test_for_result(results=results_mix, description=description)
1195 IF (dip_exists) THEN
1196 CALL get_results(results=results_mix, description=description, values=dip_mix)
1197 WRITE (unit_nr, '(/,1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE ( A.U.)|", dip_mix
1198 WRITE (unit_nr, '( 1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE (Debye)|", dip_mix*debye
1199 ELSE
1200 WRITE (unit_nr, *) "NO FORCE_EVAL section calculated the dipole"
1201 END IF
1202 END IF
1203 CALL cp_print_key_finished_output(unit_nr, logger, mixed_section, "PRINT%DIPOLE")
1204 END SUBROUTINE mixed_energy_forces
1205
1206! **************************************************************************************************
1207!> \brief Driver routine for mixed CDFT energy and force calculations
1208!> \param force_env the force_env that holds the mixed_env
1209!> \param calculate_forces if forces should be calculated
1210!> \param particles system particles
1211!> \param energies the energies of the CDFT states
1212!> \param glob_natoms the total number of particles
1213!> \param virials the virials stored in subsys
1214!> \param results results stored in subsys
1215!> \par History
1216!> 01.17 created [Nico Holmberg]
1217!> \author Nico Holmberg
1218! **************************************************************************************************
1219 SUBROUTINE mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
1220 glob_natoms, virials, results)
1221 TYPE(force_env_type), POINTER :: force_env
1222 LOGICAL, INTENT(IN) :: calculate_forces
1223 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1224 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1225 INTEGER, DIMENSION(:), POINTER :: glob_natoms
1226 TYPE(virial_p_type), DIMENSION(:), POINTER :: virials
1227 TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1228
1229 INTEGER :: iforce_eval, iparticle, jparticle, &
1230 my_group, natom, nforce_eval
1231 INTEGER, DIMENSION(:), POINTER :: map_index
1232 REAL(kind=dp) :: energy
1233 TYPE(cell_type), POINTER :: cell_mix
1234 TYPE(cp_logger_type), POINTER :: logger, my_logger
1235 TYPE(cp_result_type), POINTER :: loc_results, results_mix
1236 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1237 TYPE(cp_subsys_type), POINTER :: subsys_mix
1238 TYPE(particle_list_type), POINTER :: particles_mix
1239 TYPE(section_vals_type), POINTER :: force_env_section, mapping_section, &
1240 mixed_section, root_section
1241 TYPE(virial_type), POINTER :: loc_virial, virial_mix
1242
1243 logger => cp_get_default_logger()
1244 cpassert(ASSOCIATED(force_env))
1245 ! Get infos about the mixed subsys
1246 CALL force_env_get(force_env=force_env, &
1247 subsys=subsys_mix, &
1248 force_env_section=force_env_section, &
1249 root_section=root_section, &
1250 cell=cell_mix)
1251 CALL cp_subsys_get(subsys=subsys_mix, &
1252 particles=particles_mix, &
1253 virial=virial_mix, &
1254 results=results_mix)
1255 NULLIFY (map_index)
1256 nforce_eval = SIZE(force_env%sub_force_env)
1257 mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
1258 mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
1259 ALLOCATE (subsystems(nforce_eval))
1260 DO iforce_eval = 1, nforce_eval
1261 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1262 NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
1263 ALLOCATE (virials(iforce_eval)%virial)
1264 CALL cp_result_create(results(iforce_eval)%results)
1265 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1266 ! Get all available subsys
1267 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1268 subsys=subsystems(iforce_eval)%subsys)
1269
1270 ! all force_env share the same cell
1271 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
1272
1273 ! Get available particles
1274 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1275 particles=particles(iforce_eval)%list)
1276
1277 ! Get Mapping index array
1278 natom = SIZE(particles(iforce_eval)%list%els)
1279 ! Serial mode need to deallocate first
1280 IF (ASSOCIATED(map_index)) &
1281 DEALLOCATE (map_index)
1282 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1283 map_index)
1284
1285 ! Mapping particles from iforce_eval environment to the mixed env
1286 DO iparticle = 1, natom
1287 jparticle = map_index(iparticle)
1288 particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
1289 END DO
1290 ! Mixed CDFT + QMMM: Need to translate now
1291 IF (force_env%mixed_env%do_mixed_qmmm_cdft) &
1292 CALL apply_qmmm_translate(force_env%sub_force_env(iforce_eval)%force_env%qmmm_env)
1293 END DO
1294 ! For mixed CDFT calculations parallelized over CDFT states
1295 ! build weight and gradient on all processors before splitting into groups and
1296 ! starting energy calculation
1297 CALL mixed_cdft_build_weight(force_env, calculate_forces)
1298 DO iforce_eval = 1, nforce_eval
1299 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1300 ! From this point on the error is the sub_error
1301 IF (force_env%mixed_env%cdft_control%run_type == mixed_cdft_serial .AND. iforce_eval .GE. 2) THEN
1302 my_logger => force_env%mixed_env%cdft_control%sub_logger(iforce_eval - 1)%p
1303 ELSE
1304 my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
1305 my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
1306 END IF
1307 ! Copy iterations info (they are updated only in the main mixed_env)
1308 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1309 CALL cp_add_default_logger(my_logger)
1310 ! Serial CDFT calculation: transfer weight/gradient
1311 CALL mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
1312 ! Calculate energy and forces for each sub_force_env
1313 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1314 calc_force=calculate_forces, &
1315 skip_external_control=.true.)
1316 ! Only the rank 0 process collect info for each computation
1317 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1318 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1319 potential_energy=energy)
1320 CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1321 virial=loc_virial, results=loc_results)
1322 energies(iforce_eval) = energy
1323 glob_natoms(iforce_eval) = natom
1324 virials(iforce_eval)%virial = loc_virial
1325 CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1326 END IF
1327 ! Deallocate map_index array
1328 IF (ASSOCIATED(map_index)) THEN
1329 DEALLOCATE (map_index)
1330 END IF
1331 CALL cp_rm_default_logger()
1332 END DO
1333 DEALLOCATE (subsystems)
1334
1335 END SUBROUTINE mixed_cdft_energy_forces
1336
1337! **************************************************************************************************
1338!> \brief Perform additional tasks for mixed CDFT calculations after solving the electronic structure
1339!> of both CDFT states
1340!> \param force_env the force_env that holds the CDFT states
1341!> \par History
1342!> 01.17 created [Nico Holmberg]
1343!> \author Nico Holmberg
1344! **************************************************************************************************
1345 SUBROUTINE mixed_cdft_post_energy_forces(force_env)
1346 TYPE(force_env_type), POINTER :: force_env
1347
1348 INTEGER :: iforce_eval, nforce_eval, nvar
1349 TYPE(dft_control_type), POINTER :: dft_control
1350 TYPE(qs_environment_type), POINTER :: qs_env
1351
1352 cpassert(ASSOCIATED(force_env))
1353 NULLIFY (qs_env, dft_control)
1354 IF (force_env%mixed_env%do_mixed_cdft) THEN
1355 nforce_eval = SIZE(force_env%sub_force_env)
1356 nvar = force_env%mixed_env%cdft_control%nconstraint
1357 ! Transfer cdft strengths for writing restart
1358 IF (.NOT. ASSOCIATED(force_env%mixed_env%strength)) &
1359 ALLOCATE (force_env%mixed_env%strength(nforce_eval, nvar))
1360 force_env%mixed_env%strength = 0.0_dp
1361 DO iforce_eval = 1, nforce_eval
1362 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1363 IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1364 qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
1365 ELSE
1366 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1367 END IF
1368 CALL get_qs_env(qs_env, dft_control=dft_control)
1369 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) &
1370 force_env%mixed_env%strength(iforce_eval, :) = dft_control%qs_control%cdft_control%strength(:)
1371 END DO
1372 CALL force_env%para_env%sum(force_env%mixed_env%strength)
1373 ! Mixed CDFT: calculate ET coupling
1374 IF (force_env%mixed_env%do_mixed_et) THEN
1375 IF (modulo(force_env%mixed_env%cdft_control%sim_step, force_env%mixed_env%et_freq) == 0) &
1376 CALL mixed_cdft_calculate_coupling(force_env)
1377 END IF
1378 END IF
1379
1380 END SUBROUTINE mixed_cdft_post_energy_forces
1381
1382! **************************************************************************************************
1383!> \brief Computes the total energy for an embedded calculation
1384!> \param force_env ...
1385!> \author Vladimir Rybkin
1386! **************************************************************************************************
1387 SUBROUTINE embed_energy(force_env)
1388
1389 TYPE(force_env_type), POINTER :: force_env
1390
1391 INTEGER :: iforce_eval, iparticle, jparticle, &
1392 my_group, natom, nforce_eval
1393 INTEGER, DIMENSION(:), POINTER :: glob_natoms, map_index
1394 LOGICAL :: converged_embed
1395 REAL(kind=dp) :: energy
1396 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1397 TYPE(cell_type), POINTER :: cell_embed
1398 TYPE(cp_logger_type), POINTER :: logger, my_logger
1399 TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
1400 TYPE(cp_result_type), POINTER :: loc_results, results_embed
1401 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
1402 TYPE(cp_subsys_type), POINTER :: subsys_embed
1403 TYPE(dft_control_type), POINTER :: dft_control
1404 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1405 TYPE(particle_list_type), POINTER :: particles_embed
1406 TYPE(pw_env_type), POINTER :: pw_env
1407 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1408 TYPE(pw_r3d_rs_type), POINTER :: embed_pot, spin_embed_pot
1409 TYPE(section_vals_type), POINTER :: embed_section, force_env_section, &
1410 mapping_section, root_section
1411
1412 logger => cp_get_default_logger()
1413 cpassert(ASSOCIATED(force_env))
1414 ! Get infos about the embedding subsys
1415 CALL force_env_get(force_env=force_env, &
1416 subsys=subsys_embed, &
1417 force_env_section=force_env_section, &
1418 root_section=root_section, &
1419 cell=cell_embed)
1420 CALL cp_subsys_get(subsys=subsys_embed, &
1421 particles=particles_embed, &
1422 results=results_embed)
1423 NULLIFY (map_index, glob_natoms)
1424
1425 nforce_eval = SIZE(force_env%sub_force_env)
1426 embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1427 mapping_section => section_vals_get_subs_vals(embed_section, "MAPPING")
1428 ! Global Info
1429 ALLOCATE (subsystems(nforce_eval))
1430 ALLOCATE (particles(nforce_eval))
1431 ! Local Info to sync
1432 ALLOCATE (energies(nforce_eval))
1433 ALLOCATE (glob_natoms(nforce_eval))
1434 ALLOCATE (results(nforce_eval))
1435 energies = 0.0_dp
1436 glob_natoms = 0
1437
1438 DO iforce_eval = 1, nforce_eval
1439 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1440 NULLIFY (results(iforce_eval)%results)
1441 CALL cp_result_create(results(iforce_eval)%results)
1442 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1443 ! From this point on the error is the sub_error
1444 my_group = force_env%embed_env%group_distribution(force_env%para_env%mepos)
1445 my_logger => force_env%embed_env%sub_logger(my_group + 1)%p
1446 ! Copy iterations info (they are updated only in the main embed_env)
1447 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
1448 CALL cp_add_default_logger(my_logger)
1449
1450 ! Get all available subsys
1451 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1452 subsys=subsystems(iforce_eval)%subsys)
1453
1454 ! Check if we import density from previous force calculations
1455 ! Only for QUICKSTEP
1456 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env%qs_env)) THEN
1457 NULLIFY (dft_control)
1458 CALL get_qs_env(force_env%sub_force_env(iforce_eval)%force_env%qs_env, dft_control=dft_control)
1459 IF (dft_control%qs_control%ref_embed_subsys) THEN
1460 IF (iforce_eval .EQ. 2) cpabort("Density importing force_eval can't be the first.")
1461 END IF
1462 END IF
1463
1464 ! all force_env share the same cell
1465 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_embed)
1466
1467 ! Get available particles
1468 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
1469 particles=particles(iforce_eval)%list)
1470
1471 ! Get Mapping index array
1472 natom = SIZE(particles(iforce_eval)%list%els)
1473
1474 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
1475 map_index, .true.)
1476
1477 ! Mapping particles from iforce_eval environment to the embed env
1478 DO iparticle = 1, natom
1479 jparticle = map_index(iparticle)
1480 particles(iforce_eval)%list%els(iparticle)%r = particles_embed%els(jparticle)%r
1481 END DO
1482
1483 ! Calculate energy and forces for each sub_force_env
1484 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
1485 calc_force=.false., &
1486 skip_external_control=.true.)
1487
1488 ! Call DFT embedding
1489 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env%qs_env)) THEN
1490 NULLIFY (dft_control)
1491 CALL get_qs_env(force_env%sub_force_env(iforce_eval)%force_env%qs_env, dft_control=dft_control)
1492 IF (dft_control%qs_control%ref_embed_subsys) THEN
1493 ! Now we can optimize the embedding potential
1494 CALL dft_embedding(force_env, iforce_eval, energies, converged_embed)
1495 IF (.NOT. converged_embed) cpabort("Embedding potential optimization not converged.")
1496 END IF
1497 ! Deallocate embedding potential on the high-level subsystem
1498 IF (dft_control%qs_control%high_level_embed_subsys) THEN
1499 CALL get_qs_env(qs_env=force_env%sub_force_env(iforce_eval)%force_env%qs_env, &
1500 embed_pot=embed_pot, spin_embed_pot=spin_embed_pot, pw_env=pw_env)
1501 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1502 CALL auxbas_pw_pool%give_back_pw(embed_pot)
1503 IF (ASSOCIATED(embed_pot)) THEN
1504 CALL embed_pot%release()
1505 DEALLOCATE (embed_pot)
1506 END IF
1507 IF (ASSOCIATED(spin_embed_pot)) THEN
1508 CALL auxbas_pw_pool%give_back_pw(spin_embed_pot)
1509 CALL spin_embed_pot%release()
1510 DEALLOCATE (spin_embed_pot)
1511 END IF
1512 END IF
1513 END IF
1514
1515 ! Only the rank 0 process collect info for each computation
1516 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1517 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1518 potential_energy=energy)
1519 CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
1520 results=loc_results)
1521 energies(iforce_eval) = energy
1522 glob_natoms(iforce_eval) = natom
1523 CALL cp_result_copy(loc_results, results(iforce_eval)%results)
1524 END IF
1525 ! Deallocate map_index array
1526 IF (ASSOCIATED(map_index)) THEN
1527 DEALLOCATE (map_index)
1528 END IF
1529 CALL cp_rm_default_logger()
1530 END DO
1531
1532 ! Handling Parallel execution
1533 CALL force_env%para_env%sync()
1534 ! Let's transfer energy, natom
1535 CALL force_env%para_env%sum(energies)
1536 CALL force_env%para_env%sum(glob_natoms)
1537
1538 force_env%embed_env%energies = energies
1539
1540 !NB if the first system has fewer atoms than the second)
1541 DO iparticle = 1, SIZE(particles_embed%els)
1542 particles_embed%els(iparticle)%f(:) = 0.0_dp
1543 END DO
1544
1545 ! ONIOM type of mixing in embedding: E = E_total + E_cluster_high - E_cluster
1546 force_env%embed_env%pot_energy = energies(3) + energies(4) - energies(2)
1547
1548 !Simply deallocate and loose the pointer references..
1549 DO iforce_eval = 1, nforce_eval
1550 CALL cp_result_release(results(iforce_eval)%results)
1551 END DO
1552 DEALLOCATE (subsystems)
1553 DEALLOCATE (particles)
1554 DEALLOCATE (energies)
1555 DEALLOCATE (glob_natoms)
1556 DEALLOCATE (results)
1557
1558 END SUBROUTINE embed_energy
1559
1560! **************************************************************************************************
1561!> \brief ...
1562!> \param force_env ...
1563!> \param ref_subsys_number ...
1564!> \param energies ...
1565!> \param converged_embed ...
1566! **************************************************************************************************
1567 SUBROUTINE dft_embedding(force_env, ref_subsys_number, energies, converged_embed)
1568 TYPE(force_env_type), POINTER :: force_env
1569 INTEGER :: ref_subsys_number
1570 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1571 LOGICAL :: converged_embed
1572
1573 INTEGER :: embed_method
1574 TYPE(section_vals_type), POINTER :: embed_section, force_env_section
1575
1576 ! Find out which embedding scheme is used
1577 CALL force_env_get(force_env=force_env, &
1578 force_env_section=force_env_section)
1579 embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1580
1581 CALL section_vals_val_get(embed_section, "EMBED_METHOD", i_val=embed_method)
1582 SELECT CASE (embed_method)
1583 CASE (dfet)
1584 ! Density functional embedding
1585 CALL dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1586 CASE (dmfet)
1587 ! Density matrix embedding theory
1588 CALL dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1589 END SELECT
1590
1591 END SUBROUTINE dft_embedding
1592! **************************************************************************************************
1593!> \brief ... Main driver for DFT embedding
1594!> \param force_env ...
1595!> \param ref_subsys_number ...
1596!> \param energies ...
1597!> \param converged_embed ...
1598!> \author Vladimir Rybkin
1599! **************************************************************************************************
1600 SUBROUTINE dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1601 TYPE(force_env_type), POINTER :: force_env
1602 INTEGER :: ref_subsys_number
1603 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1604 LOGICAL :: converged_embed
1605
1606 CHARACTER(LEN=*), PARAMETER :: routinen = 'dfet_embedding'
1607
1608 INTEGER :: cluster_subsys_num, handle, &
1609 i_force_eval, i_iter, i_spin, &
1610 nforce_eval, nspins, nspins_subsys, &
1611 output_unit
1612 REAL(kind=dp) :: cluster_energy
1613 REAL(kind=dp), DIMENSION(:), POINTER :: rhs
1614 TYPE(cp_logger_type), POINTER :: logger
1615 TYPE(dft_control_type), POINTER :: dft_control
1616 TYPE(opt_embed_pot_type) :: opt_embed
1617 TYPE(pw_env_type), POINTER :: pw_env
1618 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1619 TYPE(pw_r3d_rs_type) :: diff_rho_r, diff_rho_spin
1620 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_ref, rho_r_subsys
1621 TYPE(pw_r3d_rs_type), POINTER :: embed_pot, embed_pot_subsys, &
1622 spin_embed_pot, spin_embed_pot_subsys
1623 TYPE(qs_energy_type), POINTER :: energy
1624 TYPE(qs_rho_type), POINTER :: rho, subsys_rho
1625 TYPE(section_vals_type), POINTER :: dft_section, embed_section, &
1626 force_env_section, input, &
1627 mapping_section, opt_embed_section
1628
1629 CALL timeset(routinen, handle)
1630
1631 CALL cite_reference(huang2011)
1632 CALL cite_reference(heaton_burgess2007)
1633
1634 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
1635
1636 ! Reveal input file
1637 NULLIFY (logger)
1638 logger => cp_get_default_logger()
1639 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
1640 extension=".Log")
1641
1642 NULLIFY (dft_section, input, opt_embed_section)
1643 NULLIFY (energy, dft_control)
1644
1645 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1646 pw_env=pw_env, dft_control=dft_control, rho=rho, energy=energy, &
1647 input=input)
1648 nspins = dft_control%nspins
1649
1650 dft_section => section_vals_get_subs_vals(input, "DFT")
1651 opt_embed_section => section_vals_get_subs_vals(input, &
1652 "DFT%QS%OPT_EMBED")
1653 ! Rho_r is the reference
1654 CALL qs_rho_get(rho_struct=rho, rho_r=rho_r_ref)
1655
1656 ! We need to understand how to treat spins states
1657 CALL understand_spin_states(force_env, ref_subsys_number, opt_embed%change_spin, opt_embed%open_shell_embed, &
1658 opt_embed%all_nspins)
1659
1660 ! Prepare everything for the potential maximization
1661 CALL prepare_embed_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, opt_embed, &
1662 opt_embed_section)
1663
1664 ! Initialize embedding potential
1665 CALL init_embed_pot(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, &
1666 opt_embed%add_const_pot, opt_embed%Fermi_Amaldi, opt_embed%const_pot, &
1667 opt_embed%open_shell_embed, spin_embed_pot, &
1668 opt_embed%pot_diff, opt_embed%Coulomb_guess, opt_embed%grid_opt)
1669
1670 ! Read embedding potential vector from the file
1671 IF (opt_embed%read_embed_pot .OR. opt_embed%read_embed_pot_cube) CALL read_embed_pot( &
1672 force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, spin_embed_pot, &
1673 opt_embed_section, opt_embed)
1674
1675 ! Prepare the pw object to store density differences
1676 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1677 CALL auxbas_pw_pool%create_pw(diff_rho_r)
1678 CALL pw_zero(diff_rho_r)
1679 IF (opt_embed%open_shell_embed) THEN
1680 CALL auxbas_pw_pool%create_pw(diff_rho_spin)
1681 CALL pw_zero(diff_rho_spin)
1682 END IF
1683
1684 ! Check the preliminary density differences
1685 DO i_spin = 1, nspins
1686 CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
1687 END DO
1688 IF (opt_embed%open_shell_embed) THEN ! Spin part
1689 IF (nspins .EQ. 2) THEN ! Reference systems has an open shell, else the reference diff_rho_spin is zero
1690 CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
1691 CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
1692 END IF
1693 END IF
1694
1695 DO i_force_eval = 1, ref_subsys_number - 1
1696 NULLIFY (subsys_rho, rho_r_subsys, dft_control)
1697 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, energy=energy, &
1698 dft_control=dft_control)
1699 nspins_subsys = dft_control%nspins
1700 ! Add subsystem densities
1701 CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
1702 DO i_spin = 1, nspins_subsys
1703 CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
1704 END DO
1705 IF (opt_embed%open_shell_embed) THEN ! Spin part
1706 IF (nspins_subsys .EQ. 2) THEN ! The subsystem makes contribution if it is spin-polarized
1707 ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
1708 IF ((i_force_eval .EQ. 2) .AND. (opt_embed%change_spin)) THEN
1709 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1710 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1711 ELSE
1712 ! First subsystem (always) and second subsystem (without spin change)
1713 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1714 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1715 END IF
1716 END IF
1717 END IF
1718 END DO
1719
1720 ! Print density difference
1721 CALL print_rho_diff(diff_rho_r, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1722 IF (opt_embed%open_shell_embed) THEN ! Spin part
1723 CALL print_rho_spin_diff(diff_rho_spin, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1724 END IF
1725
1726 ! Construct electrostatic guess if needed
1727 IF (opt_embed%Coulomb_guess) THEN
1728 ! Reveal resp charges for total system
1729 nforce_eval = SIZE(force_env%sub_force_env)
1730 NULLIFY (rhs)
1731 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, rhs=rhs)
1732 ! Get the mapping
1733 CALL force_env_get(force_env=force_env, &
1734 force_env_section=force_env_section)
1735 embed_section => section_vals_get_subs_vals(force_env_section, "EMBED")
1736 mapping_section => section_vals_get_subs_vals(embed_section, "MAPPING")
1737
1738 DO i_force_eval = 1, ref_subsys_number - 1
1739 IF (i_force_eval .EQ. 1) THEN
1740 CALL coulomb_guess(embed_pot, rhs, mapping_section, &
1741 force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1742 ELSE
1743 CALL coulomb_guess(opt_embed%pot_diff, rhs, mapping_section, &
1744 force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1745 END IF
1746 END DO
1747 CALL pw_axpy(opt_embed%pot_diff, embed_pot)
1748 IF (.NOT. opt_embed%grid_opt) CALL pw_copy(embed_pot, opt_embed%const_pot)
1749
1750 END IF
1751
1752 ! Difference guess
1753 IF (opt_embed%diff_guess) THEN
1754 CALL pw_copy(diff_rho_r, embed_pot)
1755 IF (.NOT. opt_embed%grid_opt) CALL pw_copy(embed_pot, opt_embed%const_pot)
1756 ! Open shell
1757 IF (opt_embed%open_shell_embed) CALL pw_copy(diff_rho_spin, spin_embed_pot)
1758 END IF
1759
1760 ! Calculate subsystems with trial embedding potential
1761 DO i_iter = 1, opt_embed%n_iter
1762 opt_embed%i_iter = i_iter
1763
1764 ! Set the density difference as the negative reference one
1765 CALL pw_zero(diff_rho_r)
1766 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, dft_control=dft_control)
1767 nspins = dft_control%nspins
1768 DO i_spin = 1, nspins
1769 CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
1770 END DO
1771 IF (opt_embed%open_shell_embed) THEN ! Spin part
1772 CALL pw_zero(diff_rho_spin)
1773 IF (nspins .EQ. 2) THEN ! Reference systems has an open shell, else the reference diff_rho_spin is zero
1774 CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
1775 CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
1776 END IF
1777 END IF
1778
1779 DO i_force_eval = 1, ref_subsys_number - 1
1780 NULLIFY (dft_control)
1781 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
1782 nspins_subsys = dft_control%nspins
1783
1784 IF ((i_force_eval .EQ. 2) .AND. (opt_embed%change_spin)) THEN
1785 ! Here we change the sign of the spin embedding potential due to spin change:
1786 ! only in spin_embed_subsys
1787 CALL make_subsys_embed_pot(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
1788 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
1789 opt_embed%open_shell_embed, .true.)
1790 ELSE ! Regular case
1791 CALL make_subsys_embed_pot(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
1792 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
1793 opt_embed%open_shell_embed, .false.)
1794 END IF
1795
1796 ! Switch on external potential in the subsystems
1797 dft_control%apply_embed_pot = .true.
1798
1799 ! Add the embedding potential
1800 CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, embed_pot=embed_pot_subsys)
1801 IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys .EQ. 2)) THEN ! Spin part
1802 CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
1803 spin_embed_pot=spin_embed_pot_subsys)
1804 END IF
1805
1806 ! Get the previous subsystem densities
1807 CALL get_prev_density(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
1808
1809 ! Calculate the new density
1810 CALL force_env_calc_energy_force(force_env=force_env%sub_force_env(i_force_eval)%force_env, &
1811 calc_force=.false., &
1812 skip_external_control=.true.)
1813
1814 CALL get_max_subsys_diff(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
1815
1816 ! Extract subsystem density and energy
1817 NULLIFY (rho_r_subsys, energy)
1818
1819 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, &
1820 energy=energy)
1821 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + energy%total
1822
1823 ! Find out which subsystem is the cluster
1824 IF (dft_control%qs_control%cluster_embed_subsys) THEN
1825 cluster_subsys_num = i_force_eval
1826 cluster_energy = energy%total
1827 END IF
1828
1829 ! Add subsystem densities
1830 CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
1831 DO i_spin = 1, nspins_subsys
1832 CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
1833 END DO
1834 IF (opt_embed%open_shell_embed) THEN ! Spin part
1835 IF (nspins_subsys .EQ. 2) THEN ! The subsystem makes contribution if it is spin-polarized
1836 ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
1837 IF ((i_force_eval .EQ. 2) .AND. (opt_embed%change_spin)) THEN
1838 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1839 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1840 ELSE
1841 ! First subsystem (always) and second subsystem (without spin change)
1842 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1843 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1844 END IF
1845 END IF
1846 END IF
1847
1848 ! Release embedding potential for subsystem
1849 CALL embed_pot_subsys%release()
1850 DEALLOCATE (embed_pot_subsys)
1851 IF (opt_embed%open_shell_embed) THEN
1852 CALL spin_embed_pot_subsys%release()
1853 DEALLOCATE (spin_embed_pot_subsys)
1854 END IF
1855
1856 END DO ! i_force_eval
1857
1858 ! Print embedding potential for restart
1859 CALL print_embed_restart(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1860 opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
1861 spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .false.)
1862 CALL print_pot_simple_grid(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1863 embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .false., &
1864 force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
1865
1866 ! Integrate the potential over density differences and add to w functional; also add regularization contribution
1867 DO i_spin = 1, nspins ! Sum over nspins for the reference system, not subsystem!
1868 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) - pw_integral_ab(embed_pot, rho_r_ref(i_spin))
1869 END DO
1870 ! Spin part
1871 IF (opt_embed%open_shell_embed) THEN
1872 ! If reference system is not spin-polarized then it does not make a contribution to W functional
1873 IF (nspins .EQ. 2) THEN
1874 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) &
1875 - pw_integral_ab(spin_embed_pot, rho_r_ref(1)) &
1876 + pw_integral_ab(spin_embed_pot, rho_r_ref(2))
1877 END IF
1878 END IF
1879 ! Finally, add the regularization term
1880 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + opt_embed%reg_term
1881
1882 ! Print information and check convergence
1883 CALL print_emb_opt_info(output_unit, i_iter, opt_embed)
1884 CALL conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit)
1885 IF (opt_embed%converged) EXIT
1886
1887 ! Update the trust radius and control the step
1888 IF ((i_iter .GT. 1) .AND. (.NOT. opt_embed%steep_desc)) CALL step_control(opt_embed)
1889
1890 ! Print density difference
1891 CALL print_rho_diff(diff_rho_r, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1892 IF (opt_embed%open_shell_embed) THEN ! Spin part
1893 CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1894 END IF
1895
1896 ! Calculate potential gradient if the step has been accepted. Otherwise, we reuse the previous one
1897
1898 IF (opt_embed%accept_step .AND. (.NOT. opt_embed%grid_opt)) &
1899 CALL calculate_embed_pot_grad(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1900 diff_rho_r, diff_rho_spin, opt_embed)
1901 ! Take the embedding step
1902 CALL opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, &
1903 force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
1904
1905 END DO ! i_iter
1906
1907 ! Print final embedding potential for restart
1908 CALL print_embed_restart(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1909 opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
1910 spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .true.)
1911 CALL print_pot_simple_grid(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1912 embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .true., &
1913 force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
1914
1915 ! Print final density difference
1916 !CALL print_rho_diff(diff_rho_r, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .TRUE.)
1917 IF (opt_embed%open_shell_embed) THEN ! Spin part
1918 CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .true.)
1919 END IF
1920
1921 ! Give away plane waves pools
1922 CALL diff_rho_r%release()
1923 IF (opt_embed%open_shell_embed) THEN
1924 CALL diff_rho_spin%release()
1925 END IF
1926
1927 CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
1928 "PRINT%PROGRAM_RUN_INFO")
1929
1930 ! If converged send the embedding potential to the higher-level calculation.
1931 IF (opt_embed%converged) THEN
1932 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, dft_control=dft_control, &
1933 pw_env=pw_env)
1934 nspins_subsys = dft_control%nspins
1935 dft_control%apply_embed_pot = .true.
1936 ! The embedded subsystem corresponds to subsystem #2, where spin change is possible
1937 CALL make_subsys_embed_pot(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, &
1938 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
1939 opt_embed%open_shell_embed, opt_embed%change_spin)
1940
1941 IF (opt_embed%Coulomb_guess) THEN
1942 CALL pw_axpy(opt_embed%pot_diff, embed_pot_subsys, -1.0_dp, allow_noncompatible_grids=.true.)
1943 END IF
1944
1945 CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, embed_pot=embed_pot_subsys)
1946
1947 IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys .EQ. 2)) THEN
1948 CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, &
1949 spin_embed_pot=spin_embed_pot_subsys)
1950 END IF
1951
1952 ! Substitute the correct energy in energies: only on rank 0
1953 IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source()) THEN
1954 energies(cluster_subsys_num) = cluster_energy
1955 END IF
1956 END IF
1957
1958 ! Deallocate and release opt_embed content
1959 CALL release_opt_embed(opt_embed)
1960
1961 ! Deallocate embedding potential
1962 CALL embed_pot%release()
1963 DEALLOCATE (embed_pot)
1964 IF (opt_embed%open_shell_embed) THEN
1965 CALL spin_embed_pot%release()
1966 DEALLOCATE (spin_embed_pot)
1967 END IF
1968
1969 converged_embed = opt_embed%converged
1970
1971 CALL timestop(handle)
1972
1973 END SUBROUTINE dfet_embedding
1974
1975! **************************************************************************************************
1976!> \brief Main driver for the DMFET embedding
1977!> \param force_env ...
1978!> \param ref_subsys_number ...
1979!> \param energies ...
1980!> \param converged_embed ...
1981!> \author Vladimir Rybkin
1982! **************************************************************************************************
1983 SUBROUTINE dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1984 TYPE(force_env_type), POINTER :: force_env
1985 INTEGER :: ref_subsys_number
1986 REAL(kind=dp), DIMENSION(:), POINTER :: energies
1987 LOGICAL :: converged_embed
1988
1989 CHARACTER(LEN=*), PARAMETER :: routinen = 'dmfet_embedding'
1990
1991 INTEGER :: cluster_subsys_num, handle, &
1992 i_force_eval, i_iter, output_unit
1993 LOGICAL :: subsys_open_shell
1994 REAL(kind=dp) :: cluster_energy
1995 TYPE(cp_logger_type), POINTER :: logger
1996 TYPE(dft_control_type), POINTER :: dft_control
1997 TYPE(mp_para_env_type), POINTER :: para_env
1998 TYPE(opt_dmfet_pot_type) :: opt_dmfet
1999 TYPE(qs_energy_type), POINTER :: energy
2000 TYPE(section_vals_type), POINTER :: dft_section, input, opt_dmfet_section
2001
2002 CALL timeset(routinen, handle)
2003
2004 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2005 para_env=para_env)
2006
2007 ! Reveal input file
2008 NULLIFY (logger)
2009 logger => cp_get_default_logger()
2010 output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
2011 extension=".Log")
2012
2013 NULLIFY (dft_section, input, opt_dmfet_section)
2014 NULLIFY (energy)
2015
2016 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2017 energy=energy, input=input)
2018
2019 dft_section => section_vals_get_subs_vals(input, "DFT")
2020 opt_dmfet_section => section_vals_get_subs_vals(input, &
2021 "DFT%QS%OPT_DMFET")
2022
2023 ! We need to understand how to treat spins states
2024 CALL understand_spin_states(force_env, ref_subsys_number, opt_dmfet%change_spin, opt_dmfet%open_shell_embed, &
2025 opt_dmfet%all_nspins)
2026
2027 ! Prepare for the potential optimization
2028 CALL prepare_dmfet_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2029 opt_dmfet, opt_dmfet_section)
2030
2031 ! Get the reference density matrix/matrices
2032 subsys_open_shell = subsys_spin(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2033 CALL build_full_dm(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2034 opt_dmfet%dm_total, subsys_open_shell, opt_dmfet%open_shell_embed, opt_dmfet%dm_total_beta)
2035
2036 ! Check the preliminary DM difference
2037 CALL cp_fm_copy_general(opt_dmfet%dm_total, opt_dmfet%dm_diff, para_env)
2038 IF (opt_dmfet%open_shell_embed) CALL cp_fm_copy_general(opt_dmfet%dm_total_beta, &
2039 opt_dmfet%dm_diff_beta, para_env)
2040
2041 DO i_force_eval = 1, ref_subsys_number - 1
2042
2043 ! Get the subsystem density matrix/matrices
2044 subsys_open_shell = subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2045
2046 CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2047 opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2048 opt_dmfet%dm_subsys_beta)
2049
2050 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2051
2052 IF (opt_dmfet%open_shell_embed) CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, &
2053 1.0_dp, opt_dmfet%dm_subsys_beta)
2054
2055 END DO
2056
2057 ! Main loop of iterative matrix potential optimization
2058 DO i_iter = 1, opt_dmfet%n_iter
2059
2060 opt_dmfet%i_iter = i_iter
2061
2062 ! Set the dm difference as the reference one
2063 CALL cp_fm_copy_general(opt_dmfet%dm_total, opt_dmfet%dm_diff, para_env)
2064
2065 IF (opt_dmfet%open_shell_embed) CALL cp_fm_copy_general(opt_dmfet%dm_total_beta, &
2066 opt_dmfet%dm_diff_beta, para_env)
2067
2068 ! Loop over force evaluations
2069 DO i_force_eval = 1, ref_subsys_number - 1
2070
2071 ! Switch on external potential in the subsystems
2072 NULLIFY (dft_control)
2073 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
2074 dft_control%apply_dmfet_pot = .true.
2075
2076 ! Calculate the new density
2077 CALL force_env_calc_energy_force(force_env=force_env%sub_force_env(i_force_eval)%force_env, &
2078 calc_force=.false., &
2079 skip_external_control=.true.)
2080
2081 ! Extract subsystem density matrix and energy
2082 NULLIFY (energy)
2083
2084 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, energy=energy)
2085 opt_dmfet%w_func(i_iter) = opt_dmfet%w_func(i_iter) + energy%total
2086
2087 ! Find out which subsystem is the cluster
2088 IF (dft_control%qs_control%cluster_embed_subsys) THEN
2089 cluster_subsys_num = i_force_eval
2090 cluster_energy = energy%total
2091 END IF
2092
2093 ! Add subsystem density matrices
2094 subsys_open_shell = subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2095
2096 CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2097 opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2098 opt_dmfet%dm_subsys_beta)
2099
2100 IF (opt_dmfet%open_shell_embed) THEN ! Open-shell embedding
2101 ! We may need to change spin ONLY FOR THE SECOND SUBSYSTEM: that's the internal convention
2102 IF ((i_force_eval .EQ. 2) .AND. (opt_dmfet%change_spin)) THEN
2103 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, 1.0_dp, opt_dmfet%dm_subsys)
2104 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys_beta)
2105 ELSE
2106 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2107 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, 1.0_dp, opt_dmfet%dm_subsys_beta)
2108 END IF
2109 ELSE ! Closed-shell embedding
2110 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff, 1.0_dp, opt_dmfet%dm_subsys)
2111 END IF
2112
2113 END DO ! i_force_eval
2114
2115 CALL check_dmfet(opt_dmfet, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2116
2117 END DO ! i_iter
2118
2119 ! Substitute the correct energy in energies: only on rank 0
2120 IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source()) THEN
2121 energies(cluster_subsys_num) = cluster_energy
2122 END IF
2123
2124 CALL release_dmfet_opt(opt_dmfet)
2125
2126 converged_embed = .false.
2127
2128 END SUBROUTINE dmfet_embedding
2129
2130END MODULE force_env_methods
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:117
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
cp_result_methods::get_results
Definition cp_result_methods.F:43
kahan_sum::accurate_sum
Definition kahan_sum.F:41
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_copy
Definition pw_methods.F:145
pw_methods::pw_integral_ab
Definition pw_methods.F:221
pw_methods::pw_zero
Definition pw_methods.F:82
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
atprop_types::atprop_init
subroutine, public atprop_init(atprop_env, natom)
...
Definition atprop_types.F:66
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:28
bibliography::huang2011
integer, save, public huang2011
Definition bibliography.F:43
bibliography::heaton_burgess2007
integer, save, public heaton_burgess2007
Definition bibliography.F:43
cell_methods
Handles all functions related to the CELL.
Definition cell_methods.F:15
cell_methods::init_cell
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
Definition cell_methods.F:117
cell_methods::cell_create
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition cell_methods.F:85
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::scaled_to_real
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:516
cell_types::cell_sym_triclinic
integer, parameter, public cell_sym_triclinic
Definition cell_types.F:29
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:486
cell_types::cell_release
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:559
cell_types::cell_clone
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition cell_types.F:107
constraint_fxd
Definition constraint_fxd.F:12
constraint_fxd::fix_atom_control
subroutine, public fix_atom_control(force_env, w)
allows for fix atom constraints
Definition constraint_fxd.F:65
constraint_vsite
Routines to handle the virtual site constraint/restraint.
Definition constraint_vsite.F:14
constraint_vsite::vsite_force_control
subroutine, public vsite_force_control(force_env)
control force distribution for virtual sites
Definition constraint_vsite.F:53
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:167
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition cp_fm_types.F:1538
cp_iter_types
Collection of routines to handle the iteration info.
Definition cp_iter_types.F:11
cp_iter_types::cp_iteration_info_copy_iter
subroutine, public cp_iteration_info_copy_iter(iteration_info_in, iteration_info_out)
Copies iterations info of an iteration info into another iteration info.
Definition cp_iter_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_rm_default_logger
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
Definition cp_log_handling.F:213
cp_log_handling::cp_add_default_logger
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
Definition cp_log_handling.F:188
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition cp_result_methods.F:17
cp_result_methods::cp_results_mp_bcast
subroutine, public cp_results_mp_bcast(results, source, para_env)
broadcast results type
Definition cp_result_methods.F:434
cp_result_methods::cp_results_erase
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
Definition cp_result_methods.F:362
cp_result_methods::test_for_result
logical function, public test_for_result(results, description)
test for a certain result in the result_list
Definition cp_result_methods.F:127
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition cp_result_types.F:18
cp_result_types::cp_result_copy
subroutine, public cp_result_copy(results_in, results_out)
Copies the cp_result type.
Definition cp_result_types.F:286
cp_result_types::cp_result_release
subroutine, public cp_result_release(results)
Releases cp_result type.
Definition cp_result_types.F:120
cp_result_types::cp_result_create
subroutine, public cp_result_create(results)
Allocates and intitializes the cp_result.
Definition cp_result_types.F:96
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_set
subroutine, public cp_subsys_set(subsys, atomic_kinds, particles, local_particles, molecules, molecule_kinds, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, results, cell)
sets various propreties of the subsys
Definition cp_subsys_types.F:205
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:345
eip_environment_types
The environment for the empirical interatomic potential methods.
Definition eip_environment_types.F:14
eip_silicon
Empirical interatomic potentials for Silicon.
Definition eip_silicon.F:17
eip_silicon::eip_lenosky
subroutine, public eip_lenosky(eip_env)
Interface routine of Goedecker's Lenosky force field to CP2K.
Definition eip_silicon.F:241
eip_silicon::eip_bazant
subroutine, public eip_bazant(eip_env)
Interface routine of Goedecker's Bazant EDIP to CP2K.
Definition eip_silicon.F:75
embed_types
Definition embed_types.F:8
external_potential_methods
Methods to include the effect of an external potential during an MD or energy calculation.
Definition external_potential_methods.F:13
external_potential_methods::add_external_potential
subroutine, public add_external_potential(force_env)
...
Definition external_potential_methods.F:55
fist_environment_types
Definition fist_environment_types.F:15
fist_force
Definition fist_force.F:16
fist_force::fist_calc_energy_force
subroutine, public fist_calc_energy_force(fist_env, debug)
Calculates the total potential energy, total force, and the total pressure tensor from the potentials...
Definition fist_force.F:112
force_env_methods
Interface for the force calculations.
Definition force_env_methods.F:17
force_env_methods::force_env_calc_energy_force
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Definition force_env_methods.F:193
force_env_methods::force_env_calc_num_pressure
subroutine, public force_env_calc_num_pressure(force_env, dx)
Evaluates the stress tensor and pressure numerically.
Definition force_env_methods.F:551
force_env_methods::force_env_create
subroutine, public force_env_create(force_env, root_section, para_env, globenv, fist_env, qs_env, meta_env, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, force_env_section, mixed_env, embed_env, nnp_env)
creates and initializes a force environment
Definition force_env_methods.F:759
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get_natom
integer function, public force_env_get_natom(force_env)
returns the number of atoms
Definition force_env_types.F:511
force_env_types::use_qmmm
integer, parameter, public use_qmmm
Definition force_env_types.F:79
force_env_types::use_prog_name
character(len=10), dimension(501:509), parameter, public use_prog_name
Definition force_env_types.F:89
force_env_types::use_mixed_force
integer, parameter, public use_mixed_force
Definition force_env_types.F:79
force_env_types::use_eip_force
integer, parameter, public use_eip_force
Definition force_env_types.F:79
force_env_types::use_embed
integer, parameter, public use_embed
Definition force_env_types.F:79
force_env_types::use_qmmmx
integer, parameter, public use_qmmmx
Definition force_env_types.F:79
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
Definition force_env_types.F:329
force_env_types::force_env_set
subroutine, public force_env_set(force_env, meta_env, fp_env, force_env_section, method_name_id, additional_potential)
changes some attributes of the force_env
Definition force_env_types.F:640
force_env_types::use_qs_force
integer, parameter, public use_qs_force
Definition force_env_types.F:79
force_env_types::use_pwdft_force
integer, parameter, public use_pwdft_force
Definition force_env_types.F:79
force_env_types::use_nnp_force
integer, parameter, public use_nnp_force
Definition force_env_types.F:79
force_env_types::use_fist_force
integer, parameter, public use_fist_force
Definition force_env_types.F:79
force_env_utils
Util force_env module.
Definition force_env_utils.F:12
force_env_utils::write_forces
subroutine, public write_forces(particles, iw, label, ndigits, total_force, grand_total_force, zero_force_core_shell_atom)
Write forces.
Definition force_env_utils.F:415
force_env_utils::write_atener
subroutine, public write_atener(iounit, particles, atener, label)
Write the atomic coordinates, types, and energies.
Definition force_env_utils.F:491
force_env_utils::rescale_forces
subroutine, public rescale_forces(force_env)
Rescale forces if requested.
Definition force_env_utils.F:369
force_fields_util
Definition force_fields_util.F:19
force_fields_util::get_generic_info
subroutine, public get_generic_info(gen_section, func_name, xfunction, parameters, values, var_values, size_variables, i_rep_sec, input_variables)
Reads from the input structure all information for generic functions.
Definition force_fields_util.F:1236
fp_methods
methods used in the flexible partitioning scheme
Definition fp_methods.F:14
fp_methods::fp_eval
subroutine, public fp_eval(fp_env, subsys, cell)
computest the forces and the energy due to the flexible potential & bias, and writes the weights file
Definition fp_methods.F:54
fparser
This public domain function parser module is intended for applications where a set of mathematical ex...
Definition fparser.F:17
fparser::parsef
subroutine, public parsef(i, funcstr, var)
Parse ith function string FuncStr and compile it into bytecode.
Definition fparser.F:148
fparser::evalf
real(rn) function, public evalf(i, val)
...
Definition fparser.F:180
fparser::evalerrtype
integer, public evalerrtype
Definition fparser.F:32
fparser::evalfd
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
Definition fparser.F:976
fparser::finalizef
subroutine, public finalizef()
...
Definition fparser.F:101
fparser::initf
subroutine, public initf(n)
...
Definition fparser.F:130
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition global_types.F:21
global_types::globenv_retain
subroutine, public globenv_retain(globenv)
Retains the global environment globenv.
Definition global_types.F:119
grrm_utils
GRRM interface.
Definition grrm_utils.F:12
grrm_utils::write_grrm
subroutine, public write_grrm(iounit, force_env, particles, energy, dipole, hessian, dipder, polar, fixed_atoms)
Write GRRM interface file.
Definition grrm_utils.F:50
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::use_bazant_eip
integer, parameter, public use_bazant_eip
Definition input_constants.F:710
input_constants::mixed_cdft_serial
integer, parameter, public mixed_cdft_serial
Definition input_constants.F:772
input_constants::mix_cdft
integer, parameter, public mix_cdft
Definition input_constants.F:781
input_constants::mix_linear_combination
integer, parameter, public mix_linear_combination
Definition input_constants.F:781
input_constants::dmfet
integer, parameter, public dmfet
Definition input_constants.F:1246
input_constants::use_lenosky_eip
integer, parameter, public use_lenosky_eip
Definition input_constants.F:710
input_constants::mix_coupled
integer, parameter, public mix_coupled
Definition input_constants.F:781
input_constants::debug_run
integer, parameter, public debug_run
Definition input_constants.F:121
input_constants::mix_restrained
integer, parameter, public mix_restrained
Definition input_constants.F:781
input_constants::mix_generic
integer, parameter, public mix_generic
Definition input_constants.F:781
input_constants::dfet
integer, parameter, public dfet
Definition input_constants.F:1246
input_constants::mix_minimum
integer, parameter, public mix_minimum
Definition input_constants.F:781
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_retain
subroutine, public section_vals_retain(section_vals)
retains the given section values (see doc/ReferenceCounting.html)
Definition input_section_types.F:631
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_memory
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition machine.F:332
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::abnormal_value
logical function, public abnormal_value(a)
determines if a value is not normal (e.g. for Inf and Nan) based on IO to work also under optimizatio...
Definition mathlib.F:150
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
metadynamics_types
defines types for metadynamics calculation
Definition metadynamics_types.F:13
mixed_cdft_methods
Methods for mixed CDFT calculations.
Definition mixed_cdft_methods.F:14
mixed_cdft_methods::mixed_cdft_calculate_coupling
subroutine, public mixed_cdft_calculate_coupling(force_env)
Driver routine to calculate the electronic coupling(s) between CDFT states.
Definition mixed_cdft_methods.F:1109
mixed_cdft_methods::mixed_cdft_build_weight
subroutine, public mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
Driver routine to handle the build of CDFT weight/gradient in parallel and serial modes.
Definition mixed_cdft_methods.F:399
mixed_cdft_methods::mixed_cdft_init
subroutine, public mixed_cdft_init(force_env, calculate_forces)
Initialize a mixed CDFT calculation.
Definition mixed_cdft_methods.F:139
mixed_energy_types
Definition mixed_energy_types.F:12
mixed_environment_types
Definition mixed_environment_types.F:11
mixed_environment_types::get_mixed_env
subroutine, public get_mixed_env(mixed_env, atomic_kind_set, particle_set, local_particles, local_molecules, molecule_kind_set, molecule_set, cell, cell_ref, mixed_energy, para_env, sub_para_env, subsys, input, results, cdft_control)
Get the MIXED environment.
Definition mixed_environment_types.F:130
mixed_environment_utils
Util mixed_environment.
Definition mixed_environment_utils.F:12
mixed_environment_utils::get_subsys_map_index
subroutine, public get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, map_index, force_eval_embed)
performs mapping of the subsystems of different force_eval
Definition mixed_environment_utils.F:141
mixed_environment_utils::mixed_map_forces
subroutine, public mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, factor, iforce_eval, nforce_eval, map_index, mapping_section, overwrite)
Maps forces between the different force_eval sections/environments.
Definition mixed_environment_utils.F:63
molecule_kind_list_types
represent a simple array based list of the given type
Definition molecule_kind_list_types.F:15
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:510
nnp_environment_types
Data types for neural network potentials.
Definition nnp_environment_types.F:13
nnp_force
Methods dealing with Neural Network potentials.
Definition nnp_force.F:13
nnp_force::nnp_calc_energy_force
subroutine, public nnp_calc_energy_force(nnp, calc_forces)
Calculate the energy and force for a given configuration with the NNP.
Definition nnp_force.F:62
optimize_dmfet_potential
Definition optimize_dmfet_potential.F:8
optimize_dmfet_potential::check_dmfet
subroutine, public check_dmfet(opt_dmfet, qs_env)
...
Definition optimize_dmfet_potential.F:306
optimize_dmfet_potential::release_dmfet_opt
subroutine, public release_dmfet_opt(opt_dmfet)
...
Definition optimize_dmfet_potential.F:188
optimize_dmfet_potential::build_full_dm
subroutine, public build_full_dm(qs_env, dm, open_shell, open_shell_embed, dm_beta)
Builds density matrices from MO coefficients in full matrix format.
Definition optimize_dmfet_potential.F:250
optimize_dmfet_potential::prepare_dmfet_opt
subroutine, public prepare_dmfet_opt(qs_env, opt_dmfet, opt_dmfet_section)
...
Definition optimize_dmfet_potential.F:112
optimize_dmfet_potential::subsys_spin
logical function, public subsys_spin(qs_env)
...
Definition optimize_dmfet_potential.F:93
optimize_embedding_potential
Definition optimize_embedding_potential.F:8
optimize_embedding_potential::print_embed_restart
subroutine, public print_embed_restart(qs_env, dimen_aux, embed_pot_coef, embed_pot, i_iter, embed_pot_spin, open_shell_embed, grid_opt, final_one)
Print embedding potential as a cube and as a binary (for restarting)
Definition optimize_embedding_potential.F:3119
optimize_embedding_potential::make_subsys_embed_pot
subroutine, public make_subsys_embed_pot(qs_env, embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, open_shell_embed, change_spin_sign)
Creates a subsystem embedding potential.
Definition optimize_embedding_potential.F:1024
optimize_embedding_potential::print_rho_spin_diff
subroutine, public print_rho_spin_diff(spin_diff_rho_r, i_iter, qs_env, final_one)
Prints a cube for the (spin_rho_A + spin_rho_B - spin_rho_ref) to be minimized in embedding.
Definition optimize_embedding_potential.F:3064
optimize_embedding_potential::get_max_subsys_diff
subroutine, public get_max_subsys_diff(opt_embed, force_env, subsys_num)
...
Definition optimize_embedding_potential.F:3446
optimize_embedding_potential::print_emb_opt_info
subroutine, public print_emb_opt_info(output_unit, step_num, opt_embed)
...
Definition optimize_embedding_potential.F:3383
optimize_embedding_potential::init_embed_pot
subroutine, public init_embed_pot(qs_env, embed_pot, add_const_pot, fermi_amaldi, const_pot, open_shell_embed, spin_embed_pot, pot_diff, coulomb_guess, grid_opt)
...
Definition optimize_embedding_potential.F:206
optimize_embedding_potential::understand_spin_states
subroutine, public understand_spin_states(force_env, ref_subsys_number, change_spin, open_shell_embed, all_nspins)
Find out whether we need to swap alpha- and beta- spind densities in the second subsystem.
Definition optimize_embedding_potential.F:128
optimize_embedding_potential::read_embed_pot
subroutine, public read_embed_pot(qs_env, embed_pot, spin_embed_pot, section, opt_embed)
...
Definition optimize_embedding_potential.F:663
optimize_embedding_potential::get_prev_density
subroutine, public get_prev_density(opt_embed, force_env, subsys_num)
...
Definition optimize_embedding_potential.F:3416
optimize_embedding_potential::coulomb_guess
subroutine, public coulomb_guess(v_rspace, rhs, mapping_section, qs_env, nforce_eval, iforce_eval, eta)
Calculates subsystem Coulomb potential from the RESP charges of the total system.
Definition optimize_embedding_potential.F:933
optimize_embedding_potential::print_rho_diff
subroutine, public print_rho_diff(diff_rho_r, i_iter, qs_env, final_one)
Prints a cube for the (rho_A + rho_B - rho_ref) to be minimized in embedding.
Definition optimize_embedding_potential.F:3013
optimize_embedding_potential::conv_check_embed
subroutine, public conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit)
...
Definition optimize_embedding_potential.F:3479
optimize_embedding_potential::step_control
subroutine, public step_control(opt_embed)
Controls the step, changes the trust radius if needed in maximization of the V_emb.
Definition optimize_embedding_potential.F:2170
optimize_embedding_potential::opt_embed_step
subroutine, public opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, qs_env)
Takes maximization step in embedding potential optimization.
Definition optimize_embedding_potential.F:1467
optimize_embedding_potential::calculate_embed_pot_grad
subroutine, public calculate_embed_pot_grad(qs_env, diff_rho_r, diff_rho_spin, opt_embed)
Calculates the derivative of the embedding potential wrt to the expansion coefficients.
Definition optimize_embedding_potential.F:1072
optimize_embedding_potential::print_pot_simple_grid
subroutine, public print_pot_simple_grid(qs_env, embed_pot, embed_pot_spin, i_iter, open_shell_embed, final_one, qs_env_cluster)
Prints a volumetric file: X Y Z value for interfacing with external programs.
Definition optimize_embedding_potential.F:3217
optimize_embedding_potential::prepare_embed_opt
subroutine, public prepare_embed_opt(qs_env, opt_embed, opt_embed_section)
Creates and allocates objects for optimization of embedding potential.
Definition optimize_embedding_potential.F:291
optimize_embedding_potential::release_opt_embed
subroutine, public release_opt_embed(opt_embed)
Deallocate stuff for optimizing embedding potential.
Definition optimize_embedding_potential.F:856
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::debye
real(kind=dp), parameter, public debye
Definition physcon.F:201
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
pwdft_environment_types
The type definitions for the PWDFT environment.
Definition pwdft_environment_types.F:14
pwdft_environment
Methods and functions on the PWDFT environment.
Definition pwdft_environment.F:14
pwdft_environment::pwdft_calc_energy_force
subroutine, public pwdft_calc_energy_force(pwdft_env, calculate_forces, calculate_stress)
Calculate energy and forces within the PWDFT/SIRIUS code.
Definition pwdft_environment.F:210
qmmm_force
Calculates QM/MM energy and forces.
Definition qmmm_force.F:14
qmmm_force::qmmm_calc_energy_force
subroutine, public qmmm_calc_energy_force(qmmm_env, calc_force, consistent_energies, linres)
calculates the qm/mm energy and forces
Definition qmmm_force.F:76
qmmm_types
Basic container type for QM/MM.
Definition qmmm_types.F:12
qmmm_util
Definition qmmm_util.F:13
qmmm_util::apply_qmmm_translate
subroutine, public apply_qmmm_translate(qmmm_env)
Apply translation to the full system in order to center the QM system into the QM box.
Definition qmmm_util.F:373
qmmmx_force
Calculates QM/MM energy and forces with Force-Mixing.
Definition qmmmx_force.F:14
qmmmx_force::qmmmx_calc_energy_force
subroutine, public qmmmx_calc_energy_force(qmmmx_env, calc_force, consistent_energies, linres, require_consistent_energy_force)
calculates the qm/mm energy and forces
Definition qmmmx_force.F:66
qmmmx_types
Basic container type for QM/MM with force mixing.
Definition qmmmx_types.F:12
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1035
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_force
Quickstep force driver routine.
Definition qs_force.F:12
qs_force::qs_calc_energy_force
subroutine, public qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres)
...
Definition qs_force.F:107
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
restraint
Handles all possible kinds of restraints in CP2K.
Definition restraint.F:14
restraint::restraint_control
subroutine, public restraint_control(force_env)
Computes restraints.
Definition restraint.F:67
scine_utils
SCINE interface.
Definition scine_utils.F:12
scine_utils::write_scine
subroutine, public write_scine(iounit, force_env, particles, energy, hessian)
Write SCINE interface file.
Definition scine_utils.F:46
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::compress
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Definition string_utilities.F:3190
virial_methods
Definition virial_methods.F:12
virial_methods::write_stress_tensor_components
subroutine, public write_stress_tensor_components(virial, iw, cell)
...
Definition virial_methods.F:163
virial_methods::write_stress_tensor
subroutine, public write_stress_tensor(pv_virial, iw, cell, numerical)
Print stress tensor to output file.
Definition virial_methods.F:254
virial_types
Definition virial_types.F:13
virial_types::zero_virial
subroutine, public zero_virial(virial, reset)
...
Definition virial_types.F:123
virial_types::symmetrize_virial
subroutine, public symmetrize_virial(virial)
Symmetrize the virial components.
Definition virial_types.F:67
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_result_types::cp_result_p_type
Definition cp_result_types.F:81
cp_result_types::cp_result_type
contains arbitrary information which need to be stored
Definition cp_result_types.F:73
cp_subsys_types::cp_subsys_p_type
represent a pointer to a subsys, to be able to create arrays of pointers
Definition cp_subsys_types.F:122
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
eip_environment_types::eip_environment_type
The empirical interatomic potential environment.
Definition eip_environment_types.F:83
embed_types::embed_env_type
Embedding environment type.
Definition embed_types.F:127
embed_types::opt_dmfet_pot_type
Type containing main data for matrix embedding potential optimization.
Definition embed_types.F:105
embed_types::opt_embed_pot_type
Type containing main data for embedding potential optimization.
Definition embed_types.F:67
fist_environment_types::fist_environment_type
Definition fist_environment_types.F:71
force_env_types::force_env_p_type
allows for the creation of an array of force_env
Definition force_env_types.F:167
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:139
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
metadynamics_types::meta_env_type
defines meta_env type
Definition metadynamics_types.F:94
mixed_energy_types::mixed_energy_type
Definition mixed_energy_types.F:22
mixed_energy_types::mixed_force_type
Definition mixed_energy_types.F:28
mixed_environment_types::mixed_environment_type
Definition mixed_environment_types.F:64
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
nnp_environment_types::nnp_type
Main data type collecting all relevant data for neural network potentials.
Definition nnp_environment_types.F:83
particle_list_types::particle_list_p_type
represents a pointer to a list
Definition particle_list_types.F:59
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
pwdft_environment_types::pwdft_environment_type
The PWDFT environment type.
Definition pwdft_environment_types.F:67
qmmm_types::qmmm_env_type
Definition qmmm_types.F:34
qmmmx_types::qmmmx_env_type
Definition qmmmx_types.F:27
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
virial_types::virial_p_type
Definition virial_types.F:52
virial_types::virial_type
Definition virial_types.F:26