80#include "./base/base_uses.f90"
86 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'almo_scf_qs'
115 name_new, size_keys, symmetry_new, &
116 spin_key, init_domains)
120 CHARACTER(len=*),
INTENT(IN) :: name_new
121 INTEGER,
DIMENSION(2),
INTENT(IN) :: size_keys
122 CHARACTER,
INTENT(IN) :: symmetry_new
123 INTEGER,
INTENT(IN) :: spin_key
124 LOGICAL,
INTENT(IN) :: init_domains
126 CHARACTER(len=*),
PARAMETER :: routinen =
'matrix_almo_create'
128 INTEGER :: dimen, handle, hold, iatom, iblock_col, &
129 iblock_row, imol, mynode, natoms, &
130 nblkrows_tot, nlength, nmols, row
131 INTEGER,
DIMENSION(:),
POINTER :: blk_distr, blk_sizes, block_sizes_new, col_blk_size, &
132 col_distr_new, col_sizes_new, distr_new_array, row_blk_size, row_distr_new, row_sizes_new
133 LOGICAL :: active, one_dim_is_mo, tr
134 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: new_block
149 CALL timeset(routinen, handle)
164 natoms = almo_scf_env%natoms
165 nmols = almo_scf_env%nmolecules
183 ALLOCATE (block_sizes_new(natoms), distr_new_array(natoms))
184 block_sizes_new(:) = blk_sizes(:)
185 distr_new_array(:) = blk_distr(:)
188 ALLOCATE (block_sizes_new(nmols), distr_new_array(nmols))
189 block_sizes_new(:) = 0
191 block_sizes_new(almo_scf_env%domain_index_of_atom(iatom)) = &
192 block_sizes_new(almo_scf_env%domain_index_of_atom(iatom)) + &
196 distr_new_array(imol) = &
197 blk_distr(almo_scf_env%first_atom_of_domain(imol))
200 cpabort(
"Illegal distribution")
213 ALLOCATE (block_sizes_new(nlength))
214 block_sizes_new(:) = 0
225 ALLOCATE (block_sizes_new(nlength))
227 block_sizes_new(:) = almo_scf_env%nocc_of_domain(:, spin_key)
230 WHERE (block_sizes_new == 0) block_sizes_new = 1
232 block_sizes_new(:) = almo_scf_env%nvirt_disc_of_domain(:, spin_key)
234 block_sizes_new(:) = almo_scf_env%nvirt_full_of_domain(:, spin_key)
236 block_sizes_new(:) = almo_scf_env%nvirt_of_domain(:, spin_key)
239 cpabort(
"Illegal distribution")
244 cpabort(
"Illegal dimension")
249 ALLOCATE (distr_new_array(nlength))
252 distr_new_array(:) = blk_distr(:)
256 distr_new_array(imol) = &
257 blk_distr(almo_scf_env%first_atom_of_domain(imol))
264 row_sizes_new => block_sizes_new
265 row_distr_new => distr_new_array
267 col_sizes_new => block_sizes_new
268 col_distr_new => distr_new_array
274 row_dist=row_distr_new, col_dist=col_distr_new, &
279 dist_new, symmetry_new, &
280 row_sizes_new, col_sizes_new, reuse_arrays=.true.)
285 IF (init_domains)
THEN
290 row_blk_size=row_blk_size, col_blk_size=col_blk_size)
293 DO row = 1, nblkrows_tot
299 IF (hold == mynode)
THEN
303 one_dim_is_mo = .false.
307 IF (one_dim_is_mo)
THEN
308 IF (almo_scf_env%nocc_of_domain(row, spin_key) == 0) active = .false.
311 one_dim_is_mo = .false.
315 IF (one_dim_is_mo)
THEN
316 IF (almo_scf_env%nvirt_of_domain(row, spin_key) == 0) active = .false.
319 one_dim_is_mo = .false.
323 IF (one_dim_is_mo)
THEN
324 IF (almo_scf_env%nvirt_disc_of_domain(row, spin_key) == 0) active = .false.
327 one_dim_is_mo = .false.
331 IF (one_dim_is_mo)
THEN
332 IF (almo_scf_env%nvirt_full_of_domain(row, spin_key) == 0) active = .false.
336 ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
337 new_block(:, :) = 1.0_dp
339 DEALLOCATE (new_block)
350 CALL timestop(handle)
366 INTEGER :: mat_distr_aos
368 CHARACTER(len=*),
PARAMETER :: routinen =
'matrix_qs_to_almo'
373 CALL timeset(routinen, handle)
376 SELECT CASE (mat_distr_aos)
382 CALL dbcsr_create(matrix_qs_nosym, template=matrix_qs, matrix_type=dbcsr_type_no_symmetry)
397 cpabort(
"Unknown mat_distr_aos for matrix_qs_to_almo")
400 CALL timestop(handle)
413 SUBROUTINE matrix_almo_to_qs(matrix_almo, matrix_qs, mat_distr_aos)
415 INTEGER,
INTENT(IN) :: mat_distr_aos
417 CHARACTER(len=*),
PARAMETER :: routinen =
'matrix_almo_to_qs'
422 CALL timeset(routinen, handle)
425 SELECT CASE (mat_distr_aos)
427 CALL dbcsr_copy(matrix_qs, matrix_almo, keep_sparsity=.true.)
429 CALL dbcsr_create(matrix_almo_redist, template=matrix_qs)
432 CALL dbcsr_copy(matrix_qs, matrix_almo_redist, keep_sparsity=.true.)
435 cpabort(
"Unknown mat_distr_aos for matrix_almo_to_qs")
438 CALL timestop(handle)
440 END SUBROUTINE matrix_almo_to_qs
456 INTEGER :: mat_distr_aos
457 REAL(kind=
dp) :: eps_filter
459 CHARACTER(len=*),
PARAMETER :: routinen =
'init_almo_ks_matrix_via_qs'
461 INTEGER :: handle, ispin, nspin
462 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_qs_ks, matrix_qs_s
468 CALL timeset(routinen, handle)
474 dft_control=dft_control, &
475 matrix_s=matrix_qs_s, &
476 matrix_ks=matrix_qs_ks, &
480 nspin = dft_control%nspins
483 IF (.NOT.
ASSOCIATED(matrix_qs_ks))
THEN
486 ALLOCATE (matrix_qs_ks(ispin)%matrix)
488 template=matrix_qs_s(1)%matrix)
490 CALL dbcsr_set(matrix_qs_ks(ispin)%matrix, 0.0_dp)
492 CALL set_ks_env(ks_env, matrix_ks=matrix_qs_ks)
497 CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
501 CALL timestop(handle)
518 CHARACTER(len=*),
PARAMETER :: routinen =
'construct_qs_mos'
520 INTEGER :: handle, ispin, ncol_fm, nrow_fm
527 CALL timeset(routinen, handle)
530 NULLIFY (mos, fm_struct_tmp, scf_env)
536 CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
538 CALL dbcsr_get_info(almo_scf_env%matrix_t(1), nfullrows_total=nrow_fm, nfullcols_total=ncol_fm)
541 DO ispin = 1, almo_scf_env%nspins
542 CALL dbcsr_get_info(almo_scf_env%matrix_t(ispin), nfullrows_total=nrow_fm, nfullcols_total=ncol_fm)
549 IF (almo_scf_env%nspins == 1)
THEN
553 nelectron=almo_scf_env%nelectrons_total, &
554 n_el_f=real(almo_scf_env%nelectrons_total,
dp), &
556 flexible_electron_count=dft_control%relax_multiplicity)
557 ELSE IF (almo_scf_env%nspins == 2)
THEN
561 nelectron=sum(almo_scf_env%nocc_of_domain(:, ispin)), &
562 n_el_f=real(sum(almo_scf_env%nocc_of_domain(:, ispin)),
dp), &
564 flexible_electron_count=dft_control%relax_multiplicity)
568 context=almo_scf_env%blacs_env, &
569 para_env=almo_scf_env%para_env)
571 CALL cp_fm_create(mo_fm_copy, fm_struct_tmp, name=
"t_orthogonal_converted_to_fm")
575 CALL init_mo_set(mos(ispin), fm_ref=mo_fm_copy, name=
'fm_mo')
581 CALL timestop(handle)
597 INTEGER,
INTENT(IN) :: mat_distr_aos
599 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_dm_to_qs_env'
601 INTEGER :: handle, ispin, nspins
605 CALL timeset(routinen, handle)
607 NULLIFY (rho, rho_ao)
608 nspins =
SIZE(matrix_p)
614 CALL matrix_almo_to_qs(matrix_p(ispin), &
615 rho_ao(ispin)%matrix, &
621 CALL timestop(handle)
639 SUBROUTINE almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, smear, kTS_sum)
642 REAL(kind=
dp) :: energy_total
643 INTEGER,
INTENT(IN) :: mat_distr_aos
644 LOGICAL,
INTENT(IN),
OPTIONAL :: smear
645 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: kts_sum
647 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_dm_to_qs_ks'
651 REAL(kind=
dp) :: entropic_term
654 CALL timeset(routinen, handle)
656 IF (
PRESENT(smear))
THEN
662 IF (
PRESENT(kts_sum))
THEN
663 entropic_term = kts_sum
665 entropic_term = 0.0_dp
677 energy%total = energy%total - energy%kTS + entropic_term
680 energy_total = energy%total
682 CALL timestop(handle)
684 END SUBROUTINE almo_dm_to_qs_ks
703 mat_distr_aos, smear, kTS_sum)
706 TYPE(
dbcsr_type),
DIMENSION(:) :: matrix_p, matrix_ks
707 REAL(kind=
dp) :: energy_total, eps_filter
708 INTEGER,
INTENT(IN) :: mat_distr_aos
709 LOGICAL,
INTENT(IN),
OPTIONAL :: smear
710 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: kts_sum
712 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_dm_to_almo_ks'
714 INTEGER :: handle, ispin, nspins
716 REAL(kind=
dp) :: entropic_term
717 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_qs_ks
719 CALL timeset(routinen, handle)
721 IF (
PRESENT(smear))
THEN
727 IF (
PRESENT(kts_sum))
THEN
728 entropic_term = kts_sum
730 entropic_term = 0.0_dp
734 CALL almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, &
736 kts_sum=entropic_term)
738 nspins =
SIZE(matrix_ks)
741 CALL get_qs_env(qs_env, matrix_ks=matrix_qs_ks)
743 CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
747 CALL timestop(handle)
763 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: energy, energy_singles_corr
769 IF (
PRESENT(energy_singles_corr))
THEN
770 qs_energy%singles_corr = energy_singles_corr
775 IF (
PRESENT(energy))
THEN
796 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_construct_quencher'
799 INTEGER :: col, contact_atom_1, contact_atom_2, domain_col, domain_map_local_entries, &
800 domain_row, global_entries, global_list_length, grid1, groupid, handle, hold, iatom, &
801 iatom2, iblock_col, iblock_row, idomain, idomain2, ientry, igrid, ineig, ineighbor, &
802 inode, inode2, ipair, ispin, jatom, jatom2, jdomain2, local_list_length, &
803 max_domain_neighbors, max_neig, mynode, nblkcols_tot, nblkrows_tot, nblks, ndomains, &
804 neig_temp, nnode2, nnodes, row, unit_nr
805 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: current_number_neighbors, domain_entries_cpu, &
806 domain_map_global, domain_map_local, first_atom_of_molecule, global_list, &
807 last_atom_of_molecule, list_length_cpu, list_offset_cpu, local_list, offset_for_cpu
808 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: domain_grid, domain_neighbor_list, &
809 domain_neighbor_list_excessive
810 INTEGER,
DIMENSION(:),
POINTER :: col_blk_size, row_blk_size
811 LOGICAL :: already_listed, block_active, &
812 delayed_increment, found, &
814 REAL(kind=
dp) :: contact1_radius, contact2_radius, &
816 r0, r1, s0, s1, trial_distance_squared
817 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: new_block
818 REAL(kind=
dp),
DIMENSION(3) :: rab
819 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: p_old_block
828 DIMENSION(:),
POINTER :: nl_iterator, nl_iterator2
833 CALL timeset(routinen, handle)
837 IF (logger%para_env%is_source())
THEN
843 ndomains = almo_scf_env%ndomains
846 particle_set=particle_set, &
847 molecule_set=molecule_set, &
855 ALLOCATE (first_atom_of_molecule(almo_scf_env%nmolecules))
856 ALLOCATE (last_atom_of_molecule(almo_scf_env%nmolecules))
858 mol_to_first_atom=first_atom_of_molecule, &
859 mol_to_last_atom=last_atom_of_molecule)
864 template=almo_scf_env%matrix_s(1), &
865 matrix_type=dbcsr_type_no_symmetry)
868 IF (sym == dbcsr_type_no_symmetry)
THEN
869 CALL dbcsr_copy(matrix_s_sym, almo_scf_env%matrix_s(1))
875 ALLOCATE (almo_scf_env%quench_t(almo_scf_env%nspins))
876 ALLOCATE (almo_scf_env%domain_map(almo_scf_env%nspins))
878 DO ispin = 1, almo_scf_env%nspins
882 matrix_qs=matrix_s(1)%matrix, &
883 almo_scf_env=almo_scf_env, &
884 name_new=
"T_QUENCHER", &
886 symmetry_new=dbcsr_type_no_symmetry, &
888 init_domains=.false.)
892 CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist, &
893 nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
895 CALL group%set_handle(groupid)
899 local_list_length = 0
905 iatom=iatom2, jatom=jatom2, inode=inode2, nnode=nnode2)
906 IF (inode2 == 1)
THEN
907 local_list_length = local_list_length + nnode2
913 ALLOCATE (local_list(2*local_list_length))
915 local_list_length = 0
919 iatom=iatom2, jatom=jatom2)
920 local_list(2*local_list_length + 1) = iatom2
921 local_list(2*local_list_length + 2) = jatom2
922 local_list_length = local_list_length + 1
927 ALLOCATE (list_length_cpu(nnodes), list_offset_cpu(nnodes))
928 CALL group%allgather(2*local_list_length, list_length_cpu)
931 list_offset_cpu(1) = 0
933 list_offset_cpu(inode) = list_offset_cpu(inode - 1) + &
934 list_length_cpu(inode - 1)
936 global_list_length = list_offset_cpu(nnodes) + list_length_cpu(nnodes)
939 ALLOCATE (global_list(global_list_length))
940 CALL group%allgatherv(local_list, global_list, &
941 list_length_cpu, list_offset_cpu)
942 DEALLOCATE (list_length_cpu, list_offset_cpu)
943 DEALLOCATE (local_list)
946 ALLOCATE (current_number_neighbors(almo_scf_env%ndomains))
947 current_number_neighbors(:) = 0
948 global_list_length = global_list_length/2
949 DO ipair = 1, global_list_length
950 iatom2 = global_list(2*(ipair - 1) + 1)
951 jatom2 = global_list(2*(ipair - 1) + 2)
952 idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
953 jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
955 current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
957 IF (idomain2 /= jdomain2)
THEN
958 current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
961 max_domain_neighbors = maxval(current_number_neighbors)
964 ALLOCATE (domain_neighbor_list_excessive(ndomains, max_domain_neighbors))
965 current_number_neighbors(:) = 1
966 DO ipair = 1, ndomains
967 domain_neighbor_list_excessive(ipair, 1) = ipair
969 DO ipair = 1, global_list_length
970 iatom2 = global_list(2*(ipair - 1) + 1)
971 jatom2 = global_list(2*(ipair - 1) + 2)
972 idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
973 jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
974 already_listed = .false.
975 DO ineighbor = 1, current_number_neighbors(idomain2)
976 IF (domain_neighbor_list_excessive(idomain2, ineighbor) == jdomain2)
THEN
977 already_listed = .true.
981 IF (.NOT. already_listed)
THEN
983 current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
984 domain_neighbor_list_excessive(idomain2, current_number_neighbors(idomain2)) = jdomain2
986 IF (idomain2 /= jdomain2)
THEN
987 current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
988 domain_neighbor_list_excessive(jdomain2, current_number_neighbors(jdomain2)) = idomain2
992 DEALLOCATE (global_list)
994 max_domain_neighbors = maxval(current_number_neighbors)
995 ALLOCATE (domain_neighbor_list(ndomains, max_domain_neighbors))
996 domain_neighbor_list(:, :) = 0
997 domain_neighbor_list(:, :) = domain_neighbor_list_excessive(:, 1:max_domain_neighbors)
998 DEALLOCATE (domain_neighbor_list_excessive)
1000 ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1001 ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(max_domain_neighbors*ndomains, 2))
1002 almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1003 almo_scf_env%domain_map(ispin)%index1(:) = 0
1004 domain_map_local_entries = 0
1010 row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1012 DO row = 1, nblkrows_tot
1013 DO col = 1, current_number_neighbors(row)
1016 iblock_col = domain_neighbor_list(row, col)
1018 iblock_row, iblock_col, hold)
1020 IF (hold == mynode)
THEN
1024 domain_row = almo_scf_env%domain_index_of_ao_block(iblock_row)
1026 domain_col = almo_scf_env%domain_index_of_mo_block(iblock_col)
1028 SELECT CASE (almo_scf_env%constraint_type)
1031 block_active = .false.
1039 IF (domain_row == domain_col)
THEN
1040 block_active = .true.
1046 cpabort(
"Illegal: atomic domains and molecular groups")
1056 cpabort(
"Illegal: molecular domains and atomic groups")
1061 IF (domain_row == domain_col)
THEN
1062 block_active = .true.
1069 IF (block_active)
THEN
1071 ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1072 new_block(:, :) = 1.0_dp
1073 CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1074 DEALLOCATE (new_block)
1076 IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains)
THEN
1077 cpabort(
"weird... max_domain_neighbors is exceeded")
1079 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1080 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1081 domain_map_local_entries = domain_map_local_entries + 1
1096 CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1098 overlap = maxval(abs(p_old_block))
1107 cpabort(
"atomic domains and molecular groups - NYI")
1118 cpabort(
"molecular domains and atomic groups - NYI")
1124 CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1126 overlap = maxval(abs(p_old_block))
1135 s0 = -log10(abs(almo_scf_env%quencher_s0))
1136 s1 = -log10(abs(almo_scf_env%quencher_s1))
1137 IF (overlap == 0.0_dp)
THEN
1138 overlap = -log10(abs(almo_scf_env%eps_filter)) + 100.0_dp
1140 overlap = -log10(overlap)
1142 IF (s0 < 0.0_dp)
THEN
1143 cpabort(
"S0 is less than zero")
1145 IF (s1 <= 0.0_dp)
THEN
1146 cpabort(
"S1 is less than or equal to zero")
1149 cpabort(
"S0 is greater than or equal to S1")
1153 IF (overlap < s1)
THEN
1154 ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1155 IF (overlap <= s0)
THEN
1156 new_block(:, :) = 1.0_dp
1158 new_block(:, :) = 1.0_dp/(1.0_dp + exp(-(s0 - s1)/(s0 - overlap) - (s0 - s1)/(overlap - s1)))
1161 IF (abs(new_block(1, 1)) > abs(almo_scf_env%eps_filter))
THEN
1162 IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains)
THEN
1163 cpabort(
"weird... max_domain_neighbors is exceeded")
1165 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1166 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1167 domain_map_local_entries = domain_map_local_entries + 1
1170 CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1171 DEALLOCATE (new_block)
1188 IF (domain_row == domain_col)
THEN
1190 contact_atom_1 = first_atom_of_molecule(domain_row)
1191 contact_atom_2 = first_atom_of_molecule(domain_col)
1196 DO iatom = first_atom_of_molecule(domain_row), last_atom_of_molecule(domain_row)
1197 DO jatom = first_atom_of_molecule(domain_col), last_atom_of_molecule(domain_col)
1198 rab(:) =
pbc(particle_set(iatom)%r(:), particle_set(jatom)%r(:), cell)
1199 trial_distance_squared = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1202 contact_atom_1 = iatom
1203 contact_atom_2 = jatom
1207 cpassert(contact_atom_1 > 0)
1215 cpabort(
"atomic domains and molecular groups - NYI")
1226 cpabort(
"molecular domains and atomic groups - NYI")
1232 rab(:) =
pbc(particle_set(domain_row)%r(:), particle_set(domain_col)%r(:), cell)
1233 distance = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1234 contact_atom_1 = domain_row
1235 contact_atom_2 = domain_col
1243 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1244 rcov=contact1_radius)
1245 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1246 rcov=contact2_radius)
1248 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1249 rvdw=contact1_radius)
1250 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1251 rvdw=contact2_radius)
1253 cpabort(
"Illegal quencher_radius_type")
1263 r0 = almo_scf_env%quencher_r0_factor*(contact1_radius + contact2_radius)
1265 r1 = almo_scf_env%quencher_r1_factor*(contact1_radius + contact2_radius)
1268 IF (r0 < 0.0_dp)
THEN
1269 cpabort(
"R0 is less than zero")
1271 IF (r1 <= 0.0_dp)
THEN
1272 cpabort(
"R1 is less than or equal to zero")
1275 cpabort(
"R0 is greater than or equal to R1")
1279 IF (distance < r1)
THEN
1280 ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1281 IF (distance <= r0)
THEN
1282 new_block(:, :) = 1.0_dp
1285 cpabort(
"distance > r0 not yet validated")
1286 new_block(:, :) = 1.0_dp/(1.0_dp + exp((r1 - r0)/(r0 - distance) + (r1 - r0)/(r1 - distance)))
1289 IF (abs(new_block(1, 1)) > abs(almo_scf_env%eps_filter))
THEN
1290 IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains)
THEN
1291 cpabort(
"weird... max_domain_neighbors is exceeded")
1293 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1294 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1295 domain_map_local_entries = domain_map_local_entries + 1
1298 CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1299 DEALLOCATE (new_block)
1303 cpabort(
"Illegal constraint type")
1311 DEALLOCATE (domain_neighbor_list)
1312 DEALLOCATE (current_number_neighbors)
1317 almo_scf_env%eps_filter)
1321 IF (nblks /= domain_map_local_entries)
THEN
1322 cpabort(
"number of blocks is wrong")
1326 ALLOCATE (domain_entries_cpu(nnodes), offset_for_cpu(nnodes))
1327 CALL group%allgather(2*domain_map_local_entries, domain_entries_cpu)
1330 offset_for_cpu(1) = 0
1331 DO inode = 2, nnodes
1332 offset_for_cpu(inode) = offset_for_cpu(inode - 1) + &
1333 domain_entries_cpu(inode - 1)
1335 global_entries = offset_for_cpu(nnodes) + domain_entries_cpu(nnodes)
1338 ALLOCATE (domain_map_global(global_entries))
1339 ALLOCATE (domain_map_local(2*domain_map_local_entries))
1340 DO ientry = 1, domain_map_local_entries
1341 domain_map_local(2*(ientry - 1) + 1) = almo_scf_env%domain_map(ispin)%pairs(ientry, 1)
1342 domain_map_local(2*ientry) = almo_scf_env%domain_map(ispin)%pairs(ientry, 2)
1344 CALL group%allgatherv(domain_map_local, domain_map_global, &
1345 domain_entries_cpu, offset_for_cpu)
1346 DEALLOCATE (domain_entries_cpu, offset_for_cpu)
1347 DEALLOCATE (domain_map_local)
1349 DEALLOCATE (almo_scf_env%domain_map(ispin)%index1)
1350 DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs)
1351 ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1352 ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(global_entries/2, 2))
1353 almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1354 almo_scf_env%domain_map(ispin)%index1(:) = 0
1360 max_neig = max_domain_neighbors
1361 max_neig_fails = .true.
1362 max_neig_loop:
DO WHILE (max_neig_fails)
1363 ALLOCATE (domain_grid(almo_scf_env%ndomains, 0:max_neig))
1364 domain_grid(:, :) = 0
1366 domain_grid(:, 0) = 1
1368 global_entries = global_entries/2
1369 DO ientry = 1, global_entries
1371 grid1 = domain_map_global(2*ientry)
1373 ineig = domain_map_global(2*(ientry - 1) + 1)
1375 IF (domain_grid(grid1, 0) > max_neig)
THEN
1378 DEALLOCATE (domain_grid)
1379 max_neig = max_neig*2
1384 delayed_increment = .false.
1385 DO igrid = 1, domain_grid(grid1, 0)
1387 IF (ineig < domain_grid(grid1, igrid))
THEN
1391 ineig = domain_grid(grid1, igrid)
1392 domain_grid(grid1, igrid) = neig_temp
1394 IF (domain_grid(grid1, igrid) == 0)
THEN
1396 domain_grid(grid1, igrid) = ineig
1398 delayed_increment = .true.
1402 IF (delayed_increment)
THEN
1403 domain_grid(grid1, 0) = domain_grid(grid1, 0) + 1
1406 cpabort(
"all records must be inserted")
1409 max_neig_fails = .false.
1410 END DO max_neig_loop
1411 DEALLOCATE (domain_map_global)
1414 DO idomain = 1, almo_scf_env%ndomains
1415 DO ineig = 1, domain_grid(idomain, 0) - 1
1416 almo_scf_env%domain_map(ispin)%pairs(ientry, 1) = domain_grid(idomain, ineig)
1417 almo_scf_env%domain_map(ispin)%pairs(ientry, 2) = idomain
1420 almo_scf_env%domain_map(ispin)%index1(idomain) = ientry
1422 DEALLOCATE (domain_grid)
1425 IF (almo_scf_env%nspins == 2)
THEN
1427 almo_scf_env%quench_t(1))
1428 almo_scf_env%domain_map(2)%pairs(:, :) = &
1429 almo_scf_env%domain_map(1)%pairs(:, :)
1430 almo_scf_env%domain_map(2)%index1(:) = &
1431 almo_scf_env%domain_map(1)%index1(:)
1438 DEALLOCATE (first_atom_of_molecule)
1439 DEALLOCATE (last_atom_of_molecule)
1442 CALL timestop(handle)
1459 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_w_matrix_almo'
1461 INTEGER :: handle, ispin
1462 REAL(kind=
dp) :: scaling
1463 TYPE(
dbcsr_type) :: tmp_nn1, tmp_no1, tmp_oo1, tmp_oo2
1465 CALL timeset(routinen, handle)
1467 IF (almo_scf_env%nspins == 1)
THEN
1473 DO ispin = 1, almo_scf_env%nspins
1475 CALL dbcsr_create(tmp_nn1, template=almo_scf_env%matrix_s(1), &
1476 matrix_type=dbcsr_type_no_symmetry)
1477 CALL dbcsr_create(tmp_no1, template=almo_scf_env%matrix_t(ispin), &
1478 matrix_type=dbcsr_type_no_symmetry)
1479 CALL dbcsr_create(tmp_oo1, template=almo_scf_env%matrix_sigma_inv(ispin), &
1480 matrix_type=dbcsr_type_no_symmetry)
1481 CALL dbcsr_create(tmp_oo2, template=almo_scf_env%matrix_sigma_inv(ispin), &
1482 matrix_type=dbcsr_type_no_symmetry)
1484 CALL dbcsr_copy(tmp_nn1, almo_scf_env%matrix_ks(ispin))
1486 CALL dbcsr_multiply(
"N",
"N", scaling, tmp_nn1, almo_scf_env%matrix_t(ispin), &
1487 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1489 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_no1, &
1490 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1492 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, tmp_oo1, almo_scf_env%matrix_sigma_inv(ispin), &
1493 0.0_dp, tmp_oo2, filter_eps=almo_scf_env%eps_filter)
1495 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_sigma_inv(ispin), tmp_oo2, &
1496 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1498 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_oo1, &
1499 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1501 CALL dbcsr_multiply(
"N",
"T", 1.0_dp, tmp_no1, almo_scf_env%matrix_t(ispin), &
1502 0.0_dp, tmp_nn1, filter_eps=almo_scf_env%eps_filter)
1503 CALL matrix_almo_to_qs(tmp_nn1, matrix_w(ispin)%matrix, almo_scf_env%mat_distr_aos)
1512 CALL timestop(handle)
double distance_squared(double *A, double *B)
Interface between ALMO SCF and QS.
subroutine, public almo_scf_update_ks_energy(qs_env, energy, energy_singles_corr)
update qs_env total energy
subroutine, public construct_qs_mos(qs_env, almo_scf_env)
Create MOs in the QS env to be able to return ALMOs to QS.
subroutine, public almo_dm_to_almo_ks(qs_env, matrix_p, matrix_ks, energy_total, eps_filter, mat_distr_aos, smear, kts_sum)
uses the ALMO density matrix to compute ALMO KS matrix and the new energy
subroutine, public calculate_w_matrix_almo(matrix_w, almo_scf_env)
Compute matrix W (energy-weighted density matrix) that is needed for the evaluation of forces.
subroutine, public matrix_almo_create(matrix_new, matrix_qs, almo_scf_env, name_new, size_keys, symmetry_new, spin_key, init_domains)
create the ALMO matrix templates
subroutine, public init_almo_ks_matrix_via_qs(qs_env, matrix_ks, mat_distr_aos, eps_filter)
Initialization of the QS and ALMO KS matrix.
subroutine, public almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
return density matrix to the qs_env
subroutine, public matrix_qs_to_almo(matrix_qs, matrix_almo, mat_distr_aos)
convert between two types of matrices: QS style to ALMO style
subroutine, public almo_scf_construct_quencher(qs_env, almo_scf_env)
Creates the matrix that imposes absolute locality on MOs.
Types for all ALMO-based methods.
integer, parameter, public almo_mat_dim_occ
integer, parameter, public almo_mat_dim_virt_full
integer, parameter, public almo_mat_dim_aobasis
integer, parameter, public almo_mat_dim_virt
integer, parameter, public almo_mat_dim_virt_disc
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_distribution_release(dist)
...
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_get_stored_coordinates(matrix, row, column, processor)
...
subroutine, public dbcsr_work_create(matrix, nblks_guess, sizedata_guess, n, work_mutable)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_finalize(matrix)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_complete_redistribute(matrix, redist)
...
integer function, public dbcsr_get_num_blocks(matrix)
...
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Routines that link DBCSR and CP2K concepts together.
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
DBCSR operations in CP2K.
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Defines the basic variable types.
integer, parameter, public dp
Interface to the message passing library MPI.
Define the data structure for the molecule information.
subroutine, public get_molecule_set_info(molecule_set, atom_to_mol, mol_to_first_atom, mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity)
returns information about molecules in the set.
Define the data structure for the particle information.
Perform a QUICKSTEP wavefunction optimization (single point)
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, mimic, 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, xcint_weights, 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, mimic, 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.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition and initialisation of the mo data type.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name, counter)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
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...
module that contains the definitions of the scf types
subroutine, public scf_env_create(scf_env)
allocates and initialize an scf_env
Type defining parameters related to the simulation cell.
keeps the information about the structure of a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.