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 Calculation of the derivative of the QMMM Hamiltonian integral
10 : !> matrix <a|\sum_i q_i|b> for semi-empirical methods
11 : !> \author Teodoro Laino - 04.2007 [tlaino]
12 : ! **************************************************************************************************
13 : MODULE qmmm_se_forces
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE cell_types, ONLY: cell_type,&
17 : pbc
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
20 : dbcsr_p_type
21 : USE input_constants, ONLY: &
22 : do_method_am1, do_method_mndo, do_method_mndod, do_method_pchg, do_method_pdg, &
23 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1
24 : USE kinds, ONLY: dp
25 : USE message_passing, ONLY: mp_para_env_type
26 : USE multipole_types, ONLY: do_multipole_none
27 : USE particle_types, ONLY: particle_type
28 : USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
29 : qmmm_pot_p_type,&
30 : qmmm_pot_type
31 : USE qmmm_util, ONLY: spherical_cutoff_factor
32 : USE qs_environment_types, ONLY: get_qs_env,&
33 : qs_environment_type
34 : USE qs_kind_types, ONLY: get_qs_kind,&
35 : qs_kind_type
36 : USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
37 : USE qs_rho_types, ONLY: qs_rho_get,&
38 : qs_rho_type
39 : USE semi_empirical_int_arrays, ONLY: se_orbital_pointer
40 : USE semi_empirical_integrals, ONLY: dcorecore,&
41 : drotnuc
42 : USE semi_empirical_types, ONLY: get_se_param,&
43 : se_int_control_type,&
44 : se_taper_type,&
45 : semi_empirical_create,&
46 : semi_empirical_release,&
47 : semi_empirical_type,&
48 : setup_se_int_control_type
49 : USE semi_empirical_utils, ONLY: get_se_type,&
50 : se_param_set_default
51 : #include "./base/base_uses.f90"
52 :
53 : IMPLICIT NONE
54 :
55 : PRIVATE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_se_forces'
58 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
59 : PUBLIC :: deriv_se_qmmm_matrix
60 :
61 : CONTAINS
62 :
63 : ! **************************************************************************************************
64 : !> \brief Constructs the derivative w.r.t. 1-el semi-empirical hamiltonian
65 : !> QMMM terms
66 : !> \param qs_env ...
67 : !> \param qmmm_env ...
68 : !> \param particles_mm ...
69 : !> \param mm_cell ...
70 : !> \param para_env ...
71 : !> \param calc_force ...
72 : !> \param Forces ...
73 : !> \param Forces_added_charges ...
74 : !> \author Teodoro Laino 04.2007 [created]
75 : ! **************************************************************************************************
76 936 : SUBROUTINE deriv_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
77 : calc_force, Forces, Forces_added_charges)
78 :
79 : TYPE(qs_environment_type), POINTER :: qs_env
80 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
81 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
82 : TYPE(cell_type), POINTER :: mm_cell
83 : TYPE(mp_para_env_type), POINTER :: para_env
84 : LOGICAL, INTENT(in), OPTIONAL :: calc_force
85 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces, Forces_added_charges
86 :
87 : CHARACTER(len=*), PARAMETER :: routineN = 'deriv_se_qmmm_matrix'
88 :
89 : INTEGER :: handle, i, iatom, ikind, iqm, ispin, &
90 : itype, natom, natorb_a, nkind, &
91 : number_qm_atoms
92 936 : INTEGER, DIMENSION(:), POINTER :: list
93 : LOGICAL :: anag, defined, found
94 : REAL(KIND=dp) :: delta, enuclear
95 936 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces_QM, p_block_a
96 936 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
97 936 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
98 : TYPE(dft_control_type), POINTER :: dft_control
99 936 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
100 936 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
101 : TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
102 : TYPE(qs_rho_type), POINTER :: rho
103 : TYPE(se_int_control_type) :: se_int_control
104 : TYPE(se_taper_type), POINTER :: se_taper
105 : TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm
106 :
107 936 : CALL timeset(routineN, handle)
108 936 : IF (calc_force) THEN
109 788 : NULLIFY (rho, atomic_kind_set, qs_kind_set, se_taper)
110 788 : NULLIFY (se_kind_a, se_kind_mm, particles_qm)
111 : CALL get_qs_env(qs_env=qs_env, &
112 : rho=rho, &
113 : se_taper=se_taper, &
114 : atomic_kind_set=atomic_kind_set, &
115 : qs_kind_set=qs_kind_set, &
116 : ks_qmmm_env=ks_qmmm_env_loc, &
117 : dft_control=dft_control, &
118 : particle_set=particles_qm, &
119 788 : natom=number_qm_atoms)
120 788 : SELECT CASE (dft_control%qs_control%method_id)
121 : CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, &
122 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
123 : ! Go on with the calculation..
124 : CASE DEFAULT
125 : ! Otherwise stop..
126 788 : CPABORT("Method not available")
127 : END SELECT
128 788 : anag = dft_control%qs_control%se_control%analytical_gradients
129 788 : delta = dft_control%qs_control%se_control%delta
130 : ! Setup SE integral control type
131 : CALL setup_se_int_control_type( &
132 : se_int_control, shortrange=.FALSE., do_ewald_r3=.FALSE., &
133 : do_ewald_gks=.FALSE., integral_screening=dft_control%qs_control%se_control%integral_screening, &
134 788 : max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
135 :
136 : ! Create a fake semi-empirical type to handle the classical atom
137 2364 : ALLOCATE (Forces_QM(3, number_qm_atoms))
138 788 : CALL semi_empirical_create(se_kind_mm)
139 788 : CALL se_param_set_default(se_kind_mm, 0, do_method_pchg)
140 788 : itype = get_se_type(se_kind_mm%typ)
141 788 : nkind = SIZE(atomic_kind_set)
142 788 : enuclear = 0.0_dp
143 17300 : Forces_QM = 0.0_dp
144 788 : CALL qs_rho_get(rho, rho_ao=matrix_p)
145 :
146 1576 : DO ispin = 1, dft_control%nspins
147 : iqm = 0
148 3204 : Kinds: DO ikind = 1, nkind
149 1628 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=list)
150 1628 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
151 : CALL get_se_param(se_kind_a, &
152 : defined=defined, &
153 1628 : natorb=natorb_a)
154 1628 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
155 9800 : Atoms: DO i = 1, SIZE(list)
156 4128 : iqm = iqm + 1
157 4128 : iatom = list(i)
158 : ! Give back block
159 4128 : NULLIFY (p_block_a)
160 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
161 4128 : row=iatom, col=iatom, BLOCK=p_block_a, found=found)
162 :
163 5756 : IF (ASSOCIATED(p_block_a)) THEN
164 : ! Expand derivative of geometrical factors
165 : CALL deriv_se_qmmm_matrix_low(p_block_a, &
166 : se_kind_a, &
167 : se_kind_mm, &
168 : qmmm_env%Potentials, &
169 : particles_mm, &
170 : qmmm_env%mm_atom_chrg, &
171 : qmmm_env%mm_atom_index, &
172 : mm_cell, &
173 : iatom, &
174 : itype, &
175 : Forces, &
176 : Forces_QM(:, iqm), &
177 : se_taper, &
178 : se_int_control, &
179 : anag, &
180 : delta, &
181 : qmmm_env%spherical_cutoff, &
182 2064 : particles_qm)
183 : ! Possibly added charges
184 2064 : IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
185 : CALL deriv_se_qmmm_matrix_low(p_block_a, &
186 : se_kind_a, &
187 : se_kind_mm, &
188 : qmmm_env%added_charges%potentials, &
189 : qmmm_env%added_charges%added_particles, &
190 : qmmm_env%added_charges%mm_atom_chrg, &
191 : qmmm_env%added_charges%mm_atom_index, &
192 : mm_cell, &
193 : iatom, &
194 : itype, &
195 : Forces_added_charges, &
196 : Forces_QM(:, iqm), &
197 : se_taper, &
198 : se_int_control, &
199 : anag, &
200 : delta, &
201 : qmmm_env%spherical_cutoff, &
202 0 : particles_qm)
203 : END IF
204 : END IF
205 : END DO Atoms
206 : END DO Kinds
207 : END DO
208 788 : CPASSERT(iqm == number_qm_atoms)
209 : ! Transfer QM gradients to the QM particles..
210 33812 : CALL para_env%sum(Forces_QM)
211 788 : iqm = 0
212 2416 : DO ikind = 1, nkind
213 1628 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
214 1628 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
215 : CALL get_se_param(se_kind_a, &
216 : defined=defined, &
217 1628 : natorb=natorb_a)
218 1628 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
219 8172 : DO i = 1, SIZE(list)
220 4128 : iqm = iqm + 1
221 4128 : iatom = qmmm_env%qm_atom_index(list(i))
222 30524 : particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
223 : END DO
224 : END DO
225 : ! MM forces will be handled directly from the QMMM module in the same way
226 : ! as for GPW/GAPW methods
227 788 : DEALLOCATE (Forces_QM)
228 788 : CALL semi_empirical_release(se_kind_mm)
229 :
230 : END IF
231 936 : CALL timestop(handle)
232 936 : END SUBROUTINE deriv_se_qmmm_matrix
233 :
234 : ! **************************************************************************************************
235 : !> \brief Low Level : Computes derivatives of the 1-el semi-empirical QMMM
236 : !> hamiltonian block w.r.t. MM and QM coordinates
237 : !> \param p_block_a ...
238 : !> \param se_kind_a ...
239 : !> \param se_kind_mm ...
240 : !> \param potentials ...
241 : !> \param particles_mm ...
242 : !> \param mm_charges ...
243 : !> \param mm_atom_index ...
244 : !> \param mm_cell ...
245 : !> \param IndQM ...
246 : !> \param itype ...
247 : !> \param forces ...
248 : !> \param forces_qm ...
249 : !> \param se_taper ...
250 : !> \param se_int_control ...
251 : !> \param anag ...
252 : !> \param delta ...
253 : !> \param qmmm_spherical_cutoff ...
254 : !> \param particles_qm ...
255 : !> \author Teodoro Laino 04.2007 [created]
256 : ! **************************************************************************************************
257 4128 : SUBROUTINE deriv_se_qmmm_matrix_low(p_block_a, se_kind_a, se_kind_mm, &
258 : potentials, particles_mm, mm_charges, mm_atom_index, &
259 2064 : mm_cell, IndQM, itype, forces, forces_qm, se_taper, &
260 : se_int_control, anag, delta, qmmm_spherical_cutoff, particles_qm)
261 :
262 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block_a
263 : TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm
264 : TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
265 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
266 : REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges
267 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
268 : TYPE(cell_type), POINTER :: mm_cell
269 : INTEGER, INTENT(IN) :: IndQM, itype
270 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces
271 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: forces_qm
272 : TYPE(se_taper_type), POINTER :: se_taper
273 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
274 : LOGICAL, INTENT(IN) :: anag
275 : REAL(KIND=dp), INTENT(IN) :: delta, qmmm_spherical_cutoff(2)
276 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
277 :
278 : CHARACTER(len=*), PARAMETER :: routineN = 'deriv_se_qmmm_matrix_low'
279 :
280 : INTEGER :: handle, i1, i1L, i2, Imm, Imp, IndMM, &
281 : Ipot, j1, j1L
282 : REAL(KIND=dp) :: rt1, rt2, rt3, sph_chrg_factor
283 : REAL(KIND=dp), DIMENSION(3) :: denuc, force_ab, r_pbc, rij
284 : REAL(KIND=dp), DIMENSION(3, 45) :: de1b
285 : TYPE(qmmm_pot_type), POINTER :: Pot
286 :
287 2064 : CALL timeset(routineN, handle)
288 : ! Loop Over MM atoms - parallelization over MM atoms...
289 : ! Loop over Pot stores atoms with the same charge
290 5802 : MainLoopPot: DO Ipot = 1, SIZE(Potentials)
291 3738 : Pot => Potentials(Ipot)%Pot
292 : ! Loop over atoms belonging to this type
293 8436363 : LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index)
294 8430561 : Imm = Pot%mm_atom_index(Imp)
295 8430561 : IndMM = mm_atom_index(Imm)
296 33722244 : r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell)
297 8430561 : rt1 = r_pbc(1)
298 8430561 : rt2 = r_pbc(2)
299 8430561 : rt3 = r_pbc(3)
300 33722244 : rij = (/rt1, rt2, rt3/)
301 8430561 : se_kind_mm%zeff = mm_charges(Imm)
302 : ! Computes the screening factor for the spherical cutoff
303 8430561 : IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
304 910272 : CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
305 910272 : se_kind_mm%zeff = se_kind_mm%zeff*sph_chrg_factor
306 : END IF
307 8430561 : IF (ABS(se_kind_mm%zeff) <= EPSILON(0.0_dp)) CYCLE
308 : ! Integrals derivatives involving QM - MM atoms
309 : CALL drotnuc(se_kind_a, se_kind_mm, rij, itype=itype, de1b=de1b, &
310 : se_int_control=se_int_control, anag=anag, delta=delta, &
311 7671782 : se_taper=se_taper)
312 : CALL dcorecore(se_kind_a, se_kind_mm, rij, itype=itype, denuc=denuc, &
313 : se_int_control=se_int_control, anag=anag, delta=delta, &
314 7671782 : se_taper=se_taper)
315 : ! Nucler - Nuclear term
316 30687128 : force_ab(1:3) = -denuc(1:3)
317 : ! Force contribution from the QMMM Hamiltonian
318 7671782 : i2 = 0
319 25586785 : DO i1L = 1, se_kind_a%natorb
320 17915003 : i1 = se_orbital_pointer(i1L)
321 38401445 : DO j1L = 1, i1L - 1
322 20486442 : j1 = se_orbital_pointer(j1L)
323 20486442 : i2 = i2 + 1
324 99860771 : force_ab = force_ab - 2.0_dp*de1b(:, i2)*p_block_a(i1, j1)
325 : END DO
326 17915003 : j1 = se_orbital_pointer(j1L)
327 17915003 : i2 = i2 + 1
328 79331794 : force_ab = force_ab - de1b(:, i2)*p_block_a(i1, j1)
329 : END DO
330 : ! The array of QM forces are really the forces
331 30687128 : forces_qm(:) = forces_qm(:) - force_ab
332 : ! The one of MM atoms are instead gradients
333 30690866 : forces(:, Imm) = forces(:, Imm) - force_ab
334 : END DO LoopMM
335 : END DO MainLoopPot
336 2064 : CALL timestop(handle)
337 2064 : END SUBROUTINE deriv_se_qmmm_matrix_low
338 :
339 : END MODULE qmmm_se_forces
|
