165#include "./base/base_uses.f90"
171 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'force_env_methods'
177 INTEGER,
SAVE,
PRIVATE :: last_force_env_id = 0
197 consistent_energies, skip_external_control, eval_energy_forces, &
198 require_consistent_energy_force, linres, calc_stress_tensor)
201 LOGICAL,
INTENT(IN),
OPTIONAL :: calc_force, consistent_energies, skip_external_control, &
202 eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor
204 REAL(kind=
dp),
PARAMETER :: ateps = 1.0e-6_dp
206 CHARACTER(LEN=default_string_length) :: unit_string
207 INTEGER :: ikind, nat, ndigits, nfixed_atoms, &
208 nfixed_atoms_total, nkind, &
209 output_unit, print_forces, print_grrm, &
211 LOGICAL :: calculate_forces, calculate_stress_tensor, do_apt_fd, energy_consistency, &
212 eval_ef, linres_run, my_skip, print_components
213 REAL(kind=
dp) :: checksum, e_entropy, e_gap, e_pot, &
215 REAL(kind=
dp),
DIMENSION(3) :: grand_total_force, total_force
228 NULLIFY (logger, virial, subsys, atprop_env, cell)
232 calculate_forces = .true.
233 energy_consistency = .false.
239 IF (
PRESENT(eval_energy_forces)) eval_ef = eval_energy_forces
240 IF (
PRESENT(skip_external_control)) my_skip = skip_external_control
241 IF (
PRESENT(calc_force)) calculate_forces = calc_force
242 IF (
PRESENT(calc_stress_tensor))
THEN
243 calculate_stress_tensor = calc_stress_tensor
245 calculate_stress_tensor = calculate_forces
247 IF (
PRESENT(consistent_energies)) energy_consistency = consistent_energies
248 IF (
PRESENT(linres)) linres_run = linres
250 cpassert(
ASSOCIATED(force_env))
251 cpassert(force_env%ref_count > 0)
254 CALL cp_subsys_get(subsys, virial=virial, atprop=atprop_env, cell=cell)
255 IF (virial%pv_availability)
CALL zero_virial(virial, reset=.false.)
260 SELECT CASE (force_env%in_use)
266 IF (virial%pv_availability .AND. calculate_stress_tensor)
THEN
271 e_gap = force_env%pwdft_env%energy%band_gap
272 e_entropy = force_env%pwdft_env%energy%entropy
281 calculate_forces, energy_consistency, linres=linres_run)
284 calculate_forces, energy_consistency, linres=linres_run, &
285 require_consistent_energy_force=require_consistent_energy_force)
287 CALL mixed_energy_forces(force_env, calculate_forces)
292 CALL embed_energy(force_env)
300 IF (virial%pv_availability)
THEN
301 IF (virial%pv_numer .AND. calculate_stress_tensor)
THEN
305 IF (calculate_forces)
THEN
309 IF (calculate_stress_tensor)
THEN
310 CALL cp_warn(__location__,
"The calculation of the stress tensor "// &
311 "requires the calculation of the forces")
320 CALL section_vals_val_get(force_env%qs_env%input,
"PROPERTIES%LINRES%DCDR%APT_FD", l_val=do_apt_fd)
323 subsection_name=
"PROPERTIES%LINRES%DCDR%PRINT%APT")
334 IF (.NOT. my_skip)
THEN
336 IF (
ASSOCIATED(force_env%fp_env))
THEN
337 IF (force_env%fp_env%use_fp)
THEN
338 CALL fp_eval(force_env%fp_env, subsys, cell)
356 output_unit =
cp_print_key_unit_nr(logger, force_env%force_env_section,
"PRINT%PROGRAM_RUN_INFO", &
358 IF (output_unit > 0)
THEN
362 WRITE (unit=output_unit, fmt=
"(/,T2,A,T55,F26.15)") &
363 "ENERGY| Total FORCE_EVAL ( "//trim(adjustl(
use_prog_name(force_env%in_use)))// &
364 " ) energy ["//trim(adjustl(unit_string))//
"]", e_pot*fconv
365 IF (e_gap > -0.1_dp)
THEN
366 WRITE (unit=output_unit, fmt=
"(/,T2,A,T55,F26.15)") &
367 "ENERGY| Total FORCE_EVAL ( "//trim(adjustl(
use_prog_name(force_env%in_use)))// &
368 " ) gap ["//trim(adjustl(unit_string))//
"]", e_gap*fconv
370 IF (e_entropy > -0.1_dp)
THEN
371 WRITE (unit=output_unit, fmt=
"(/,T2,A,T55,F26.15)") &
372 "ENERGY| Total FORCE_EVAL ( "//trim(adjustl(
use_prog_name(force_env%in_use)))// &
373 " ) free energy ["//trim(adjustl(unit_string))//
"]", (e_pot - e_entropy)*fconv
377 "PRINT%PROGRAM_RUN_INFO")
381 cpabort(
"Potential energy is an abnormal value (NaN/Inf).")
386 IF ((print_forces > 0) .AND. calculate_forces)
THEN
389 core_particles=core_particles, &
390 particles=particles, &
391 shell_particles=shell_particles)
397 IF (
ASSOCIATED(core_particles) .OR.
ASSOCIATED(shell_particles))
THEN
398 CALL write_forces(particles, print_forces,
"Atomic", ndigits, unit_string, &
399 total_force, zero_force_core_shell_atom=.true.)
400 grand_total_force(1:3) = total_force(1:3)
401 IF (
ASSOCIATED(core_particles))
THEN
402 CALL write_forces(core_particles, print_forces,
"Core particle", ndigits, &
403 unit_string, total_force, zero_force_core_shell_atom=.false.)
404 grand_total_force(:) = grand_total_force(:) + total_force(:)
406 IF (
ASSOCIATED(shell_particles))
THEN
407 CALL write_forces(shell_particles, print_forces,
"Shell particle", ndigits, &
408 unit_string, total_force, zero_force_core_shell_atom=.false., &
409 grand_total_force=grand_total_force)
412 CALL write_forces(particles, print_forces,
"Atomic", ndigits, unit_string, total_force)
418 IF (virial%pv_availability)
THEN
421 IF (calculate_forces .AND. calculate_stress_tensor)
THEN
422 output_unit =
cp_print_key_unit_nr(logger, force_env%force_env_section,
"PRINT%STRESS_TENSOR", &
423 extension=
".stress_tensor")
424 IF (output_unit > 0)
THEN
426 l_val=print_components)
429 IF (print_components)
THEN
430 IF ((.NOT. virial%pv_numer) .AND. (force_env%in_use ==
use_qs_force))
THEN
434 CALL write_stress_tensor(virial%pv_virial, output_unit, cell, unit_string, virial%pv_numer)
437 "PRINT%STRESS_TENSOR")
442 output_unit =
cp_print_key_unit_nr(logger, force_env%force_env_section,
"PRINT%STRESS_TENSOR", &
443 extension=
".stress_tensor")
444 IF (output_unit > 0)
THEN
445 CALL cp_warn(__location__,
"To print the stress tensor switch on the "// &
446 "virial evaluation with the keyword: STRESS_TENSOR")
449 "PRINT%STRESS_TENSOR")
453 output_unit =
cp_print_key_unit_nr(logger, force_env%force_env_section,
"PRINT%PROGRAM_RUN_INFO", &
455 IF (atprop_env%energy)
THEN
456 CALL force_env%para_env%sum(atprop_env%atener)
458 IF (output_unit > 0)
THEN
461 CALL write_atener(output_unit, particles, atprop_env%atener,
"Mulliken Atomic Energies")
464 checksum = abs(e_pot - sum_energy)
465 WRITE (unit=output_unit, fmt=
"(/,(T2,A,T56,F25.13))") &
466 "Potential energy (Atomic):", sum_energy, &
467 "Potential energy (Total) :", e_pot, &
468 "Difference :", checksum
469 cpassert((checksum < ateps*abs(e_pot)))
472 "PRINT%PROGRAM_RUN_INFO")
477 file_position=
"REWIND", extension=
".rrm")
478 IF (print_grrm > 0)
THEN
481 molecule_kinds=molecule_kinds)
483 nfixed_atoms_total = 0
484 nkind = molecule_kinds%n_els
485 molecule_kind_set => molecule_kinds%els
487 molecule_kind => molecule_kind_set(ikind)
489 nfixed_atoms_total = nfixed_atoms_total + nfixed_atoms
492 CALL write_grrm(print_grrm, force_env, particles%els, e_pot, fixed_atoms=nfixed_atoms_total)
497 file_position=
"REWIND", extension=
".scine")
498 IF (print_scine > 0)
THEN
502 CALL write_scine(print_scine, force_env, particles%els, e_pot)
521 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: dx
523 REAL(kind=
dp),
PARAMETER :: default_dx = 0.001_dp
525 CHARACTER(LEN=default_string_length) :: unit_string
526 INTEGER :: i, ip, iq, j, k, natom, ncore, nshell, &
527 output_unit, symmetry_id
528 REAL(kind=
dp) :: dx_w
529 REAL(kind=
dp),
DIMENSION(2) :: numer_energy
530 REAL(kind=
dp),
DIMENSION(3) :: s
531 REAL(kind=
dp),
DIMENSION(3, 3) :: numer_stress
532 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: ref_pos_atom, ref_pos_core, ref_pos_shell
533 TYPE(
cell_type),
POINTER :: cell, cell_local
542 NULLIFY (core_particles)
544 NULLIFY (shell_particles)
545 NULLIFY (ref_pos_atom)
546 NULLIFY (ref_pos_core)
547 NULLIFY (ref_pos_shell)
551 numer_stress = 0.0_dp
556 IF (
PRESENT(dx)) dx_w = dx
557 CALL force_env_get(force_env, subsys=subsys, globenv=globenv)
559 core_particles=core_particles, &
560 particles=particles, &
561 shell_particles=shell_particles, &
563 output_unit =
cp_print_key_unit_nr(logger, force_env%force_env_section,
"PRINT%STRESS_TENSOR", &
564 extension=
".stress_tensor")
565 IF (output_unit > 0)
THEN
566 WRITE (output_unit,
"(/A,A/)")
" **************************** ", &
567 "NUMERICAL STRESS ********************************"
571 natom = particles%n_els
572 ALLOCATE (ref_pos_atom(natom, 3))
574 ref_pos_atom(i, :) = particles%els(i)%r
576 IF (
ASSOCIATED(core_particles))
THEN
577 ncore = core_particles%n_els
578 ALLOCATE (ref_pos_core(ncore, 3))
580 ref_pos_core(i, :) = core_particles%els(i)%r
583 IF (
ASSOCIATED(shell_particles))
THEN
584 nshell = shell_particles%n_els
585 ALLOCATE (ref_pos_shell(nshell, 3))
587 ref_pos_shell(i, :) = shell_particles%els(i)%r
592 symmetry_id = cell%symmetry_id
600 IF (virial%pv_diagonal .AND. (ip /= iq)) cycle
602 cell%hmat(ip, iq) = cell_local%hmat(ip, iq) - (-1.0_dp)**k*dx_w
619 calc_force=.false., &
620 consistent_energies=.true., &
621 calc_stress_tensor=.false.)
622 CALL force_env_get(force_env, potential_energy=numer_energy(k))
624 cell%hmat(ip, iq) = cell_local%hmat(ip, iq)
627 numer_stress(ip, iq) = 0.5_dp*(numer_energy(1) - numer_energy(2))/dx_w
628 IF (output_unit > 0)
THEN
629 IF (globenv%run_type_id ==
debug_run)
THEN
630 WRITE (unit=output_unit, fmt=
"(/,T2,A,T19,A,F7.4,A,T44,A,F7.4,A,T69,A)") &
631 "DEBUG|",
"E("//achar(119 + ip)//achar(119 + iq)//
" +", dx_w,
")", &
632 "E("//achar(119 + ip)//achar(119 + iq)//
" -", dx_w,
")", &
634 WRITE (unit=output_unit, fmt=
"(T2,A,2(1X,F24.8),1X,F22.8)") &
635 "DEBUG|", numer_energy(1:2), numer_stress(ip, iq)
637 WRITE (unit=output_unit, fmt=
"(/,T7,A,F7.4,A,T27,A,F7.4,A,T49,A)") &
638 "E("//achar(119 + ip)//achar(119 + iq)//
" +", dx_w,
")", &
639 "E("//achar(119 + ip)//achar(119 + iq)//
" -", dx_w,
")", &
641 WRITE (unit=output_unit, fmt=
"(3(1X,F19.8))") &
642 numer_energy(1:2), numer_stress(ip, iq)
649 cell%symmetry_id = symmetry_id
652 particles%els(i)%r = ref_pos_atom(i, :)
655 core_particles%els(i)%r = ref_pos_core(i, :)
658 shell_particles%els(i)%r = ref_pos_shell(i, :)
661 calc_force=.false., &
662 consistent_energies=.true., &
663 calc_stress_tensor=.false.)
666 virial%pv_virial = 0.0_dp
670 virial%pv_virial(i, j) = virial%pv_virial(i, j) - &
671 0.5_dp*(numer_stress(i, k)*cell_local%hmat(j, k) + &
672 numer_stress(j, k)*cell_local%hmat(i, k))
677 IF (output_unit > 0)
THEN
678 IF (globenv%run_type_id ==
debug_run)
THEN
681 CALL write_stress_tensor(virial%pv_virial, output_unit, cell, unit_string, virial%pv_numer)
683 WRITE (output_unit,
"(/,A,/)")
" **************************** "// &
684 "NUMERICAL STRESS END *****************************"
688 "PRINT%STRESS_TENSOR")
691 IF (
ASSOCIATED(ref_pos_atom))
THEN
692 DEALLOCATE (ref_pos_atom)
694 IF (
ASSOCIATED(ref_pos_core))
THEN
695 DEALLOCATE (ref_pos_core)
697 IF (
ASSOCIATED(ref_pos_shell))
THEN
698 DEALLOCATE (ref_pos_shell)
700 IF (
ASSOCIATED(cell_local))
CALL cell_release(cell_local)
729 qs_env, meta_env, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, force_env_section, &
730 mixed_env, embed_env, nnp_env, ipi_env)
740 POINTER :: sub_force_env
748 TYPE(
nnp_type),
OPTIONAL,
POINTER :: nnp_env
752 NULLIFY (force_env%fist_env, force_env%qs_env, &
753 force_env%para_env, force_env%globenv, &
754 force_env%meta_env, force_env%sub_force_env, &
755 force_env%qmmm_env, force_env%qmmmx_env, force_env%fp_env, &
756 force_env%force_env_section, force_env%eip_env, force_env%mixed_env, &
757 force_env%embed_env, force_env%pwdft_env, force_env%nnp_env, &
758 force_env%root_section)
759 last_force_env_id = last_force_env_id + 1
760 force_env%ref_count = 1
762 force_env%additional_potential = 0.0_dp
764 force_env%globenv => globenv
767 force_env%root_section => root_section
770 force_env%para_env => para_env
771 CALL force_env%para_env%retain()
774 force_env%force_env_section => force_env_section
776 IF (
PRESENT(fist_env))
THEN
777 cpassert(
ASSOCIATED(fist_env))
778 cpassert(force_env%in_use == 0)
780 force_env%fist_env => fist_env
782 IF (
PRESENT(eip_env))
THEN
783 cpassert(
ASSOCIATED(eip_env))
784 cpassert(force_env%in_use == 0)
786 force_env%eip_env => eip_env
788 IF (
PRESENT(pwdft_env))
THEN
789 cpassert(
ASSOCIATED(pwdft_env))
790 cpassert(force_env%in_use == 0)
792 force_env%pwdft_env => pwdft_env
794 IF (
PRESENT(qs_env))
THEN
795 cpassert(
ASSOCIATED(qs_env))
796 cpassert(force_env%in_use == 0)
798 force_env%qs_env => qs_env
800 IF (
PRESENT(qmmm_env))
THEN
801 cpassert(
ASSOCIATED(qmmm_env))
802 cpassert(force_env%in_use == 0)
804 force_env%qmmm_env => qmmm_env
806 IF (
PRESENT(qmmmx_env))
THEN
807 cpassert(
ASSOCIATED(qmmmx_env))
808 cpassert(force_env%in_use == 0)
810 force_env%qmmmx_env => qmmmx_env
812 IF (
PRESENT(mixed_env))
THEN
813 cpassert(
ASSOCIATED(mixed_env))
814 cpassert(force_env%in_use == 0)
816 force_env%mixed_env => mixed_env
818 IF (
PRESENT(embed_env))
THEN
819 cpassert(
ASSOCIATED(embed_env))
820 cpassert(force_env%in_use == 0)
822 force_env%embed_env => embed_env
824 IF (
PRESENT(nnp_env))
THEN
825 cpassert(
ASSOCIATED(nnp_env))
826 cpassert(force_env%in_use == 0)
828 force_env%nnp_env => nnp_env
830 IF (
PRESENT(ipi_env))
THEN
831 cpassert(
ASSOCIATED(ipi_env))
832 cpassert(force_env%in_use == 0)
834 force_env%ipi_env => ipi_env
836 cpassert(force_env%in_use /= 0)
838 IF (
PRESENT(sub_force_env))
THEN
839 force_env%sub_force_env => sub_force_env
842 IF (
PRESENT(meta_env))
THEN
843 force_env%meta_env => meta_env
845 NULLIFY (force_env%meta_env)
866 SUBROUTINE mixed_energy_forces(force_env, calculate_forces)
869 LOGICAL,
INTENT(IN) :: calculate_forces
871 CHARACTER(LEN=default_path_length) :: coupling_function
872 CHARACTER(LEN=default_string_length) :: def_error, description, this_error
873 INTEGER :: iforce_eval, iparticle, istate(2), &
874 jparticle, mixing_type, my_group, &
875 natom, nforce_eval, source, unit_nr
876 INTEGER,
DIMENSION(:),
POINTER :: glob_natoms, itmplist, map_index
877 LOGICAL :: dip_exists
878 REAL(kind=
dp) :: coupling_parameter, dedf, der_1, der_2, &
879 dx, energy, err, lambda, lerr, &
880 restraint_strength, restraint_target, &
882 REAL(kind=
dp),
DIMENSION(3) :: dip_mix
883 REAL(kind=
dp),
DIMENSION(:),
POINTER :: energies
895 mapping_section, mixed_section, &
898 TYPE(
virial_type),
POINTER :: loc_virial, virial_mix
901 cpassert(
ASSOCIATED(force_env))
905 force_env_section=force_env_section, &
906 root_section=root_section, &
909 particles=particles_mix, &
912 NULLIFY (map_index, glob_natoms, global_forces, itmplist)
914 nforce_eval =
SIZE(force_env%sub_force_env)
918 ALLOCATE (subsystems(nforce_eval))
919 ALLOCATE (particles(nforce_eval))
921 ALLOCATE (global_forces(nforce_eval))
922 ALLOCATE (energies(nforce_eval))
923 ALLOCATE (glob_natoms(nforce_eval))
924 ALLOCATE (virials(nforce_eval))
925 ALLOCATE (results(nforce_eval))
932 IF (.NOT. force_env%mixed_env%do_mixed_cdft)
THEN
933 DO iforce_eval = 1, nforce_eval
934 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
935 NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
936 ALLOCATE (virials(iforce_eval)%virial)
938 IF (.NOT.
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
940 my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
941 my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
947 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
948 subsys=subsystems(iforce_eval)%subsys)
951 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
955 particles=particles(iforce_eval)%list)
958 natom =
SIZE(particles(iforce_eval)%list%els)
964 DO iparticle = 1, natom
965 jparticle = map_index(iparticle)
966 particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
971 calc_force=calculate_forces, &
972 skip_external_control=.true.)
975 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
976 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
977 potential_energy=energy)
979 virial=loc_virial, results=loc_results)
980 energies(iforce_eval) = energy
981 glob_natoms(iforce_eval) = natom
982 virials(iforce_eval)%virial = loc_virial
986 IF (
ASSOCIATED(map_index))
THEN
987 DEALLOCATE (map_index)
992 CALL mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
993 glob_natoms, virials, results)
996 CALL force_env%para_env%sync()
998 CALL mixed_cdft_post_energy_forces(force_env)
1000 CALL force_env%para_env%sum(energies)
1001 CALL force_env%para_env%sum(glob_natoms)
1003 DO iforce_eval = 1, nforce_eval
1004 ALLOCATE (global_forces(iforce_eval)%forces(3, glob_natoms(iforce_eval)))
1005 global_forces(iforce_eval)%forces = 0.0_dp
1006 IF (
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env))
THEN
1007 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
1009 DO iparticle = 1, glob_natoms(iforce_eval)
1010 global_forces(iforce_eval)%forces(:, iparticle) = &
1011 particles(iforce_eval)%list%els(iparticle)%f
1015 CALL force_env%para_env%sum(global_forces(iforce_eval)%forces)
1017 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_total)
1018 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_kinetic)
1019 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_virial)
1020 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_xc)
1021 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_fock_4c)
1022 CALL force_env%para_env%sum(virials(iforce_eval)%virial%pv_constraint)
1025 IF (
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env))
THEN
1026 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) &
1027 source = force_env%para_env%mepos
1029 CALL force_env%para_env%sum(source)
1033 force_env%mixed_env%energies = energies
1036 mixed_energy=mixed_energy)
1040 DO iparticle = 1,
SIZE(particles_mix%els)
1041 particles_mix%els(iparticle)%f(:) = 0.0_dp
1045 SELECT CASE (mixing_type)
1048 cpassert(nforce_eval == 2)
1050 mixed_energy%pot = lambda*energies(1) + (1.0_dp - lambda)*energies(2)
1052 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1053 lambda, 1, nforce_eval, map_index, mapping_section, .true.)
1054 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1055 (1.0_dp - lambda), 2, nforce_eval, map_index, mapping_section, .false.)
1058 cpassert(nforce_eval == 2)
1059 IF (energies(1) < energies(2))
THEN
1060 mixed_energy%pot = energies(1)
1061 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1062 1.0_dp, 1, nforce_eval, map_index, mapping_section, .true.)
1064 mixed_energy%pot = energies(2)
1065 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1066 1.0_dp, 2, nforce_eval, map_index, mapping_section, .true.)
1070 cpassert(nforce_eval == 2)
1072 r_val=coupling_parameter)
1073 sd = sqrt((energies(1) - energies(2))**2 + 4.0_dp*coupling_parameter**2)
1074 der_1 = (1.0_dp - (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1075 der_2 = (1.0_dp + (1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1) - energies(2)))/2.0_dp
1076 mixed_energy%pot = (energies(1) + energies(2) - sd)/2.0_dp
1078 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1079 der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1080 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1081 der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1084 cpassert(nforce_eval == 2)
1086 r_val=restraint_target)
1088 r_val=restraint_strength)
1089 mixed_energy%pot = energies(1) + restraint_strength*(energies(1) - energies(2) - restraint_target)**2
1090 der_2 = -2.0_dp*restraint_strength*(energies(1) - energies(2) - restraint_target)
1091 der_1 = 1.0_dp - der_2
1093 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1094 der_1, 1, nforce_eval, map_index, mapping_section, .true.)
1095 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1096 der_2, 2, nforce_eval, map_index, mapping_section, .false.)
1100 CALL get_generic_info(gen_section,
"MIXING_FUNCTION", coupling_function, force_env%mixed_env%par, &
1101 force_env%mixed_env%val, energies)
1103 CALL parsef(1, trim(coupling_function), force_env%mixed_env%par)
1105 mixed_energy%pot =
evalf(1, force_env%mixed_env%val)
1109 DO iforce_eval = 1, nforce_eval
1112 dedf =
evalfd(1, iforce_eval, force_env%mixed_env%val, dx, err)
1113 IF (abs(err) > lerr)
THEN
1114 WRITE (this_error,
"(A,G12.6,A)")
"(", err,
")"
1115 WRITE (def_error,
"(A,G12.6,A)")
"(", lerr,
")"
1118 CALL cp_warn(__location__, &
1119 'ASSERTION (cond) failed at line '//
cp_to_string(__line__)// &
1120 ' Error '//trim(this_error)//
' in computing numerical derivatives larger then'// &
1121 trim(def_error)//
' .')
1124 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1125 dedf, iforce_eval, nforce_eval, map_index, mapping_section, .false.)
1126 force_env%mixed_env%val(iforce_eval) = energies(iforce_eval)
1129 force_env%mixed_env%dx = dx
1130 force_env%mixed_env%lerr = lerr
1131 force_env%mixed_env%coupling_function = trim(coupling_function)
1138 IF (
SIZE(itmplist) /= 2) &
1139 CALL cp_abort(__location__, &
1140 "Keyword FORCE_STATES takes exactly two input values.")
1141 IF (any(itmplist .LT. 0)) &
1142 cpabort(
"Invalid force_eval index.")
1144 IF (istate(1) > nforce_eval .OR. istate(2) > nforce_eval) &
1145 cpabort(
"Invalid force_eval index.")
1146 mixed_energy%pot = lambda*energies(istate(1)) + (1.0_dp - lambda)*energies(istate(2))
1148 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1149 lambda, istate(1), nforce_eval, map_index, mapping_section, .true.)
1150 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, &
1151 (1.0_dp - lambda), istate(2), nforce_eval, map_index, mapping_section, .false.)
1156 DO iforce_eval = 1, nforce_eval
1157 DEALLOCATE (global_forces(iforce_eval)%forces)
1158 IF (
ASSOCIATED(virials(iforce_eval)%virial))
DEALLOCATE (virials(iforce_eval)%virial)
1161 DEALLOCATE (global_forces)
1162 DEALLOCATE (subsystems)
1163 DEALLOCATE (particles)
1164 DEALLOCATE (energies)
1165 DEALLOCATE (glob_natoms)
1166 DEALLOCATE (virials)
1167 DEALLOCATE (results)
1170 extension=
".data", middle_name=
"MIXED_DIPOLE", log_filename=.false.)
1171 IF (unit_nr > 0)
THEN
1172 description =
'[DIPOLE]'
1173 dip_exists =
test_for_result(results=results_mix, description=description)
1174 IF (dip_exists)
THEN
1175 CALL get_results(results=results_mix, description=description, values=dip_mix)
1176 WRITE (unit_nr,
'(/,1X,A,T48,3F21.16)')
"MIXED ENV| DIPOLE ( A.U.)|", dip_mix
1177 WRITE (unit_nr,
'( 1X,A,T48,3F21.16)')
"MIXED ENV| DIPOLE (Debye)|", dip_mix*
debye
1179 WRITE (unit_nr, *)
"NO FORCE_EVAL section calculated the dipole"
1183 END SUBROUTINE mixed_energy_forces
1198 SUBROUTINE mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
1199 glob_natoms, virials, results)
1201 LOGICAL,
INTENT(IN) :: calculate_forces
1203 REAL(kind=
dp),
DIMENSION(:),
POINTER :: energies
1204 INTEGER,
DIMENSION(:),
POINTER :: glob_natoms
1208 INTEGER :: iforce_eval, iparticle, jparticle, &
1209 my_group, natom, nforce_eval
1210 INTEGER,
DIMENSION(:),
POINTER :: map_index
1211 REAL(kind=
dp) :: energy
1219 mixed_section, root_section
1220 TYPE(
virial_type),
POINTER :: loc_virial, virial_mix
1223 cpassert(
ASSOCIATED(force_env))
1226 subsys=subsys_mix, &
1227 force_env_section=force_env_section, &
1228 root_section=root_section, &
1231 particles=particles_mix, &
1232 virial=virial_mix, &
1233 results=results_mix)
1235 nforce_eval =
SIZE(force_env%sub_force_env)
1238 ALLOCATE (subsystems(nforce_eval))
1239 DO iforce_eval = 1, nforce_eval
1240 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1241 NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
1242 ALLOCATE (virials(iforce_eval)%virial)
1244 IF (.NOT.
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1246 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1247 subsys=subsystems(iforce_eval)%subsys)
1250 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
1254 particles=particles(iforce_eval)%list)
1257 natom =
SIZE(particles(iforce_eval)%list%els)
1259 IF (
ASSOCIATED(map_index)) &
1260 DEALLOCATE (map_index)
1265 DO iparticle = 1, natom
1266 jparticle = map_index(iparticle)
1267 particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
1270 IF (force_env%mixed_env%do_mixed_qmmm_cdft) &
1277 DO iforce_eval = 1, nforce_eval
1278 IF (.NOT.
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1280 IF (force_env%mixed_env%cdft_control%run_type ==
mixed_cdft_serial .AND. iforce_eval .GE. 2)
THEN
1281 my_logger => force_env%mixed_env%cdft_control%sub_logger(iforce_eval - 1)%p
1283 my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
1284 my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
1293 calc_force=calculate_forces, &
1294 skip_external_control=.true.)
1296 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
1297 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1298 potential_energy=energy)
1300 virial=loc_virial, results=loc_results)
1301 energies(iforce_eval) = energy
1302 glob_natoms(iforce_eval) = natom
1303 virials(iforce_eval)%virial = loc_virial
1307 IF (
ASSOCIATED(map_index))
THEN
1308 DEALLOCATE (map_index)
1312 DEALLOCATE (subsystems)
1314 END SUBROUTINE mixed_cdft_energy_forces
1324 SUBROUTINE mixed_cdft_post_energy_forces(force_env)
1327 INTEGER :: iforce_eval, nforce_eval, nvar
1331 cpassert(
ASSOCIATED(force_env))
1332 NULLIFY (qs_env, dft_control)
1333 IF (force_env%mixed_env%do_mixed_cdft)
THEN
1334 nforce_eval =
SIZE(force_env%sub_force_env)
1335 nvar = force_env%mixed_env%cdft_control%nconstraint
1337 IF (.NOT.
ASSOCIATED(force_env%mixed_env%strength)) &
1338 ALLOCATE (force_env%mixed_env%strength(nforce_eval, nvar))
1339 force_env%mixed_env%strength = 0.0_dp
1340 DO iforce_eval = 1, nforce_eval
1341 IF (.NOT.
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1342 IF (force_env%mixed_env%do_mixed_qmmm_cdft)
THEN
1343 qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
1345 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1347 CALL get_qs_env(qs_env, dft_control=dft_control)
1348 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) &
1349 force_env%mixed_env%strength(iforce_eval, :) = dft_control%qs_control%cdft_control%strength(:)
1351 CALL force_env%para_env%sum(force_env%mixed_env%strength)
1353 IF (force_env%mixed_env%do_mixed_et)
THEN
1354 IF (
modulo(force_env%mixed_env%cdft_control%sim_step, force_env%mixed_env%et_freq) == 0) &
1359 END SUBROUTINE mixed_cdft_post_energy_forces
1366 SUBROUTINE embed_energy(force_env)
1370 INTEGER :: iforce_eval, iparticle, jparticle, &
1371 my_group, natom, nforce_eval
1372 INTEGER,
DIMENSION(:),
POINTER :: glob_natoms, map_index
1373 LOGICAL :: converged_embed
1374 REAL(kind=
dp) :: energy
1375 REAL(kind=
dp),
DIMENSION(:),
POINTER :: energies
1389 mapping_section, root_section
1392 cpassert(
ASSOCIATED(force_env))
1395 subsys=subsys_embed, &
1396 force_env_section=force_env_section, &
1397 root_section=root_section, &
1400 particles=particles_embed, &
1401 results=results_embed)
1402 NULLIFY (map_index, glob_natoms)
1404 nforce_eval =
SIZE(force_env%sub_force_env)
1408 ALLOCATE (subsystems(nforce_eval))
1409 ALLOCATE (particles(nforce_eval))
1411 ALLOCATE (energies(nforce_eval))
1412 ALLOCATE (glob_natoms(nforce_eval))
1413 ALLOCATE (results(nforce_eval))
1417 DO iforce_eval = 1, nforce_eval
1418 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
1419 NULLIFY (results(iforce_eval)%results)
1421 IF (.NOT.
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1423 my_group = force_env%embed_env%group_distribution(force_env%para_env%mepos)
1424 my_logger => force_env%embed_env%sub_logger(my_group + 1)%p
1430 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
1431 subsys=subsystems(iforce_eval)%subsys)
1435 IF (
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env%qs_env))
THEN
1436 NULLIFY (dft_control)
1437 CALL get_qs_env(force_env%sub_force_env(iforce_eval)%force_env%qs_env, dft_control=dft_control)
1438 IF (dft_control%qs_control%ref_embed_subsys)
THEN
1439 IF (iforce_eval .EQ. 2) cpabort(
"Density importing force_eval can't be the first.")
1444 CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_embed)
1448 particles=particles(iforce_eval)%list)
1451 natom =
SIZE(particles(iforce_eval)%list%els)
1457 DO iparticle = 1, natom
1458 jparticle = map_index(iparticle)
1459 particles(iforce_eval)%list%els(iparticle)%r = particles_embed%els(jparticle)%r
1464 calc_force=.false., &
1465 skip_external_control=.true.)
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
1473 CALL dft_embedding(force_env, iforce_eval, energies, converged_embed)
1474 IF (.NOT. converged_embed) cpabort(
"Embedding potential optimization not converged.")
1477 IF (dft_control%qs_control%high_level_embed_subsys)
THEN
1478 CALL get_qs_env(qs_env=force_env%sub_force_env(iforce_eval)%force_env%qs_env, &
1479 embed_pot=embed_pot, spin_embed_pot=spin_embed_pot, pw_env=pw_env)
1480 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1481 CALL auxbas_pw_pool%give_back_pw(embed_pot)
1482 IF (
ASSOCIATED(embed_pot))
THEN
1483 CALL embed_pot%release()
1484 DEALLOCATE (embed_pot)
1486 IF (
ASSOCIATED(spin_embed_pot))
THEN
1487 CALL auxbas_pw_pool%give_back_pw(spin_embed_pot)
1488 CALL spin_embed_pot%release()
1489 DEALLOCATE (spin_embed_pot)
1495 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
1496 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
1497 potential_energy=energy)
1499 results=loc_results)
1500 energies(iforce_eval) = energy
1501 glob_natoms(iforce_eval) = natom
1505 IF (
ASSOCIATED(map_index))
THEN
1506 DEALLOCATE (map_index)
1512 CALL force_env%para_env%sync()
1514 CALL force_env%para_env%sum(energies)
1515 CALL force_env%para_env%sum(glob_natoms)
1517 force_env%embed_env%energies = energies
1520 DO iparticle = 1,
SIZE(particles_embed%els)
1521 particles_embed%els(iparticle)%f(:) = 0.0_dp
1525 force_env%embed_env%pot_energy = energies(3) + energies(4) - energies(2)
1528 DO iforce_eval = 1, nforce_eval
1531 DEALLOCATE (subsystems)
1532 DEALLOCATE (particles)
1533 DEALLOCATE (energies)
1534 DEALLOCATE (glob_natoms)
1535 DEALLOCATE (results)
1537 END SUBROUTINE embed_energy
1546 SUBROUTINE dft_embedding(force_env, ref_subsys_number, energies, converged_embed)
1548 INTEGER :: ref_subsys_number
1549 REAL(kind=
dp),
DIMENSION(:),
POINTER :: energies
1550 LOGICAL :: converged_embed
1552 INTEGER :: embed_method
1557 force_env_section=force_env_section)
1561 SELECT CASE (embed_method)
1564 CALL dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1567 CALL dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1570 END SUBROUTINE dft_embedding
1579 SUBROUTINE dfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1581 INTEGER :: ref_subsys_number
1582 REAL(kind=
dp),
DIMENSION(:),
POINTER :: energies
1583 LOGICAL :: converged_embed
1585 CHARACTER(LEN=*),
PARAMETER :: routinen =
'dfet_embedding'
1587 INTEGER :: cluster_subsys_num, handle, &
1588 i_force_eval, i_iter, i_spin, &
1589 nforce_eval, nspins, nspins_subsys, &
1591 REAL(kind=
dp) :: cluster_energy
1592 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rhs
1599 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r_ref, rho_r_subsys
1601 spin_embed_pot, spin_embed_pot_subsys
1605 force_env_section, input, &
1606 mapping_section, opt_embed_section
1608 CALL timeset(routinen, handle)
1613 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
1618 output_unit =
cp_print_key_unit_nr(logger, force_env%force_env_section,
"PRINT%PROGRAM_RUN_INFO", &
1621 NULLIFY (dft_section, input, opt_embed_section)
1622 NULLIFY (energy, dft_control)
1624 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1625 pw_env=pw_env, dft_control=dft_control, rho=rho, energy=energy, &
1627 nspins = dft_control%nspins
1633 CALL qs_rho_get(rho_struct=rho, rho_r=rho_r_ref)
1636 CALL understand_spin_states(force_env, ref_subsys_number, opt_embed%change_spin, opt_embed%open_shell_embed, &
1637 opt_embed%all_nspins)
1640 CALL prepare_embed_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, opt_embed, &
1644 CALL init_embed_pot(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, &
1645 opt_embed%add_const_pot, opt_embed%Fermi_Amaldi, opt_embed%const_pot, &
1646 opt_embed%open_shell_embed, spin_embed_pot, &
1647 opt_embed%pot_diff, opt_embed%Coulomb_guess, opt_embed%grid_opt)
1650 IF (opt_embed%read_embed_pot .OR. opt_embed%read_embed_pot_cube)
CALL read_embed_pot( &
1651 force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, embed_pot, spin_embed_pot, &
1652 opt_embed_section, opt_embed)
1655 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1656 CALL auxbas_pw_pool%create_pw(diff_rho_r)
1658 IF (opt_embed%open_shell_embed)
THEN
1659 CALL auxbas_pw_pool%create_pw(diff_rho_spin)
1664 DO i_spin = 1, nspins
1665 CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
1667 IF (opt_embed%open_shell_embed)
THEN
1668 IF (nspins .EQ. 2)
THEN
1669 CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
1670 CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
1674 DO i_force_eval = 1, ref_subsys_number - 1
1675 NULLIFY (subsys_rho, rho_r_subsys, dft_control)
1676 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, energy=energy, &
1677 dft_control=dft_control)
1678 nspins_subsys = dft_control%nspins
1680 CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
1681 DO i_spin = 1, nspins_subsys
1682 CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
1684 IF (opt_embed%open_shell_embed)
THEN
1685 IF (nspins_subsys .EQ. 2)
THEN
1687 IF ((i_force_eval .EQ. 2) .AND. (opt_embed%change_spin))
THEN
1688 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1689 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1692 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1693 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1700 CALL print_rho_diff(diff_rho_r, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1701 IF (opt_embed%open_shell_embed)
THEN
1702 CALL print_rho_spin_diff(diff_rho_spin, 0, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1706 IF (opt_embed%Coulomb_guess)
THEN
1708 nforce_eval =
SIZE(force_env%sub_force_env)
1710 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, rhs=rhs)
1713 force_env_section=force_env_section)
1717 DO i_force_eval = 1, ref_subsys_number - 1
1718 IF (i_force_eval .EQ. 1)
THEN
1720 force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1722 CALL coulomb_guess(opt_embed%pot_diff, rhs, mapping_section, &
1723 force_env%sub_force_env(i_force_eval)%force_env%qs_env, nforce_eval, i_force_eval, opt_embed%eta)
1726 CALL pw_axpy(opt_embed%pot_diff, embed_pot)
1727 IF (.NOT. opt_embed%grid_opt)
CALL pw_copy(embed_pot, opt_embed%const_pot)
1732 IF (opt_embed%diff_guess)
THEN
1733 CALL pw_copy(diff_rho_r, embed_pot)
1734 IF (.NOT. opt_embed%grid_opt)
CALL pw_copy(embed_pot, opt_embed%const_pot)
1736 IF (opt_embed%open_shell_embed)
CALL pw_copy(diff_rho_spin, spin_embed_pot)
1740 DO i_iter = 1, opt_embed%n_iter
1741 opt_embed%i_iter = i_iter
1745 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, dft_control=dft_control)
1746 nspins = dft_control%nspins
1747 DO i_spin = 1, nspins
1748 CALL pw_axpy(rho_r_ref(i_spin), diff_rho_r, -1.0_dp)
1750 IF (opt_embed%open_shell_embed)
THEN
1752 IF (nspins .EQ. 2)
THEN
1753 CALL pw_axpy(rho_r_ref(1), diff_rho_spin, -1.0_dp)
1754 CALL pw_axpy(rho_r_ref(2), diff_rho_spin, 1.0_dp)
1758 DO i_force_eval = 1, ref_subsys_number - 1
1759 NULLIFY (dft_control)
1760 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
1761 nspins_subsys = dft_control%nspins
1763 IF ((i_force_eval .EQ. 2) .AND. (opt_embed%change_spin))
THEN
1767 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
1768 opt_embed%open_shell_embed, .true.)
1771 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
1772 opt_embed%open_shell_embed, .false.)
1776 dft_control%apply_embed_pot = .true.
1779 CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, embed_pot=embed_pot_subsys)
1780 IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys .EQ. 2))
THEN
1781 CALL set_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
1782 spin_embed_pot=spin_embed_pot_subsys)
1786 CALL get_prev_density(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
1790 calc_force=.false., &
1791 skip_external_control=.true.)
1793 CALL get_max_subsys_diff(opt_embed, force_env%sub_force_env(i_force_eval)%force_env, i_force_eval)
1796 NULLIFY (rho_r_subsys, energy)
1798 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, rho=subsys_rho, &
1800 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + energy%total
1803 IF (dft_control%qs_control%cluster_embed_subsys)
THEN
1804 cluster_subsys_num = i_force_eval
1805 cluster_energy = energy%total
1809 CALL qs_rho_get(rho_struct=subsys_rho, rho_r=rho_r_subsys)
1810 DO i_spin = 1, nspins_subsys
1811 CALL pw_axpy(rho_r_subsys(i_spin), diff_rho_r, allow_noncompatible_grids=.true.)
1813 IF (opt_embed%open_shell_embed)
THEN
1814 IF (nspins_subsys .EQ. 2)
THEN
1816 IF ((i_force_eval .EQ. 2) .AND. (opt_embed%change_spin))
THEN
1817 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1818 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1821 CALL pw_axpy(rho_r_subsys(1), diff_rho_spin, 1.0_dp, allow_noncompatible_grids=.true.)
1822 CALL pw_axpy(rho_r_subsys(2), diff_rho_spin, -1.0_dp, allow_noncompatible_grids=.true.)
1828 CALL embed_pot_subsys%release()
1829 DEALLOCATE (embed_pot_subsys)
1830 IF (opt_embed%open_shell_embed)
THEN
1831 CALL spin_embed_pot_subsys%release()
1832 DEALLOCATE (spin_embed_pot_subsys)
1839 opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
1840 spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .false.)
1842 embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .false., &
1843 force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
1846 DO i_spin = 1, nspins
1847 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) -
pw_integral_ab(embed_pot, rho_r_ref(i_spin))
1850 IF (opt_embed%open_shell_embed)
THEN
1852 IF (nspins .EQ. 2)
THEN
1853 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) &
1859 opt_embed%w_func(i_iter) = opt_embed%w_func(i_iter) + opt_embed%reg_term
1864 IF (opt_embed%converged)
EXIT
1867 IF ((i_iter .GT. 1) .AND. (.NOT. opt_embed%steep_desc))
CALL step_control(opt_embed)
1870 CALL print_rho_diff(diff_rho_r, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1871 IF (opt_embed%open_shell_embed)
THEN
1872 CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .false.)
1877 IF (opt_embed%accept_step .AND. (.NOT. opt_embed%grid_opt)) &
1879 diff_rho_r, diff_rho_spin, opt_embed)
1881 CALL opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, &
1882 force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
1888 opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, &
1889 spin_embed_pot, opt_embed%open_shell_embed, opt_embed%grid_opt, .true.)
1891 embed_pot, spin_embed_pot, i_iter, opt_embed%open_shell_embed, .true., &
1892 force_env%sub_force_env(cluster_subsys_num)%force_env%qs_env)
1896 IF (opt_embed%open_shell_embed)
THEN
1897 CALL print_rho_spin_diff(diff_rho_spin, i_iter, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, .true.)
1901 CALL diff_rho_r%release()
1902 IF (opt_embed%open_shell_embed)
THEN
1903 CALL diff_rho_spin%release()
1907 "PRINT%PROGRAM_RUN_INFO")
1910 IF (opt_embed%converged)
THEN
1911 CALL get_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, dft_control=dft_control, &
1913 nspins_subsys = dft_control%nspins
1914 dft_control%apply_embed_pot = .true.
1917 embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, &
1918 opt_embed%open_shell_embed, opt_embed%change_spin)
1920 IF (opt_embed%Coulomb_guess)
THEN
1921 CALL pw_axpy(opt_embed%pot_diff, embed_pot_subsys, -1.0_dp, allow_noncompatible_grids=.true.)
1924 CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, embed_pot=embed_pot_subsys)
1926 IF ((opt_embed%open_shell_embed) .AND. (nspins_subsys .EQ. 2))
THEN
1927 CALL set_qs_env(force_env%sub_force_env(ref_subsys_number + 1)%force_env%qs_env, &
1928 spin_embed_pot=spin_embed_pot_subsys)
1932 IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source())
THEN
1933 energies(cluster_subsys_num) = cluster_energy
1941 CALL embed_pot%release()
1942 DEALLOCATE (embed_pot)
1943 IF (opt_embed%open_shell_embed)
THEN
1944 CALL spin_embed_pot%release()
1945 DEALLOCATE (spin_embed_pot)
1948 converged_embed = opt_embed%converged
1950 CALL timestop(handle)
1952 END SUBROUTINE dfet_embedding
1962 SUBROUTINE dmfet_embedding(force_env, ref_subsys_number, energies, converged_embed)
1964 INTEGER :: ref_subsys_number
1965 REAL(kind=
dp),
DIMENSION(:),
POINTER :: energies
1966 LOGICAL :: converged_embed
1968 CHARACTER(LEN=*),
PARAMETER :: routinen =
'dmfet_embedding'
1970 INTEGER :: cluster_subsys_num, handle, &
1971 i_force_eval, i_iter, output_unit
1972 LOGICAL :: subsys_open_shell
1973 REAL(kind=
dp) :: cluster_energy
1981 CALL timeset(routinen, handle)
1983 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1989 output_unit =
cp_print_key_unit_nr(logger, force_env%force_env_section,
"PRINT%PROGRAM_RUN_INFO", &
1992 NULLIFY (dft_section, input, opt_dmfet_section)
1995 CALL get_qs_env(qs_env=force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
1996 energy=energy, input=input)
2003 CALL understand_spin_states(force_env, ref_subsys_number, opt_dmfet%change_spin, opt_dmfet%open_shell_embed, &
2004 opt_dmfet%all_nspins)
2007 CALL prepare_dmfet_opt(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2008 opt_dmfet, opt_dmfet_section)
2011 subsys_open_shell =
subsys_spin(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2012 CALL build_full_dm(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, &
2013 opt_dmfet%dm_total, subsys_open_shell, opt_dmfet%open_shell_embed, opt_dmfet%dm_total_beta)
2018 opt_dmfet%dm_diff_beta, para_env)
2020 DO i_force_eval = 1, ref_subsys_number - 1
2023 subsys_open_shell =
subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2025 CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2026 opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2027 opt_dmfet%dm_subsys_beta)
2031 IF (opt_dmfet%open_shell_embed)
CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, &
2032 1.0_dp, opt_dmfet%dm_subsys_beta)
2037 DO i_iter = 1, opt_dmfet%n_iter
2039 opt_dmfet%i_iter = i_iter
2045 opt_dmfet%dm_diff_beta, para_env)
2048 DO i_force_eval = 1, ref_subsys_number - 1
2051 NULLIFY (dft_control)
2052 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, dft_control=dft_control)
2053 dft_control%apply_dmfet_pot = .true.
2057 calc_force=.false., &
2058 skip_external_control=.true.)
2063 CALL get_qs_env(force_env%sub_force_env(i_force_eval)%force_env%qs_env, energy=energy)
2064 opt_dmfet%w_func(i_iter) = opt_dmfet%w_func(i_iter) + energy%total
2067 IF (dft_control%qs_control%cluster_embed_subsys)
THEN
2068 cluster_subsys_num = i_force_eval
2069 cluster_energy = energy%total
2073 subsys_open_shell =
subsys_spin(force_env%sub_force_env(i_force_eval)%force_env%qs_env)
2075 CALL build_full_dm(force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
2076 opt_dmfet%dm_subsys, subsys_open_shell, opt_dmfet%open_shell_embed, &
2077 opt_dmfet%dm_subsys_beta)
2079 IF (opt_dmfet%open_shell_embed)
THEN
2081 IF ((i_force_eval .EQ. 2) .AND. (opt_dmfet%change_spin))
THEN
2086 CALL cp_fm_scale_and_add(-1.0_dp, opt_dmfet%dm_diff_beta, 1.0_dp, opt_dmfet%dm_subsys_beta)
2094 CALL check_dmfet(opt_dmfet, force_env%sub_force_env(ref_subsys_number)%force_env%qs_env)
2099 IF (force_env%sub_force_env(cluster_subsys_num)%force_env%para_env%is_source())
THEN
2100 energies(cluster_subsys_num) = cluster_energy
2105 converged_embed = .false.
2107 END SUBROUTINE dmfet_embedding
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Holds information on atomic properties.
subroutine, public atprop_init(atprop_env, natom)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public huang2011
integer, save, public heaton_burgess2007
Handles all functions related to the CELL.
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Handles all functions related to the CELL.
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
integer, parameter, public cell_sym_triclinic
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
subroutine, public fix_atom_control(force_env, w)
allows for fix atom constraints
Routines to handle the virtual site constraint/restraint.
subroutine, public vsite_force_control(force_env)
control force distribution for virtual sites
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
represent a full matrix distributed on many processors
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Collection of routines to handle the iteration info.
subroutine, public cp_iteration_info_copy_iter(iteration_info_in, iteration_info_out)
Copies iterations info of an iteration info into another iteration info.
various routines to log and control the output. The idea is that decisions about where to log should ...
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
integer, parameter, public low_print_level
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
set of type/routines to handle the storage of results in force_envs
subroutine, public cp_results_mp_bcast(results, source, para_env)
broadcast results type
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
logical function, public test_for_result(results, description)
test for a certain result in the result_list
set of type/routines to handle the storage of results in force_envs
subroutine, public cp_result_copy(results_in, results_out)
Copies the cp_result type.
subroutine, public cp_result_release(results)
Releases cp_result type.
subroutine, public cp_result_create(results)
Allocates and intitializes the cp_result.
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_set(subsys, atomic_kinds, particles, local_particles, molecules, molecule_kinds, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, results, cell)
sets various propreties of the subsys
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
The environment for the empirical interatomic potential methods.
Empirical interatomic potentials for Silicon.
subroutine, public eip_lenosky(eip_env)
Interface routine of Goedecker's Lenosky force field to CP2K.
subroutine, public eip_bazant(eip_env)
Interface routine of Goedecker's Bazant EDIP to CP2K.
Methods to include the effect of an external potential during an MD or energy calculation.
subroutine, public add_external_potential(force_env)
...
subroutine, public fist_calc_energy_force(fist_env, debug)
Calculates the total potential energy, total force, and the total pressure tensor from the potentials...
Interface for the force calculations.
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
subroutine, public force_env_calc_num_pressure(force_env, dx)
Evaluates the stress tensor and pressure numerically.
subroutine, public force_env_create(force_env, root_section, para_env, globenv, fist_env, qs_env, meta_env, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, force_env_section, mixed_env, embed_env, nnp_env, ipi_env)
creates and initializes a force environment
Interface for the force calculations.
integer function, public force_env_get_natom(force_env)
returns the number of atoms
integer, parameter, public use_qmmm
integer, parameter, public use_mixed_force
character(len=10), dimension(501:510), parameter, public use_prog_name
integer, parameter, public use_eip_force
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
integer, parameter, public use_embed
integer, parameter, public use_qmmmx
subroutine, public force_env_set(force_env, meta_env, fp_env, force_env_section, method_name_id, additional_potential)
changes some attributes of the force_env
integer, parameter, public use_qs_force
integer, parameter, public use_pwdft_force
integer, parameter, public use_nnp_force
integer, parameter, public use_ipi
integer, parameter, public use_fist_force
subroutine, public write_forces(particles, iw, label, ndigits, unit_string, total_force, grand_total_force, zero_force_core_shell_atom)
Write forces either to the screen or to a file.
subroutine, public write_atener(iounit, particles, atener, label)
Write the atomic coordinates, types, and energies.
subroutine, public rescale_forces(force_env)
Rescale forces if requested.
subroutine, public get_generic_info(gen_section, func_name, xfunction, parameters, values, var_values, size_variables, i_rep_sec, input_variables)
Reads from the input structure all information for generic functions.
methods used in the flexible partitioning scheme
subroutine, public fp_eval(fp_env, subsys, cell)
Computes the forces and the energy due to the flexible potential & bias, and writes the weights file.
This public domain function parser module is intended for applications where a set of mathematical ex...
subroutine, public parsef(i, funcstr, var)
Parse ith function string FuncStr and compile it into bytecode.
real(rn) function, public evalf(i, val)
...
integer, public evalerrtype
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
subroutine, public finalizef()
...
subroutine, public initf(n)
...
Define type storing the global information of a run. Keep the amount of stored data small....
subroutine, public globenv_retain(globenv)
Retains the global environment globenv.
subroutine, public write_grrm(iounit, force_env, particles, energy, dipole, hessian, dipder, polar, fixed_atoms)
Write GRRM interface file.
The environment for the empirical interatomic potential methods.
i–PI server mode: Communication with i–PI clients
subroutine, public request_forces(ipi_env)
Send atomic positions to a client and retrieve forces.
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Collection of simple mathematical functions and subroutines.
logical function, public abnormal_value(a)
determines if a value is not normal (e.g. for Inf and Nan) based on IO to work also under optimizatio...
Interface to the message passing library MPI.
Methods for mixed CDFT calculations.
subroutine, public mixed_cdft_calculate_coupling(force_env)
Driver routine to calculate the electronic coupling(s) between CDFT states.
subroutine, public mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
Driver routine to handle the build of CDFT weight/gradient in parallel and serial modes.
subroutine, public mixed_cdft_init(force_env, calculate_forces)
Initialize a mixed CDFT calculation.
subroutine, public get_mixed_env(mixed_env, atomic_kind_set, particle_set, local_particles, local_molecules, molecule_kind_set, molecule_set, cell, cell_ref, mixed_energy, para_env, sub_para_env, subsys, input, results, cdft_control)
Get the MIXED environment.
subroutine, public get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, map_index, force_eval_embed)
performs mapping of the subsystems of different force_eval
subroutine, public mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results, factor, iforce_eval, nforce_eval, map_index, mapping_section, overwrite)
Maps forces between the different force_eval sections/environments.
represent a simple array based list of the given type
Define the molecule kind structure types and the corresponding functionality.
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Data types for neural network potentials.
Methods dealing with Neural Network potentials.
subroutine, public nnp_calc_energy_force(nnp, calc_forces)
Calculate the energy and force for a given configuration with the NNP.
subroutine, public check_dmfet(opt_dmfet, qs_env)
...
subroutine, public release_dmfet_opt(opt_dmfet)
...
subroutine, public build_full_dm(qs_env, dm, open_shell, open_shell_embed, dm_beta)
Builds density matrices from MO coefficients in full matrix format.
subroutine, public prepare_dmfet_opt(qs_env, opt_dmfet, opt_dmfet_section)
...
logical function, public subsys_spin(qs_env)
...
subroutine, public print_embed_restart(qs_env, dimen_aux, embed_pot_coef, embed_pot, i_iter, embed_pot_spin, open_shell_embed, grid_opt, final_one)
Print embedding potential as a cube and as a binary (for restarting)
subroutine, public make_subsys_embed_pot(qs_env, embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, open_shell_embed, change_spin_sign)
Creates a subsystem embedding potential.
subroutine, public print_rho_spin_diff(spin_diff_rho_r, i_iter, qs_env, final_one)
Prints a cube for the (spin_rho_A + spin_rho_B - spin_rho_ref) to be minimized in embedding.
subroutine, public get_max_subsys_diff(opt_embed, force_env, subsys_num)
...
subroutine, public print_emb_opt_info(output_unit, step_num, opt_embed)
...
subroutine, public init_embed_pot(qs_env, embed_pot, add_const_pot, fermi_amaldi, const_pot, open_shell_embed, spin_embed_pot, pot_diff, coulomb_guess, grid_opt)
...
subroutine, public understand_spin_states(force_env, ref_subsys_number, change_spin, open_shell_embed, all_nspins)
Find out whether we need to swap alpha- and beta- spind densities in the second subsystem.
subroutine, public read_embed_pot(qs_env, embed_pot, spin_embed_pot, section, opt_embed)
...
subroutine, public get_prev_density(opt_embed, force_env, subsys_num)
...
subroutine, public coulomb_guess(v_rspace, rhs, mapping_section, qs_env, nforce_eval, iforce_eval, eta)
Calculates subsystem Coulomb potential from the RESP charges of the total system.
subroutine, public print_rho_diff(diff_rho_r, i_iter, qs_env, final_one)
Prints a cube for the (rho_A + rho_B - rho_ref) to be minimized in embedding.
subroutine, public conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit)
...
subroutine, public step_control(opt_embed)
Controls the step, changes the trust radius if needed in maximization of the V_emb.
subroutine, public opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, qs_env)
Takes maximization step in embedding potential optimization.
subroutine, public calculate_embed_pot_grad(qs_env, diff_rho_r, diff_rho_spin, opt_embed)
Calculates the derivative of the embedding potential wrt to the expansion coefficients.
subroutine, public print_pot_simple_grid(qs_env, embed_pot, embed_pot_spin, i_iter, open_shell_embed, final_one, qs_env_cluster)
Prints a volumetric file: X Y Z value for interfacing with external programs.
subroutine, public prepare_embed_opt(qs_env, opt_embed, opt_embed_section)
Creates and allocates objects for optimization of embedding potential.
subroutine, public release_opt_embed(opt_embed)
Deallocate stuff for optimizing embedding potential.
represent a simple array based list of the given type
Definition of physical constants:
real(kind=dp), parameter, public debye
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
The type definitions for the PWDFT environment.
Methods and functions on the PWDFT environment.
subroutine, public pwdft_calc_energy_force(pwdft_env, calculate_forces, calculate_stress)
Calculate energy and forces within the PWDFT/SIRIUS code.
Calculates QM/MM energy and forces.
subroutine, public qmmm_calc_energy_force(qmmm_env, calc_force, consistent_energies, linres)
calculates the qm/mm energy and forces
Basic container type for QM/MM.
subroutine, public apply_qmmm_translate(qmmm_env)
Apply translation to the full system in order to center the QM system into the QM box.
Calculates QM/MM energy and forces with Force-Mixing.
subroutine, public qmmmx_calc_energy_force(qmmmx_env, calc_force, consistent_energies, linres, require_consistent_energy_force)
calculates the qm/mm energy and forces
Basic container type for QM/MM with force mixing.
Atomic Polarization Tensor calculation by dF/d(E-field) finite differences.
subroutine, public apt_fdiff(force_env)
Calculate Atomic Polarization Tensors by dF/d(E-field) finite differences.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Quickstep force driver routine.
subroutine, public qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Handles all possible kinds of restraints in CP2K.
subroutine, public restraint_control(force_env)
Computes restraints.
subroutine, public write_scine(iounit, force_env, particles, energy, hessian)
Write SCINE interface file.
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
subroutine, public write_stress_tensor(pv_virial, iw, cell, unit_string, numerical)
Print stress tensor to output file.
subroutine, public write_stress_tensor_components(virial, iw, cell, unit_string)
...
subroutine, public zero_virial(virial, reset)
...
subroutine, public symmetrize_virial(virial)
Symmetrize the virial components.
type for the atomic properties
Type defining parameters related to the simulation cell.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
contains arbitrary information which need to be stored
represent a pointer to a subsys, to be able to create arrays of pointers
represents a system: atoms, molecules, their pos,vel,...
The empirical interatomic potential environment.
Embedding environment type.
Type containing main data for matrix embedding potential optimization.
Type containing main data for embedding potential optimization.
allows for the creation of an array of force_env
wrapper to abstract the force evaluation of the various methods
contains the initially parsed file and the initial parallel environment
stores all the informations relevant to an mpi environment
represent a list of objects
Main data type collecting all relevant data for neural network potentials.
represents a pointer to a list
represent a list of objects
contained for different pw related things
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
The PWDFT environment type.
keeps the density in various representations, keeping track of which ones are valid.