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