68#include "./base/base_uses.f90"
73 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .true.
74 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qmmm_image_charge'
103 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_image_pot'
107 CALL timeset(routinen, handle)
109 IF (qmmm_env%image_charge_pot%coeff_iterative)
THEN
111 IF (qs_env%calc_image_preconditioner)
THEN
112 IF (qmmm_env%image_charge_pot%image_restart)
THEN
113 CALL restart_image_matrix(image_matrix=qs_env%image_matrix, &
114 qs_env=qs_env, qmmm_env=qmmm_env)
116 CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
117 qs_env=qs_env, qmmm_env=qmmm_env)
120 CALL calc_image_coeff_iterative(v_hartree_rspace=v_hartree_rspace, &
121 coeff=qs_env%image_coeff, qmmm_env=qmmm_env, &
125 CALL calc_image_coeff_gaussalgorithm(v_hartree_rspace=v_hartree_rspace, &
126 coeff=qs_env%image_coeff, qmmm_env=qmmm_env, &
131 ALLOCATE (qs_env%ks_qmmm_env%v_metal_rspace)
132 CALL calculate_potential_metal(v_metal_rspace= &
133 qs_env%ks_qmmm_env%v_metal_rspace, coeff=qs_env%image_coeff, &
134 rho_hartree_gspace=rho_hartree_gspace, &
135 energy=energy, qs_env=qs_env)
137 CALL timestop(handle)
150 SUBROUTINE calc_image_coeff_gaussalgorithm(v_hartree_rspace, coeff, qmmm_env, &
154 REAL(kind=
dp),
DIMENSION(:),
POINTER :: coeff
158 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_image_coeff_gaussalgorithm'
160 INTEGER :: handle, info, natom
161 REAL(kind=
dp) :: eta, v0
162 REAL(kind=
dp),
DIMENSION(:),
POINTER :: pot_const
164 CALL timeset(routinen, handle)
169 v0 = -qmmm_env%image_charge_pot%V0
170 eta = qmmm_env%image_charge_pot%eta
171 natom =
SIZE(qmmm_env%image_charge_pot%image_mm_list)
173 ALLOCATE (pot_const(natom))
174 IF (.NOT.
ASSOCIATED(coeff))
THEN
175 ALLOCATE (coeff(natom))
179 CALL integrate_potential_ga_rspace(v_hartree_rspace, qmmm_env, qs_env, &
182 pot_const(:) = -pot_const(:) + v0*sqrt((
pi/eta)**3)
186 CALL dgetrs(
'N', natom, 1, qs_env%image_matrix, natom, qs_env%ipiv, &
187 pot_const, natom, info)
192 DEALLOCATE (pot_const)
194 CALL timestop(handle)
196 END SUBROUTINE calc_image_coeff_gaussalgorithm
206 SUBROUTINE calc_image_coeff_iterative(v_hartree_rspace, coeff, qmmm_env, &
210 REAL(kind=
dp),
DIMENSION(:),
POINTER :: coeff
214 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_image_coeff_iterative'
216 INTEGER :: handle, iter_steps, natom, output_unit
217 REAL(kind=
dp) :: alpha, eta, rsnew, rsold, v0
218 REAL(kind=
dp),
DIMENSION(:),
POINTER :: ad, d, pot_const, r, vmetal_const, z
223 CALL timeset(routinen, handle)
225 NULLIFY (pot_const, vmetal_const, logger, input)
229 v0 = -qmmm_env%image_charge_pot%V0
230 eta = qmmm_env%image_charge_pot%eta
231 natom =
SIZE(qmmm_env%image_charge_pot%image_mm_list)
233 ALLOCATE (pot_const(natom))
234 ALLOCATE (vmetal_const(natom))
239 IF (.NOT.
ASSOCIATED(coeff))
THEN
240 ALLOCATE (coeff(natom))
243 CALL integrate_potential_ga_rspace(v_hartree_rspace, qmmm_env, qs_env, &
247 pot_const(:) = -pot_const(:) + v0*sqrt((
pi/eta)**3)
259 CALL calculate_potential_metal(v_metal_rspace=v_metal_rspace_guess, &
260 coeff=coeff, qs_env=qs_env)
261 CALL integrate_potential_ga_rspace(potential=v_metal_rspace_guess, &
262 qmmm_env=qmmm_env, qs_env=qs_env, int_res=vmetal_const)
265 r = pot_const - vmetal_const
266 z = matmul(qs_env%image_matrix, r)
268 rsold = dot_product(r, z)
273 CALL calculate_potential_metal(v_metal_rspace= &
274 auxpot_ad_rspace, coeff=d, qs_env=qs_env)
275 CALL integrate_potential_ga_rspace(potential= &
276 auxpot_ad_rspace, qmmm_env=qmmm_env, &
277 qs_env=qs_env, int_res=ad)
279 alpha = rsold/dot_product(d, ad)
280 coeff = coeff + alpha*d
283 z = matmul(qs_env%image_matrix, r)
284 rsnew = dot_product(r, z)
285 iter_steps = iter_steps + 1
287 IF (rsnew < 1.0e-16)
THEN
288 CALL auxpot_ad_rspace%release()
291 d = z + rsnew/rsold*d
293 CALL auxpot_ad_rspace%release()
300 "QMMM%PRINT%PROGRAM_RUN_INFO", &
301 extension=
".qmmmLog")
302 IF (output_unit > 0)
WRITE (unit=output_unit, fmt=
"(T3,A,T74,I7)") &
303 "Number of iteration steps for determination of image coefficients:", iter_steps
305 "QMMM%PRINT%PROGRAM_RUN_INFO")
307 IF (iter_steps .LT. 25)
THEN
308 qs_env%calc_image_preconditioner = .false.
310 qs_env%calc_image_preconditioner = .true.
313 CALL v_metal_rspace_guess%release()
314 DEALLOCATE (pot_const)
315 DEALLOCATE (vmetal_const)
320 CALL timestop(handle)
322 END SUBROUTINE calc_image_coeff_iterative
334 SUBROUTINE integrate_potential_ga_rspace(potential, qmmm_env, qs_env, int_res, &
335 atom_num, atom_num_ref)
340 REAL(kind=
dp),
DIMENSION(:),
POINTER :: int_res
341 INTEGER,
INTENT(IN),
OPTIONAL :: atom_num, atom_num_ref
343 CHARACTER(LEN=*),
PARAMETER :: routinen =
'integrate_potential_ga_rspace'
345 INTEGER :: atom_a, atom_b, atom_ref, handle, iatom, &
347 INTEGER,
DIMENSION(:),
POINTER :: cores
348 REAL(kind=
dp) :: eps_rho_rspace, radius
349 REAL(kind=
dp),
DIMENSION(3) :: ra
350 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: hab
358 CALL timeset(routinen, handle)
360 NULLIFY (cores, hab, cell, auxbas_rs_desc, pw_env, para_env, &
365 CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
371 dft_control=dft_control, &
372 para_env=para_env, pw_env=pw_env)
374 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
376 natom =
SIZE(qmmm_env%image_charge_pot%image_mm_list)
378 IF (
PRESENT(atom_num)) k = atom_num
386 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed)
THEN
388 IF (
modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos)
THEN
401 atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
403 IF (
PRESENT(atom_num) .AND.
PRESENT(atom_num_ref))
THEN
405 atom_b = qmmm_env%image_charge_pot%image_mm_list(k)
406 atom_ref = qmmm_env%image_charge_pot%image_mm_list(atom_num_ref)
407 ra(:) =
pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell) &
408 -
pbc(qmmm_env%image_charge_pot%particles_all(atom_b)%r, cell) &
409 +
pbc(qmmm_env%image_charge_pot%particles_all(atom_ref)%r, cell)
412 ra(:) =
pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell)
418 ra=ra, rb=ra, rp=ra, &
419 zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, &
420 prefactor=1.0_dp, cutoff=1.0_dp)
422 CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, &
423 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
424 rs_v, hab, o1=0, o2=0, &
425 radius=radius, calculate_forces=.false., &
426 use_subpatch=.true., subpatch_pattern=0)
428 int_res(iatom) = hab(1, 1)
432 CALL para_env%sum(int_res)
434 DEALLOCATE (hab, cores)
436 CALL timestop(handle)
438 END SUBROUTINE integrate_potential_ga_rspace
454 REAL(kind=
dp),
DIMENSION(:),
POINTER :: coeff
455 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: forces
459 CHARACTER(LEN=*),
PARAMETER :: routinen =
'integrate_potential_devga_rspace'
461 INTEGER :: atom_a, handle, iatom, j, natom, npme
462 INTEGER,
DIMENSION(:),
POINTER :: cores
463 LOGICAL :: use_virial
464 REAL(kind=
dp) :: eps_rho_rspace, radius
465 REAL(kind=
dp),
DIMENSION(3) :: force_a, force_b, ra
466 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: hab, pab
475 CALL timeset(routinen, handle)
477 NULLIFY (cores, hab, pab, cell, auxbas_rs_desc, pw_env, para_env, &
478 dft_control, rs_v, virial)
485 CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
491 dft_control=dft_control, &
492 para_env=para_env, pw_env=pw_env, &
495 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
498 cpabort(
"Virial not implemented for image charge method")
501 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
503 natom =
SIZE(qmmm_env%image_charge_pot%image_mm_list)
505 IF (.NOT.
ASSOCIATED(forces))
THEN
506 ALLOCATE (forces(3, natom))
509 forces(:, :) = 0.0_dp
516 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed)
THEN
518 IF (
modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos)
THEN
531 atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
532 ra(:) =
pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell)
539 ra=ra, rb=ra, rp=ra, &
540 zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, &
541 pab=pab, o1=0, o2=0, &
542 prefactor=1.0_dp, cutoff=1.0_dp)
544 CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, &
545 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
546 rs_v, hab, pab, o1=0, o2=0, &
547 radius=radius, calculate_forces=.true., &
548 force_a=force_a, force_b=force_b, use_subpatch=.true., &
551 force_a(:) = coeff(iatom)*force_a(:)
552 forces(:, iatom) = force_a(:)
556 CALL para_env%sum(forces)
558 DEALLOCATE (hab, pab, cores)
561 CALL print_gradients_image_atoms(forces, qs_env)
563 CALL timestop(handle)
578 IF (.NOT. qmmm_env%image_charge_pot%coeff_iterative)
THEN
579 SELECT CASE (qmmm_env%image_charge_pot%state_image_matrix)
581 CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
582 ipiv=qs_env%ipiv, qs_env=qs_env, qmmm_env=qmmm_env)
586 CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
587 ipiv=qs_env%ipiv, qs_env=qs_env, qmmm_env=qmmm_env)
589 IF (qmmm_env%center_qm_subsys0) &
590 CALL cp_warn(__location__, &
591 "The image atoms are fully "// &
592 "constrained and the image matrix is only calculated once. "// &
593 "To be safe, set CENTER to NEVER ")
597 cpabort(
"No initialization for image charges available?")
610 SUBROUTINE calculate_image_matrix(image_matrix, ipiv, qs_env, qmmm_env)
612 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: image_matrix
613 INTEGER,
DIMENSION(:),
OPTIONAL,
POINTER :: ipiv
617 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_image_matrix'
619 INTEGER :: handle, j, k, natom, output_unit, stat
623 CALL timeset(routinen, handle)
624 NULLIFY (input, logger)
628 natom =
SIZE(qmmm_env%image_charge_pot%image_mm_list)
630 IF (.NOT.
ASSOCIATED(image_matrix))
THEN
631 ALLOCATE (image_matrix(natom, natom))
633 IF (
PRESENT(ipiv))
THEN
634 IF (.NOT.
ASSOCIATED(ipiv))
THEN
635 ALLOCATE (ipiv(natom))
643 "QMMM%PRINT%PROGRAM_RUN_INFO", &
644 extension=
".qmmmLog")
645 IF (qmmm_env%image_charge_pot%coeff_iterative)
THEN
646 IF (output_unit > 0)
WRITE (unit=output_unit, fmt=
"(T3,A)") &
647 "Calculating image matrix"
649 IF (output_unit > 0)
WRITE (unit=output_unit, fmt=
"(T2,A)") &
650 "Calculating image matrix"
653 "QMMM%PRINT%PROGRAM_RUN_INFO")
656 SELECT CASE (qmmm_env%image_charge_pot%image_matrix_method)
658 CALL calculate_image_matrix_gpw(image_matrix, qs_env, qmmm_env)
660 CALL calculate_image_matrix_mme(image_matrix, qs_env, qmmm_env)
662 cpabort(
"Unknown method for calculating image matrix")
665 IF (qmmm_env%image_charge_pot%coeff_iterative)
THEN
667 CALL dpotrf(
'L', natom, qs_env%image_matrix, natom, stat)
669 CALL dpotri(
'L', natom, qs_env%image_matrix, natom, stat)
673 qs_env%image_matrix(j, k) = qs_env%image_matrix(k, j)
676 CALL write_image_matrix(qs_env%image_matrix, qs_env)
679 IF (
PRESENT(ipiv))
THEN
680 CALL dgetrf(natom, natom, image_matrix, natom, ipiv, stat)
685 CALL timestop(handle)
687 END SUBROUTINE calculate_image_matrix
695 SUBROUTINE calculate_image_matrix_gpw(image_matrix, qs_env, qmmm_env)
696 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: image_matrix
700 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_image_matrix_gpw'
702 INTEGER :: handle, iatom, iatom_ref, natom
703 REAL(kind=
dp),
DIMENSION(:),
POINTER :: int_res
711 CALL timeset(routinen, handle)
712 NULLIFY (pw_env, auxbas_pw_pool, poisson_env, para_env, int_res)
714 natom =
SIZE(qmmm_env%image_charge_pot%image_mm_list)
715 ALLOCATE (int_res(natom))
717 image_matrix = 0.0_dp
719 CALL get_qs_env(qs_env, pw_env=pw_env, para_env=para_env)
721 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
722 poisson_env=poisson_env)
723 CALL auxbas_pw_pool%create_pw(rho_gb)
724 CALL auxbas_pw_pool%create_pw(vb_gspace)
725 CALL auxbas_pw_pool%create_pw(vb_rspace)
737 CALL pw_scale(vb_rspace, vb_rspace%pw_grid%dvol)
742 CALL integrate_potential_ga_rspace(vb_rspace, qs_env%qmmm_env_qm, &
743 qs_env, int_res, atom_num=iatom, &
744 atom_num_ref=iatom_ref)
745 image_matrix(iatom, iatom:natom) = int_res(iatom:natom)
746 image_matrix(iatom + 1:natom, iatom) = int_res(iatom + 1:natom)
749 CALL vb_gspace%release()
750 CALL vb_rspace%release()
751 CALL rho_gb%release()
755 CALL timestop(handle)
756 END SUBROUTINE calculate_image_matrix_gpw
764 SUBROUTINE calculate_image_matrix_mme(image_matrix, qs_env, qmmm_env)
765 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: image_matrix
769 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_image_matrix_mme'
771 INTEGER :: atom_a, handle, iatom, natom
772 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: zeta
773 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ra
776 CALL timeset(routinen, handle)
779 natom =
SIZE(qmmm_env%image_charge_pot%image_mm_list)
780 ALLOCATE (zeta(natom), ra(3, natom))
782 zeta(:) = qmmm_env%image_charge_pot%eta
785 atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
786 ra(:, iatom) = qmmm_env%image_charge_pot%particles_all(atom_a)%r(:)
791 CALL integrate_s_mme(qmmm_env%image_charge_pot%eri_mme_param, &
792 zeta, zeta, ra, ra, image_matrix, para_env)
794 CALL timestop(handle)
795 END SUBROUTINE calculate_image_matrix_mme
808 SUBROUTINE integrate_s_mme(param, zeta, zetb, ra, rb, hab, para_env)
810 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: zeta, zetb
811 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: ra, rb
812 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: hab
815 CHARACTER(len=*),
PARAMETER :: routinen =
'integrate_s_mme'
817 INTEGER :: g_count, handle, ipgf, ipgf_prod, jpgf, &
818 npgf_prod, npgfa, npgfb, r_count
819 INTEGER,
DIMENSION(2) :: limits
820 REAL(kind=
dp),
DIMENSION(3) :: rab
822 CALL timeset(routinen, handle)
823 g_count = 0; r_count = 0
829 npgf_prod = npgfa*npgfb
831 limits =
get_limit(npgf_prod, para_env%num_pe, para_env%mepos)
833 DO ipgf_prod = limits(1), limits(2)
834 ipgf = (ipgf_prod - 1)/npgfb + 1
835 jpgf = mod(ipgf_prod - 1, npgfb) + 1
836 rab(:) = ra(:, ipgf) - rb(:, jpgf)
838 zetb(jpgf), rab, hab, ipgf - 1, jpgf - 1, g_count=g_count, r_count=r_count)
842 CALL para_env%sum(hab)
843 CALL timestop(handle)
845 END SUBROUTINE integrate_s_mme
857 SUBROUTINE calculate_potential_metal(v_metal_rspace, coeff, rho_hartree_gspace, energy, &
861 REAL(kind=
dp),
DIMENSION(:),
POINTER :: coeff
866 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_potential_metal'
869 REAL(kind=
dp) :: en_external, en_vmetal_rhohartree, &
876 CALL timeset(routinen, handle)
878 NULLIFY (pw_env, auxbas_pw_pool, poisson_env)
879 en_vmetal_rhohartree = 0.0_dp
883 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
884 poisson_env=poisson_env)
886 CALL auxbas_pw_pool%create_pw(rho_metal)
888 CALL auxbas_pw_pool%create_pw(v_metal_gspace)
890 CALL auxbas_pw_pool%create_pw(v_metal_rspace)
898 vhartree=v_metal_gspace)
900 IF (
PRESENT(rho_hartree_gspace))
THEN
903 en_external = qs_env%qmmm_env_qm%image_charge_pot%V0*total_rho_metal
904 energy%image_charge = en_vmetal_rhohartree - 0.5_dp*en_external
905 CALL print_image_energy_terms(en_vmetal_rhohartree, en_external, &
906 total_rho_metal, qs_env)
911 CALL pw_scale(v_metal_rspace, v_metal_rspace%pw_grid%dvol)
912 CALL v_metal_gspace%release()
913 CALL rho_metal%release()
915 CALL timestop(handle)
917 END SUBROUTINE calculate_potential_metal
931 CHARACTER(len=*),
PARAMETER :: routinen =
'add_image_pot_to_hartree_pot'
933 INTEGER :: handle, output_unit
937 CALL timeset(routinen, handle)
939 NULLIFY (input, logger)
943 CALL pw_axpy(v_metal, v_hartree)
949 "QMMM%PRINT%PROGRAM_RUN_INFO", &
950 extension=
".qmmmLog")
951 IF (output_unit > 0)
WRITE (unit=output_unit, fmt=
"(T3,A)") &
952 "Adding image charge potential to the Hartree potential."
954 "QMMM%PRINT%PROGRAM_RUN_INFO")
956 CALL timestop(handle)
966 SUBROUTINE write_image_matrix(image_matrix, qs_env)
968 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: image_matrix
971 CHARACTER(LEN=*),
PARAMETER :: routinen =
'write_image_matrix'
973 CHARACTER(LEN=default_path_length) :: filename
974 INTEGER :: handle, rst_unit
979 CALL timeset(routinen, handle)
981 NULLIFY (qmmm_section, print_key, logger, para_env)
985 CALL get_qs_env(qs_env=qs_env, para_env=para_env, &
989 "QMMM%PRINT%IMAGE_CHARGE_RESTART")
992 qmmm_section,
"QMMM%PRINT%IMAGE_CHARGE_RESTART"), &
996 "QMMM%PRINT%IMAGE_CHARGE_RESTART", &
997 extension=
".Image", &
998 file_status=
"REPLACE", &
999 file_action=
"WRITE", &
1000 file_form=
"UNFORMATTED")
1003 print_key, extension=
".IMAGE", &
1006 IF (rst_unit > 0)
THEN
1007 WRITE (rst_unit) image_matrix
1011 "QMMM%PRINT%IMAGE_CHARGE_RESTART")
1014 CALL timestop(handle)
1016 END SUBROUTINE write_image_matrix
1025 SUBROUTINE restart_image_matrix(image_matrix, qs_env, qmmm_env)
1027 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: image_matrix
1031 CHARACTER(LEN=*),
PARAMETER :: routinen =
'restart_image_matrix'
1033 CHARACTER(LEN=default_path_length) :: image_filename
1034 INTEGER :: handle, natom, output_unit, rst_unit
1040 CALL timeset(routinen, handle)
1042 NULLIFY (qmmm_section, logger, para_env)
1047 natom =
SIZE(qmmm_env%image_charge_pot%image_mm_list)
1049 IF (.NOT.
ASSOCIATED(image_matrix))
THEN
1050 ALLOCATE (image_matrix(natom, natom))
1053 image_matrix = 0.0_dp
1055 CALL get_qs_env(qs_env=qs_env, para_env=para_env, &
1059 c_val=image_filename)
1061 INQUIRE (file=image_filename, exist=exist)
1064 IF (para_env%is_source())
THEN
1065 CALL open_file(file_name=image_filename, &
1066 file_status=
"OLD", &
1067 file_form=
"UNFORMATTED", &
1068 file_action=
"READ", &
1069 unit_number=rst_unit)
1071 READ (rst_unit) qs_env%image_matrix
1074 CALL para_env%bcast(qs_env%image_matrix)
1076 IF (para_env%is_source())
CALL close_file(unit_number=rst_unit)
1079 "QMMM%PRINT%PROGRAM_RUN_INFO", &
1080 extension=
".qmmmLog")
1081 IF (output_unit > 0)
WRITE (unit=output_unit, fmt=
"(T3,A)") &
1082 "Restarted image matrix"
1084 cpabort(
"Restart file for image matrix not found")
1087 qmmm_env%image_charge_pot%image_restart = .false.
1089 CALL timestop(handle)
1091 END SUBROUTINE restart_image_matrix
1099 SUBROUTINE print_gradients_image_atoms(forces, qs_env)
1101 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: forces
1104 INTEGER :: atom_a, iatom, natom, output_unit
1105 REAL(kind=
dp),
DIMENSION(3) :: sum_gradients
1109 NULLIFY (input, logger)
1112 sum_gradients = 0.0_dp
1113 natom =
SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
1116 sum_gradients(:) = sum_gradients(:) + forces(:, iatom)
1122 "QMMM%PRINT%DERIVATIVES", extension=
".Log")
1123 IF (output_unit > 0)
THEN
1124 WRITE (unit=output_unit, fmt=
"(/1X,A)") &
1125 "Image gradients [a.u.] on MM image charge atoms after QMMM calculation: "
1126 WRITE (unit=output_unit, fmt=
"(T4,A4,T27,A1,T50,A1,T74,A1)") &
1127 "Atom",
"X",
"Y",
"Z"
1129 atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1130 WRITE (unit=output_unit, fmt=
"(T2,I6,T22,ES12.5,T45,ES12.5,T69,ES12.5)") &
1131 atom_a, forces(:, iatom)
1134 WRITE (unit=output_unit, fmt=
"(T2,A)") repeat(
"-", 79)
1135 WRITE (unit=output_unit, fmt=
"(T2,A,T22,ES12.5,T45,ES12.5,T69,ES12.5)") &
1136 "sum gradients:", sum_gradients
1137 WRITE (unit=output_unit, fmt=
"(/)")
1141 "QMMM%PRINT%DERIVATIVES")
1143 END SUBROUTINE print_gradients_image_atoms
1152 REAL(kind=
dp),
DIMENSION(:),
POINTER :: image_coeff
1155 INTEGER :: atom_a, iatom, natom, output_unit
1156 REAL(kind=
dp) :: normalize_factor, sum_coeff
1160 NULLIFY (input, logger)
1164 natom =
SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
1165 normalize_factor = sqrt((qs_env%qmmm_env_qm%image_charge_pot%eta/
pi)**3)
1168 sum_coeff = sum_coeff + image_coeff(iatom)
1174 "QMMM%PRINT%IMAGE_CHARGE_INFO", extension=
".Log")
1175 IF (output_unit > 0)
THEN
1176 WRITE (unit=output_unit, fmt=
"(/)")
1177 WRITE (unit=output_unit, fmt=
"(T2,A)") &
1178 "Image charges [a.u.] after QMMM calculation: "
1179 WRITE (unit=output_unit, fmt=
"(T4,A4,T67,A)")
"Atom",
"Image charge"
1180 WRITE (unit=output_unit, fmt=
"(T4,A,T67,A)") repeat(
"-", 4), repeat(
"-", 12)
1183 atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1187 WRITE (unit=output_unit, fmt=
"(T2,I6,T65,ES16.9)") &
1188 atom_a, -image_coeff(iatom)/normalize_factor
1191 WRITE (unit=output_unit, fmt=
"(T2,A)") repeat(
"-", 79)
1192 WRITE (unit=output_unit, fmt=
"(T2,A,T65,ES16.9)") &
1193 "sum image charges:", -sum_coeff/normalize_factor
1197 "QMMM%PRINT%IMAGE_CHARGE_INFO")
1210 SUBROUTINE print_image_energy_terms(en_vmetal_rhohartree, en_external, &
1211 total_rho_metal, qs_env)
1213 REAL(kind=
dp),
INTENT(IN) :: en_vmetal_rhohartree, en_external, &
1217 INTEGER :: output_unit
1221 NULLIFY (input, logger)
1227 "QMMM%PRINT%IMAGE_CHARGE_INFO", extension=
".Log")
1229 IF (output_unit > 0)
THEN
1230 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
1231 "Total induced charge density [a.u.]:", total_rho_metal
1232 WRITE (unit=output_unit, fmt=
"(T3,A)")
"Image energy terms: "
1233 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
1234 "Coulomb energy of QM and image charge density [a.u.]:", en_vmetal_rhohartree
1235 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
1236 "External potential energy term [a.u.]:", -0.5_dp*en_external
1237 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
1238 "Total image charge energy [a.u.]:", en_vmetal_rhohartree - 0.5_dp*en_external
1242 "QMMM%PRINT%IMAGE_CHARGE_INFO")
1244 END SUBROUTINE print_image_energy_terms
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
All kind of helpful little routines.
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Handles all functions related to the CELL.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
subroutine, public cp_eri_mme_update_local_counts(param, para_env, g_count_2c, r_count_2c, gg_count_3c, gr_count_3c, rr_count_3c)
Update local counters to gather statistics on different paths taken in MME algorithm (each Ewald sum ...
Utility routines to open and close files. Tracking of preconnections.
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
various routines to log and control the output. The idea is that decisions about where to log should ...
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)
...
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
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...
Minimax-Ewald (MME) method for calculating 2-center and 3-center electron repulsion integrals (ERI) o...
subroutine, public eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, hab, o1, o2, g_count, r_count, normalize, potential, pot_par)
Low-level integration routine for 2-center ERIs.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_path_length
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Utility routines for the memory handling.
Interface to the message passing library MPI.
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
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Routines for image charge calculation within QM/MM.
subroutine, public calculate_image_pot(v_hartree_rspace, rho_hartree_gspace, energy, qmmm_env, qs_env)
determines coefficients by solving image_matrix*coeff=-pot_const by Gaussian elimination or in an ite...
subroutine, public print_image_coefficients(image_coeff, qs_env)
Print image coefficients.
subroutine, public integrate_potential_devga_rspace(potential, coeff, forces, qmmm_env, qs_env)
calculates the image forces on the MM atoms
subroutine, public conditional_calc_image_matrix(qs_env, qmmm_env)
calculate image matrix T depending on constraints on image atoms in case coefficients are estimated n...
subroutine, public add_image_pot_to_hartree_pot(v_hartree, v_metal, qs_env)
Add potential of metal (image charge pot) to Hartree Potential.
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_metal(rho_metal, coeff, total_rho_metal, qs_env)
computes the image charge density on the grid (including coeffcients)
subroutine, public calculate_rho_single_gaussian(rho_gb, qs_env, iatom_in)
collocate a single Gaussian on the grid
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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Integrate single or product functions over a potential on a RS grid.
subroutine, public transfer_pw2rs(rs, pw)
...
All kind of helpful little routines.
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
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...
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...