LCOV - code coverage report
Current view: top level - src - force_env_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 833 909 91.6 %
Date: 2024-03-28 07:31:50 Functions: 9 10 90.0 %

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

Generated by: LCOV version 1.15