Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Generate an initial guess (dm and orb) from EHT calculation
10 : ! **************************************************************************************************
11 : MODULE qs_eht_guess
12 : USE basis_set_types, ONLY: gto_basis_set_p_type
13 : USE cp_blacs_env, ONLY: cp_blacs_env_type
14 : USE cp_control_types, ONLY: dft_control_type
15 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
16 : dbcsr_desymmetrize,&
17 : dbcsr_get_info,&
18 : dbcsr_p_type,&
19 : dbcsr_release,&
20 : dbcsr_type,&
21 : dbcsr_type_no_symmetry
22 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
23 : cp_dbcsr_sm_fm_multiply,&
24 : dbcsr_deallocate_matrix_set
25 : USE cp_fm_basic_linalg, ONLY: cp_fm_invert
26 : USE cp_fm_diag, ONLY: cp_fm_geeig
27 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
28 : cp_fm_struct_release,&
29 : cp_fm_struct_type
30 : USE cp_fm_types, ONLY: cp_fm_create,&
31 : cp_fm_release,&
32 : cp_fm_to_fm,&
33 : cp_fm_type
34 : USE cp_subsys_types, ONLY: cp_subsys_type
35 : USE input_constants, ONLY: do_method_xtb
36 : USE input_section_types, ONLY: section_vals_duplicate,&
37 : section_vals_get_subs_vals,&
38 : section_vals_release,&
39 : section_vals_type,&
40 : section_vals_val_set
41 : USE kinds, ONLY: dp
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE parallel_gemm_api, ONLY: parallel_gemm
44 : USE qs_energy_init, ONLY: qs_energies_init
45 : USE qs_environment, ONLY: qs_init
46 : USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env
47 : USE qs_environment_types, ONLY: get_qs_env,&
48 : qs_env_create,&
49 : qs_env_release,&
50 : qs_environment_type
51 : USE qs_integral_utils, ONLY: basis_set_list_setup
52 : USE qs_kind_types, ONLY: qs_kind_type
53 : USE qs_ks_types, ONLY: qs_ks_env_type
54 : USE qs_mo_occupation, ONLY: set_mo_occupation
55 : USE qs_mo_types, ONLY: get_mo_set,&
56 : mo_set_type
57 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
58 : USE qs_overlap, ONLY: build_overlap_matrix_simple
59 : USE tblite_ks_matrix, ONLY: build_tblite_ks_matrix
60 : USE xtb_ks_matrix, ONLY: build_xtb_ks_matrix
61 : #include "./base/base_uses.f90"
62 :
63 : IMPLICIT NONE
64 :
65 : PRIVATE
66 :
67 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_eht_guess'
68 :
69 : PUBLIC :: calculate_eht_guess
70 :
71 : ! **************************************************************************************************
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief EHT MO guess calclulation
77 : !> \param qs_env ...
78 : !> \param mo_array ...
79 : ! **************************************************************************************************
80 4 : SUBROUTINE calculate_eht_guess(qs_env, mo_array)
81 : TYPE(qs_environment_type), POINTER :: qs_env
82 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
83 :
84 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_eht_guess'
85 :
86 : INTEGER :: handle, ispin, nao, nbas, neeht, neorb, &
87 : nkind, nmo, nspins, zero
88 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
89 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigval
90 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
91 : TYPE(cp_fm_struct_type), POINTER :: mstruct_ee, mstruct_oe, mstruct_oo
92 : TYPE(cp_fm_type) :: fmksmat, fmorb, fmscr, fmsmat, fmvec, &
93 : fmwork, sfull, sinv
94 : TYPE(cp_fm_type), POINTER :: mo_coeff
95 : TYPE(cp_subsys_type), POINTER :: cp_subsys
96 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_s, matrix_t, smat
97 : TYPE(dbcsr_type) :: tempmat, tmat
98 : TYPE(dft_control_type), POINTER :: dft_control
99 4 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
100 : TYPE(mp_para_env_type), POINTER :: para_env
101 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
102 4 : POINTER :: sab_nl
103 : TYPE(qs_environment_type), POINTER :: eht_env
104 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
105 : TYPE(qs_ks_env_type), POINTER :: ks_env
106 : TYPE(section_vals_type), POINTER :: dft_section, eht_force_env_section, &
107 : force_env_section, qs_section, &
108 : subsys_section, xtb_section
109 :
110 4 : CALL timeset(routineN, handle)
111 :
112 4 : NULLIFY (subsys_section)
113 : CALL get_qs_env(qs_env, &
114 : ks_env=ks_env, &
115 : para_env=para_env, &
116 : input=force_env_section, &
117 : cp_subsys=cp_subsys, &
118 4 : dft_control=dft_control)
119 :
120 4 : NULLIFY (eht_force_env_section)
121 4 : CALL section_vals_duplicate(force_env_section, eht_force_env_section)
122 4 : dft_section => section_vals_get_subs_vals(eht_force_env_section, "DFT")
123 4 : qs_section => section_vals_get_subs_vals(dft_section, "QS")
124 4 : CALL section_vals_val_set(qs_section, "METHOD", i_val=do_method_xtb)
125 4 : xtb_section => section_vals_get_subs_vals(qs_section, "xTB")
126 4 : zero = 0
127 4 : CALL section_vals_val_set(xtb_section, "GFN_TYPE", i_val=zero)
128 : !
129 4 : ALLOCATE (eht_env)
130 4 : CALL qs_env_create(eht_env)
131 : CALL qs_init(eht_env, para_env, cp_subsys=cp_subsys, &
132 : force_env_section=eht_force_env_section, &
133 : subsys_section=subsys_section, &
134 4 : use_motion_section=.FALSE., silent=.TRUE.)
135 : !
136 4 : CALL get_qs_env(qs_env, nelectron_total=neorb)
137 4 : CALL get_qs_env(eht_env, nelectron_total=neeht)
138 4 : IF (neorb /= neeht) THEN
139 0 : CPWARN("EHT has different number of electrons than calculation method.")
140 0 : CPABORT("EHT Initial Guess")
141 : END IF
142 : !
143 4 : CALL qs_env_rebuild_pw_env(eht_env)
144 4 : CALL qs_energies_init(eht_env, calc_forces=.FALSE.)
145 4 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
146 0 : CALL build_tblite_ks_matrix(eht_env, .FALSE., .FALSE.)
147 : ELSE
148 4 : CALL build_xtb_ks_matrix(eht_env, .FALSE., .FALSE.)
149 : END IF
150 : !
151 : CALL get_qs_env(eht_env, &
152 4 : matrix_s=smat, matrix_ks=ksmat)
153 4 : nspins = SIZE(ksmat, 1)
154 4 : CALL get_qs_env(eht_env, para_env=para_env, blacs_env=blacs_env)
155 4 : CALL dbcsr_get_info(smat(1)%matrix, nfullrows_total=nao)
156 : CALL cp_fm_struct_create(fmstruct=mstruct_ee, context=blacs_env, &
157 4 : nrow_global=nao, ncol_global=nao)
158 4 : CALL cp_fm_create(fmksmat, mstruct_ee)
159 4 : CALL cp_fm_create(fmsmat, mstruct_ee)
160 4 : CALL cp_fm_create(fmvec, mstruct_ee)
161 4 : CALL cp_fm_create(fmwork, mstruct_ee)
162 12 : ALLOCATE (eigenvalues(nao))
163 :
164 : ! DBCSR matrix
165 4 : CALL dbcsr_create(tempmat, template=smat(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
166 :
167 : ! transfer to FM
168 4 : CALL dbcsr_desymmetrize(smat(1)%matrix, tempmat)
169 4 : CALL copy_dbcsr_to_fm(tempmat, fmsmat)
170 :
171 : !SINV of origianl basis
172 4 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
173 4 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
174 4 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nbas)
175 4 : CALL dbcsr_create(tmat, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
176 : CALL cp_fm_struct_create(fmstruct=mstruct_oo, context=blacs_env, &
177 4 : nrow_global=nbas, ncol_global=nbas)
178 4 : CALL cp_fm_create(sfull, mstruct_oo)
179 4 : CALL cp_fm_create(sinv, mstruct_oo)
180 4 : CALL dbcsr_desymmetrize(matrix_s(1)%matrix, tmat)
181 4 : CALL copy_dbcsr_to_fm(tmat, sfull)
182 4 : CALL cp_fm_invert(sfull, sinv)
183 4 : CALL dbcsr_release(tmat)
184 4 : CALL cp_fm_release(sfull)
185 : !TMAT(bas1, bas2)
186 4 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, sab_all=sab_nl, nkind=nkind)
187 4 : IF (.NOT. ASSOCIATED(sab_nl)) THEN
188 0 : CPWARN("Full neighborlist not available for this method. EHT initial guess not possible.")
189 0 : CPABORT("EHT Initial Guess")
190 : END IF
191 36 : ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
192 4 : CALL basis_set_list_setup(basis_set_list_a, "ORB", qs_kind_set)
193 4 : CALL get_qs_env(eht_env, qs_kind_set=qs_kind_set)
194 4 : CALL basis_set_list_setup(basis_set_list_b, "ORB", qs_kind_set)
195 : !
196 4 : NULLIFY (matrix_t)
197 : CALL build_overlap_matrix_simple(ks_env, matrix_t, &
198 4 : basis_set_list_a, basis_set_list_b, sab_nl)
199 4 : DEALLOCATE (basis_set_list_a, basis_set_list_b)
200 :
201 : ! KS matrix is not spin dependent!
202 4 : CALL dbcsr_desymmetrize(ksmat(1)%matrix, tempmat)
203 4 : CALL copy_dbcsr_to_fm(tempmat, fmksmat)
204 : ! diagonalize
205 4 : CALL cp_fm_geeig(fmksmat, fmsmat, fmvec, eigenvalues, fmwork)
206 : ! Sinv*T*d
207 : CALL cp_fm_struct_create(fmstruct=mstruct_oe, context=blacs_env, &
208 4 : nrow_global=nbas, ncol_global=nao)
209 4 : CALL cp_fm_create(fmscr, mstruct_oe)
210 4 : CALL cp_fm_create(fmorb, mstruct_oe)
211 4 : CALL cp_dbcsr_sm_fm_multiply(matrix_t(1)%matrix, fmvec, fmscr, ncol=nao)
212 4 : CALL parallel_gemm('N', 'N', nbas, nao, nbas, 1.0_dp, sinv, fmscr, 0.0_dp, fmorb)
213 : !
214 8 : DO ispin = 1, nspins
215 4 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, nmo=nmo)
216 4 : CALL cp_fm_to_fm(fmorb, mo_coeff, nmo, 1, 1)
217 4 : NULLIFY (eigval)
218 4 : CALL get_mo_set(mo_set=mo_array(ispin), eigenvalues=eigval)
219 12 : IF (ASSOCIATED(eigval)) THEN
220 20 : eigval(1:nmo) = eigenvalues(1:nmo)
221 : END IF
222 : END DO
223 4 : CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear)
224 :
225 4 : DEALLOCATE (eigenvalues)
226 4 : CALL dbcsr_release(tempmat)
227 4 : CALL dbcsr_deallocate_matrix_set(matrix_t)
228 4 : CALL cp_fm_release(fmksmat)
229 4 : CALL cp_fm_release(fmsmat)
230 4 : CALL cp_fm_release(fmvec)
231 4 : CALL cp_fm_release(fmwork)
232 4 : CALL cp_fm_release(fmscr)
233 4 : CALL cp_fm_release(fmorb)
234 4 : CALL cp_fm_release(sinv)
235 4 : CALL cp_fm_struct_release(mstruct_ee)
236 4 : CALL cp_fm_struct_release(mstruct_oe)
237 4 : CALL cp_fm_struct_release(mstruct_oo)
238 : !
239 4 : CALL qs_env_release(eht_env)
240 4 : DEALLOCATE (eht_env)
241 4 : CALL section_vals_release(eht_force_env_section)
242 :
243 4 : CALL timestop(handle)
244 :
245 28 : END SUBROUTINE calculate_eht_guess
246 :
247 : END MODULE qs_eht_guess
|
