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