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 tblite matrix build
10 : !> \author JVP
11 : !> \history creation 09.2024
12 : ! **************************************************************************************************
13 :
14 : MODULE tblite_ks_matrix
15 :
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
18 : dbcsr_copy,&
19 : dbcsr_multiply,&
20 : dbcsr_p_type,&
21 : dbcsr_type
22 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
23 : USE kinds, ONLY: dp
24 : USE message_passing, ONLY: mp_para_env_type
25 : USE qs_energy_types, ONLY: qs_energy_type
26 : USE qs_environment_types, ONLY: get_qs_env,&
27 : qs_environment_type
28 : USE qs_ks_types, ONLY: qs_ks_env_type
29 : USE qs_mo_types, ONLY: get_mo_set,&
30 : mo_set_type
31 : USE qs_rho_types, ONLY: qs_rho_get,&
32 : qs_rho_type
33 : USE tblite_interface, ONLY: tb_derive_dH_diag,&
34 : tb_derive_dH_off,&
35 : tb_ham_add_coulomb,&
36 : tb_update_charges
37 : #include "./base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tblite_ks_matrix'
44 :
45 : PUBLIC :: build_tblite_ks_matrix
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief ...
51 : !> \param qs_env ...
52 : !> \param calculate_forces ...
53 : !> \param just_energy ...
54 : !> \param ext_ks_matrix ...
55 : ! **************************************************************************************************
56 0 : SUBROUTINE build_tblite_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
57 : TYPE(qs_environment_type), POINTER :: qs_env
58 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
59 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
60 : POINTER :: ext_ks_matrix
61 :
62 : CHARACTER(len=*), PARAMETER :: routineN = 'build_tblite_ks_matrix'
63 :
64 : INTEGER :: handle, img, ispin, nimg, ns, nspins
65 : LOGICAL :: do_efield
66 : REAL(KIND=dp) :: pc_ener, qmmm_el
67 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs
68 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h
69 : TYPE(dbcsr_type), POINTER :: mo_coeff
70 : TYPE(dft_control_type), POINTER :: dft_control
71 : TYPE(mp_para_env_type), POINTER :: para_env
72 : TYPE(qs_energy_type), POINTER :: energy
73 : TYPE(qs_ks_env_type), POINTER :: ks_env
74 : TYPE(qs_rho_type), POINTER :: rho
75 :
76 0 : CALL timeset(routineN, handle)
77 :
78 0 : NULLIFY (dft_control, ks_env, ks_matrix, rho, energy)
79 0 : CPASSERT(ASSOCIATED(qs_env))
80 :
81 : CALL get_qs_env(qs_env, &
82 : dft_control=dft_control, &
83 : matrix_h_kp=matrix_h, &
84 : para_env=para_env, &
85 : ks_env=ks_env, &
86 : matrix_ks_kp=ks_matrix, &
87 : rho=rho, &
88 0 : energy=energy)
89 :
90 0 : IF (PRESENT(ext_ks_matrix)) THEN
91 : ! remap pointer to allow for non-kpoint external ks matrix
92 : ! ext_ks_matrix is used in linear response code
93 0 : ns = SIZE(ext_ks_matrix)
94 0 : ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns)
95 : END IF
96 :
97 0 : energy%qmmm_el = 0.0_dp
98 :
99 0 : nspins = dft_control%nspins
100 0 : nimg = dft_control%nimages
101 0 : CPASSERT(ASSOCIATED(matrix_h))
102 0 : CPASSERT(ASSOCIATED(rho))
103 0 : CPASSERT(SIZE(ks_matrix) > 0)
104 :
105 0 : DO ispin = 1, nspins
106 0 : DO img = 1, nimg
107 : ! copy the core matrix into the fock matrix
108 0 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix)
109 : END DO
110 : END DO
111 :
112 0 : IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
113 : dft_control%apply_efield_field) THEN
114 0 : do_efield = .TRUE.
115 0 : CPABORT("Not implemented yet. Use CP2K routines for GFN1")
116 : ELSE
117 0 : do_efield = .FALSE.
118 : END IF
119 :
120 0 : CALL tb_update_charges(qs_env, dft_control, qs_env%tb_tblite, calculate_forces, .TRUE.)
121 :
122 0 : CALL tb_ham_add_coulomb(qs_env, qs_env%tb_tblite, dft_control)
123 :
124 0 : IF (qs_env%qmmm) THEN
125 0 : CPASSERT(SIZE(ks_matrix, 2) == 1)
126 0 : DO ispin = 1, nspins
127 : ! If QM/MM sumup the 1el Hamiltonian
128 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
129 0 : 1.0_dp, 1.0_dp)
130 0 : CALL qs_rho_get(rho, rho_ao=matrix_p1)
131 : ! Compute QM/MM Energy
132 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
133 0 : matrix_p1(ispin)%matrix, qmmm_el)
134 0 : energy%qmmm_el = energy%qmmm_el + qmmm_el
135 : END DO
136 0 : pc_ener = qs_env%ks_qmmm_env%pc_ener
137 0 : energy%qmmm_el = energy%qmmm_el + pc_ener
138 : END IF
139 :
140 0 : IF (calculate_forces) THEN
141 0 : CALL tb_derive_dH_diag(qs_env, .TRUE., nimg)
142 0 : CALL tb_derive_dH_off(qs_env, .TRUE., nimg)
143 : END IF
144 :
145 : ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C
146 0 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN
147 0 : CPASSERT(SIZE(ks_matrix, 2) == 1)
148 : BLOCK
149 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
150 0 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array)
151 0 : DO ispin = 1, SIZE(mo_derivs)
152 0 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
153 0 : IF (.NOT. mo_array(ispin)%use_mo_coeff_b) THEN
154 0 : CPABORT("")
155 : END IF
156 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, &
157 0 : 0.0_dp, mo_derivs(ispin)%matrix)
158 : END DO
159 : END BLOCK
160 : END IF
161 :
162 0 : CALL timestop(handle)
163 :
164 0 : END SUBROUTINE build_tblite_ks_matrix
165 :
166 : END MODULE tblite_ks_matrix
|
