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)
350 DO row = 1, nblkrows_tot
356 IF (hold == mynode)
THEN
360 one_dim_is_mo = .false.
364 IF (one_dim_is_mo)
THEN
365 IF (almo_scf_env%nocc_of_domain(row, spin_key) == 0) active = .false.
368 one_dim_is_mo = .false.
372 IF (one_dim_is_mo)
THEN
373 IF (almo_scf_env%nvirt_of_domain(row, spin_key) == 0) active = .false.
376 one_dim_is_mo = .false.
380 IF (one_dim_is_mo)
THEN
381 IF (almo_scf_env%nvirt_disc_of_domain(row, spin_key) == 0) active = .false.
384 one_dim_is_mo = .false.
388 IF (one_dim_is_mo)
THEN
389 IF (almo_scf_env%nvirt_full_of_domain(row, spin_key) == 0) active = .false.
393 ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
394 new_block(:, :) = 1.0_dp
396 DEALLOCATE (new_block)
407 CALL timestop(handle)
423 INTEGER :: mat_distr_aos
425 CHARACTER(len=*),
PARAMETER :: routinen =
'matrix_qs_to_almo'
430 CALL timeset(routinen, handle)
433 SELECT CASE (mat_distr_aos)
439 CALL dbcsr_create(matrix_qs_nosym, template=matrix_qs, matrix_type=dbcsr_type_no_symmetry)
457 CALL timestop(handle)
470 SUBROUTINE matrix_almo_to_qs(matrix_almo, matrix_qs, mat_distr_aos)
472 INTEGER,
INTENT(IN) :: mat_distr_aos
474 CHARACTER(len=*),
PARAMETER :: routinen =
'matrix_almo_to_qs'
479 CALL timeset(routinen, handle)
482 SELECT CASE (mat_distr_aos)
484 CALL dbcsr_copy(matrix_qs, matrix_almo, keep_sparsity=.true.)
486 CALL dbcsr_create(matrix_almo_redist, template=matrix_qs)
489 CALL dbcsr_copy(matrix_qs, matrix_almo_redist, keep_sparsity=.true.)
495 CALL timestop(handle)
497 END SUBROUTINE matrix_almo_to_qs
513 INTEGER :: mat_distr_aos
514 REAL(kind=
dp) :: eps_filter
516 CHARACTER(len=*),
PARAMETER :: routinen =
'init_almo_ks_matrix_via_qs'
518 INTEGER :: handle, ispin, nspin
519 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_qs_ks, matrix_qs_s
525 CALL timeset(routinen, handle)
531 dft_control=dft_control, &
532 matrix_s=matrix_qs_s, &
533 matrix_ks=matrix_qs_ks, &
537 nspin = dft_control%nspins
540 IF (.NOT.
ASSOCIATED(matrix_qs_ks))
THEN
543 ALLOCATE (matrix_qs_ks(ispin)%matrix)
545 template=matrix_qs_s(1)%matrix)
547 CALL dbcsr_set(matrix_qs_ks(ispin)%matrix, 0.0_dp)
549 CALL set_ks_env(ks_env, matrix_ks=matrix_qs_ks)
554 CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
558 CALL timestop(handle)
575 CHARACTER(len=*),
PARAMETER :: routinen =
'construct_qs_mos'
577 INTEGER :: handle, ispin, ncol_fm, nrow_fm
584 CALL timeset(routinen, handle)
587 NULLIFY (mos, fm_struct_tmp, scf_env)
593 CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
595 CALL dbcsr_get_info(almo_scf_env%matrix_t(1), nfullrows_total=nrow_fm, nfullcols_total=ncol_fm)
598 DO ispin = 1, almo_scf_env%nspins
599 CALL dbcsr_get_info(almo_scf_env%matrix_t(ispin), nfullrows_total=nrow_fm, nfullcols_total=ncol_fm)
606 IF (almo_scf_env%nspins == 1)
THEN
610 nelectron=almo_scf_env%nelectrons_total, &
611 n_el_f=real(almo_scf_env%nelectrons_total,
dp), &
613 flexible_electron_count=dft_control%relax_multiplicity)
614 ELSEIF (almo_scf_env%nspins == 2)
THEN
618 nelectron=sum(almo_scf_env%nocc_of_domain(:, ispin)), &
619 n_el_f=real(sum(almo_scf_env%nocc_of_domain(:, ispin)),
dp), &
621 flexible_electron_count=dft_control%relax_multiplicity)
625 context=almo_scf_env%blacs_env, &
626 para_env=almo_scf_env%para_env)
628 CALL cp_fm_create(mo_fm_copy, fm_struct_tmp, name=
"t_orthogonal_converted_to_fm")
632 CALL init_mo_set(mos(ispin), fm_ref=mo_fm_copy, name=
'fm_mo')
638 CALL timestop(handle)
654 INTEGER,
INTENT(IN) :: mat_distr_aos
656 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_dm_to_qs_env'
658 INTEGER :: handle, ispin, nspins
662 CALL timeset(routinen, handle)
664 NULLIFY (rho, rho_ao)
665 nspins =
SIZE(matrix_p)
671 CALL matrix_almo_to_qs(matrix_p(ispin), &
672 rho_ao(ispin)%matrix, &
678 CALL timestop(handle)
696 SUBROUTINE almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, smear, kTS_sum)
699 REAL(kind=
dp) :: energy_total
700 INTEGER,
INTENT(IN) :: mat_distr_aos
701 LOGICAL,
INTENT(IN),
OPTIONAL :: smear
702 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: kts_sum
704 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_dm_to_qs_ks'
708 REAL(kind=
dp) :: entropic_term
711 CALL timeset(routinen, handle)
713 IF (
PRESENT(smear))
THEN
719 IF (
PRESENT(kts_sum))
THEN
720 entropic_term = kts_sum
722 entropic_term = 0.0_dp
734 energy%total = energy%total - energy%kTS + entropic_term
737 energy_total = energy%total
739 CALL timestop(handle)
741 END SUBROUTINE almo_dm_to_qs_ks
760 mat_distr_aos, smear, kTS_sum)
763 TYPE(
dbcsr_type),
DIMENSION(:) :: matrix_p, matrix_ks
764 REAL(kind=
dp) :: energy_total, eps_filter
765 INTEGER,
INTENT(IN) :: mat_distr_aos
766 LOGICAL,
INTENT(IN),
OPTIONAL :: smear
767 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: kts_sum
769 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_dm_to_almo_ks'
771 INTEGER :: handle, ispin, nspins
773 REAL(kind=
dp) :: entropic_term
774 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_qs_ks
776 CALL timeset(routinen, handle)
778 IF (
PRESENT(smear))
THEN
784 IF (
PRESENT(kts_sum))
THEN
785 entropic_term = kts_sum
787 entropic_term = 0.0_dp
791 CALL almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, &
793 kts_sum=entropic_term)
795 nspins =
SIZE(matrix_ks)
798 CALL get_qs_env(qs_env, matrix_ks=matrix_qs_ks)
800 CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
804 CALL timestop(handle)
820 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: energy, energy_singles_corr
826 IF (
PRESENT(energy_singles_corr))
THEN
827 qs_energy%singles_corr = energy_singles_corr
832 IF (
PRESENT(energy))
THEN
853 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_construct_quencher'
856 INTEGER :: col, contact_atom_1, contact_atom_2, domain_col, domain_map_local_entries, &
857 domain_row, global_entries, global_list_length, grid1, groupid, handle, hold, iatom, &
858 iatom2, iblock_col, iblock_row, idomain, idomain2, ientry, igrid, ineig, ineighbor, &
859 inode, inode2, ipair, ispin, jatom, jatom2, jdomain2, local_list_length, &
860 max_domain_neighbors, max_neig, mynode, nblkcols_tot, nblkrows_tot, nblks, ndomains, &
861 neig_temp, nnode2, nnodes, row, unit_nr
862 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: current_number_neighbors, domain_entries_cpu, &
863 domain_map_global, domain_map_local, first_atom_of_molecule, global_list, &
864 last_atom_of_molecule, list_length_cpu, list_offset_cpu, local_list, offset_for_cpu
865 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: domain_grid, domain_neighbor_list, &
866 domain_neighbor_list_excessive
867 INTEGER,
DIMENSION(:),
POINTER :: col_blk_size, row_blk_size
868 LOGICAL :: already_listed, block_active, &
869 delayed_increment, found, &
871 REAL(kind=
dp) :: contact1_radius, contact2_radius, &
873 r0, r1, s0, s1, trial_distance_squared
874 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: new_block
875 REAL(kind=
dp),
DIMENSION(3) :: rab
876 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: p_old_block
885 DIMENSION(:),
POINTER :: nl_iterator, nl_iterator2
890 CALL timeset(routinen, handle)
894 IF (logger%para_env%is_source())
THEN
900 ndomains = almo_scf_env%ndomains
903 particle_set=particle_set, &
904 molecule_set=molecule_set, &
912 ALLOCATE (first_atom_of_molecule(almo_scf_env%nmolecules))
913 ALLOCATE (last_atom_of_molecule(almo_scf_env%nmolecules))
915 mol_to_first_atom=first_atom_of_molecule, &
916 mol_to_last_atom=last_atom_of_molecule)
921 template=almo_scf_env%matrix_s(1), &
922 matrix_type=dbcsr_type_no_symmetry)
925 IF (sym == dbcsr_type_no_symmetry)
THEN
926 CALL dbcsr_copy(matrix_s_sym, almo_scf_env%matrix_s(1))
932 ALLOCATE (almo_scf_env%quench_t(almo_scf_env%nspins))
933 ALLOCATE (almo_scf_env%domain_map(almo_scf_env%nspins))
935 DO ispin = 1, almo_scf_env%nspins
939 matrix_qs=matrix_s(1)%matrix, &
940 almo_scf_env=almo_scf_env, &
941 name_new=
"T_QUENCHER", &
943 symmetry_new=dbcsr_type_no_symmetry, &
945 init_domains=.false.)
949 CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist, &
950 nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
952 CALL group%set_handle(groupid)
956 local_list_length = 0
962 iatom=iatom2, jatom=jatom2, inode=inode2, nnode=nnode2)
963 IF (inode2 == 1)
THEN
964 local_list_length = local_list_length + nnode2
970 ALLOCATE (local_list(2*local_list_length))
972 local_list_length = 0
976 iatom=iatom2, jatom=jatom2)
977 local_list(2*local_list_length + 1) = iatom2
978 local_list(2*local_list_length + 2) = jatom2
979 local_list_length = local_list_length + 1
984 ALLOCATE (list_length_cpu(nnodes), list_offset_cpu(nnodes))
985 CALL group%allgather(2*local_list_length, list_length_cpu)
988 list_offset_cpu(1) = 0
990 list_offset_cpu(inode) = list_offset_cpu(inode - 1) + &
991 list_length_cpu(inode - 1)
993 global_list_length = list_offset_cpu(nnodes) + list_length_cpu(nnodes)
996 ALLOCATE (global_list(global_list_length))
997 CALL group%allgatherv(local_list, global_list, &
998 list_length_cpu, list_offset_cpu)
999 DEALLOCATE (list_length_cpu, list_offset_cpu)
1000 DEALLOCATE (local_list)
1003 ALLOCATE (current_number_neighbors(almo_scf_env%ndomains))
1004 current_number_neighbors(:) = 0
1005 global_list_length = global_list_length/2
1006 DO ipair = 1, global_list_length
1007 iatom2 = global_list(2*(ipair - 1) + 1)
1008 jatom2 = global_list(2*(ipair - 1) + 2)
1009 idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
1010 jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
1012 current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
1014 IF (idomain2 /= jdomain2)
THEN
1015 current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
1018 max_domain_neighbors = maxval(current_number_neighbors)
1021 ALLOCATE (domain_neighbor_list_excessive(ndomains, max_domain_neighbors))
1022 current_number_neighbors(:) = 1
1023 DO ipair = 1, ndomains
1024 domain_neighbor_list_excessive(ipair, 1) = ipair
1026 DO ipair = 1, global_list_length
1027 iatom2 = global_list(2*(ipair - 1) + 1)
1028 jatom2 = global_list(2*(ipair - 1) + 2)
1029 idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
1030 jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
1031 already_listed = .false.
1032 DO ineighbor = 1, current_number_neighbors(idomain2)
1033 IF (domain_neighbor_list_excessive(idomain2, ineighbor) == jdomain2)
THEN
1034 already_listed = .true.
1038 IF (.NOT. already_listed)
THEN
1040 current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
1041 domain_neighbor_list_excessive(idomain2, current_number_neighbors(idomain2)) = jdomain2
1043 IF (idomain2 /= jdomain2)
THEN
1044 current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
1045 domain_neighbor_list_excessive(jdomain2, current_number_neighbors(jdomain2)) = idomain2
1049 DEALLOCATE (global_list)
1051 max_domain_neighbors = maxval(current_number_neighbors)
1052 ALLOCATE (domain_neighbor_list(ndomains, max_domain_neighbors))
1053 domain_neighbor_list(:, :) = 0
1054 domain_neighbor_list(:, :) = domain_neighbor_list_excessive(:, 1:max_domain_neighbors)
1055 DEALLOCATE (domain_neighbor_list_excessive)
1057 ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1058 ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(max_domain_neighbors*ndomains, 2))
1059 almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1060 almo_scf_env%domain_map(ispin)%index1(:) = 0
1061 domain_map_local_entries = 0
1067 row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1069 DO row = 1, nblkrows_tot
1070 DO col = 1, current_number_neighbors(row)
1073 iblock_col = domain_neighbor_list(row, col)
1075 iblock_row, iblock_col, hold)
1077 IF (hold == mynode)
THEN
1081 domain_row = almo_scf_env%domain_index_of_ao_block(iblock_row)
1083 domain_col = almo_scf_env%domain_index_of_mo_block(iblock_col)
1085 SELECT CASE (almo_scf_env%constraint_type)
1088 block_active = .false.
1096 IF (domain_row == domain_col)
THEN
1097 block_active = .true.
1103 cpabort(
"Illegal: atomic domains and molecular groups")
1113 cpabort(
"Illegal: molecular domains and atomic groups")
1118 IF (domain_row == domain_col)
THEN
1119 block_active = .true.
1126 IF (block_active)
THEN
1128 ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1129 new_block(:, :) = 1.0_dp
1130 CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1131 DEALLOCATE (new_block)
1133 IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains)
THEN
1134 cpabort(
"weird... max_domain_neighbors is exceeded")
1136 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1137 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1138 domain_map_local_entries = domain_map_local_entries + 1
1153 CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1155 overlap = maxval(abs(p_old_block))
1164 cpabort(
"atomic domains and molecular groups - NYI")
1175 cpabort(
"molecular domains and atomic groups - NYI")
1181 CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1183 overlap = maxval(abs(p_old_block))
1192 s0 = -log10(abs(almo_scf_env%quencher_s0))
1193 s1 = -log10(abs(almo_scf_env%quencher_s1))
1194 IF (overlap == 0.0_dp)
THEN
1195 overlap = -log10(abs(almo_scf_env%eps_filter)) + 100.0_dp
1197 overlap = -log10(overlap)
1199 IF (s0 < 0.0_dp)
THEN
1200 cpabort(
"S0 is less than zero")
1202 IF (s1 <= 0.0_dp)
THEN
1203 cpabort(
"S1 is less than or equal to zero")
1206 cpabort(
"S0 is greater than or equal to S1")
1210 IF (overlap < s1)
THEN
1211 ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1212 IF (overlap <= s0)
THEN
1213 new_block(:, :) = 1.0_dp
1215 new_block(:, :) = 1.0_dp/(1.0_dp + exp(-(s0 - s1)/(s0 - overlap) - (s0 - s1)/(overlap - s1)))
1218 IF (abs(new_block(1, 1)) > abs(almo_scf_env%eps_filter))
THEN
1219 IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains)
THEN
1220 cpabort(
"weird... max_domain_neighbors is exceeded")
1222 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1223 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1224 domain_map_local_entries = domain_map_local_entries + 1
1227 CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1228 DEALLOCATE (new_block)
1245 IF (domain_row == domain_col)
THEN
1247 contact_atom_1 = first_atom_of_molecule(domain_row)
1248 contact_atom_2 = first_atom_of_molecule(domain_col)
1253 DO iatom = first_atom_of_molecule(domain_row), last_atom_of_molecule(domain_row)
1254 DO jatom = first_atom_of_molecule(domain_col), last_atom_of_molecule(domain_col)
1255 rab(:) =
pbc(particle_set(iatom)%r(:), particle_set(jatom)%r(:), cell)
1256 trial_distance_squared = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1259 contact_atom_1 = iatom
1260 contact_atom_2 = jatom
1264 cpassert(contact_atom_1 > 0)
1272 cpabort(
"atomic domains and molecular groups - NYI")
1283 cpabort(
"molecular domains and atomic groups - NYI")
1289 rab(:) =
pbc(particle_set(domain_row)%r(:), particle_set(domain_col)%r(:), cell)
1290 distance = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1291 contact_atom_1 = domain_row
1292 contact_atom_2 = domain_col
1300 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1301 rcov=contact1_radius)
1302 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1303 rcov=contact2_radius)
1305 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1306 rvdw=contact1_radius)
1307 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1308 rvdw=contact2_radius)
1310 cpabort(
"Illegal quencher_radius_type")
1320 r0 = almo_scf_env%quencher_r0_factor*(contact1_radius + contact2_radius)
1322 r1 = almo_scf_env%quencher_r1_factor*(contact1_radius + contact2_radius)
1325 IF (r0 < 0.0_dp)
THEN
1326 cpabort(
"R0 is less than zero")
1328 IF (r1 <= 0.0_dp)
THEN
1329 cpabort(
"R1 is less than or equal to zero")
1332 cpabort(
"R0 is greater than or equal to R1")
1336 IF (distance < r1)
THEN
1337 ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1338 IF (distance <= r0)
THEN
1339 new_block(:, :) = 1.0_dp
1343 new_block(:, :) = 1.0_dp/(1.0_dp + exp((r1 - r0)/(r0 - distance) + (r1 - r0)/(r1 - distance)))
1346 IF (abs(new_block(1, 1)) > abs(almo_scf_env%eps_filter))
THEN
1347 IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains)
THEN
1348 cpabort(
"weird... max_domain_neighbors is exceeded")
1350 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1351 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1352 domain_map_local_entries = domain_map_local_entries + 1
1355 CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1356 DEALLOCATE (new_block)
1360 cpabort(
"Illegal constraint type")
1368 DEALLOCATE (domain_neighbor_list)
1369 DEALLOCATE (current_number_neighbors)
1374 almo_scf_env%eps_filter)
1378 IF (nblks /= domain_map_local_entries)
THEN
1379 cpabort(
"number of blocks is wrong")
1383 ALLOCATE (domain_entries_cpu(nnodes), offset_for_cpu(nnodes))
1384 CALL group%allgather(2*domain_map_local_entries, domain_entries_cpu)
1387 offset_for_cpu(1) = 0
1388 DO inode = 2, nnodes
1389 offset_for_cpu(inode) = offset_for_cpu(inode - 1) + &
1390 domain_entries_cpu(inode - 1)
1392 global_entries = offset_for_cpu(nnodes) + domain_entries_cpu(nnodes)
1395 ALLOCATE (domain_map_global(global_entries))
1396 ALLOCATE (domain_map_local(2*domain_map_local_entries))
1397 DO ientry = 1, domain_map_local_entries
1398 domain_map_local(2*(ientry - 1) + 1) = almo_scf_env%domain_map(ispin)%pairs(ientry, 1)
1399 domain_map_local(2*ientry) = almo_scf_env%domain_map(ispin)%pairs(ientry, 2)
1401 CALL group%allgatherv(domain_map_local, domain_map_global, &
1402 domain_entries_cpu, offset_for_cpu)
1403 DEALLOCATE (domain_entries_cpu, offset_for_cpu)
1404 DEALLOCATE (domain_map_local)
1406 DEALLOCATE (almo_scf_env%domain_map(ispin)%index1)
1407 DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs)
1408 ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1409 ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(global_entries/2, 2))
1410 almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1411 almo_scf_env%domain_map(ispin)%index1(:) = 0
1417 max_neig = max_domain_neighbors
1418 max_neig_fails = .true.
1419 max_neig_loop:
DO WHILE (max_neig_fails)
1420 ALLOCATE (domain_grid(almo_scf_env%ndomains, 0:max_neig))
1421 domain_grid(:, :) = 0
1423 domain_grid(:, 0) = 1
1425 global_entries = global_entries/2
1426 DO ientry = 1, global_entries
1428 grid1 = domain_map_global(2*ientry)
1430 ineig = domain_map_global(2*(ientry - 1) + 1)
1432 IF (domain_grid(grid1, 0) > max_neig)
THEN
1435 DEALLOCATE (domain_grid)
1436 max_neig = max_neig*2
1441 delayed_increment = .false.
1442 DO igrid = 1, domain_grid(grid1, 0)
1444 IF (ineig < domain_grid(grid1, igrid))
THEN
1448 ineig = domain_grid(grid1, igrid)
1449 domain_grid(grid1, igrid) = neig_temp
1451 IF (domain_grid(grid1, igrid) == 0)
THEN
1453 domain_grid(grid1, igrid) = ineig
1455 delayed_increment = .true.
1459 IF (delayed_increment)
THEN
1460 domain_grid(grid1, 0) = domain_grid(grid1, 0) + 1
1463 cpabort(
"all records must be inserted")
1466 max_neig_fails = .false.
1467 END DO max_neig_loop
1468 DEALLOCATE (domain_map_global)
1471 DO idomain = 1, almo_scf_env%ndomains
1472 DO ineig = 1, domain_grid(idomain, 0) - 1
1473 almo_scf_env%domain_map(ispin)%pairs(ientry, 1) = domain_grid(idomain, ineig)
1474 almo_scf_env%domain_map(ispin)%pairs(ientry, 2) = idomain
1477 almo_scf_env%domain_map(ispin)%index1(idomain) = ientry
1479 DEALLOCATE (domain_grid)
1482 IF (almo_scf_env%nspins == 2)
THEN
1484 almo_scf_env%quench_t(1))
1485 almo_scf_env%domain_map(2)%pairs(:, :) = &
1486 almo_scf_env%domain_map(1)%pairs(:, :)
1487 almo_scf_env%domain_map(2)%index1(:) = &
1488 almo_scf_env%domain_map(1)%index1(:)
1495 DEALLOCATE (first_atom_of_molecule)
1496 DEALLOCATE (last_atom_of_molecule)
1516 CALL timestop(handle)
1533 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_w_matrix_almo'
1535 INTEGER :: handle, ispin
1536 REAL(kind=
dp) :: scaling
1537 TYPE(
dbcsr_type) :: tmp_nn1, tmp_no1, tmp_oo1, tmp_oo2
1539 CALL timeset(routinen, handle)
1541 IF (almo_scf_env%nspins == 1)
THEN
1547 DO ispin = 1, almo_scf_env%nspins
1549 CALL dbcsr_create(tmp_nn1, template=almo_scf_env%matrix_s(1), &
1550 matrix_type=dbcsr_type_no_symmetry)
1551 CALL dbcsr_create(tmp_no1, template=almo_scf_env%matrix_t(ispin), &
1552 matrix_type=dbcsr_type_no_symmetry)
1553 CALL dbcsr_create(tmp_oo1, template=almo_scf_env%matrix_sigma_inv(ispin), &
1554 matrix_type=dbcsr_type_no_symmetry)
1555 CALL dbcsr_create(tmp_oo2, template=almo_scf_env%matrix_sigma_inv(ispin), &
1556 matrix_type=dbcsr_type_no_symmetry)
1558 CALL dbcsr_copy(tmp_nn1, almo_scf_env%matrix_ks(ispin))
1560 CALL dbcsr_multiply(
"N",
"N", scaling, tmp_nn1, almo_scf_env%matrix_t(ispin), &
1561 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1563 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_no1, &
1564 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1566 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, tmp_oo1, almo_scf_env%matrix_sigma_inv(ispin), &
1567 0.0_dp, tmp_oo2, filter_eps=almo_scf_env%eps_filter)
1569 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_sigma_inv(ispin), tmp_oo2, &
1570 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1572 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_oo1, &
1573 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1575 CALL dbcsr_multiply(
"N",
"T", 1.0_dp, tmp_no1, almo_scf_env%matrix_t(ispin), &
1576 0.0_dp, tmp_nn1, filter_eps=almo_scf_env%eps_filter)
1577 CALL matrix_almo_to_qs(tmp_nn1, matrix_w(ispin)%matrix, almo_scf_env%mat_distr_aos)
1586 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.