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