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