(git:5e7fe52)
Loading...
Searching...
No Matches
qs_efield_berry.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculates the energy contribution and the mo_derivative of
10!> a static periodic electric field
11!> \par History
12!> none
13!> \author fschiff (06.2010)
14! **************************************************************************************************
15MODULE qs_efield_berry
16 USE ai_moments, ONLY: cossin
17 USE atomic_kind_types, ONLY: atomic_kind_type,&
18 get_atomic_kind,&
19 get_atomic_kind_set
20 USE basis_set_types, ONLY: gto_basis_set_p_type,&
21 gto_basis_set_type
22 USE block_p_types, ONLY: block_p_type
23 USE cell_types, ONLY: cell_type,&
24 pbc
25 USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm,&
26 cp_cfm_solve
27 USE cp_cfm_types, ONLY: cp_cfm_create,&
28 cp_cfm_release,&
29 cp_cfm_set_all,&
30 cp_cfm_type
31 USE cp_control_types, ONLY: dft_control_type
32 USE cp_dbcsr_api, ONLY: dbcsr_copy,&
33 dbcsr_get_block_p,&
34 dbcsr_p_type,&
35 dbcsr_set,&
36 dbcsr_type
37 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
38 copy_fm_to_dbcsr,&
39 cp_dbcsr_plus_fm_fm_t,&
40 cp_dbcsr_sm_fm_multiply,&
41 dbcsr_deallocate_matrix_set
42 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
43 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
44 cp_fm_struct_release,&
45 cp_fm_struct_type
46 USE cp_fm_types, ONLY: cp_fm_create,&
47 cp_fm_release,&
48 cp_fm_set_all,&
49 cp_fm_type
50 USE kinds, ONLY: dp
51 USE mathconstants, ONLY: gaussi,&
52 pi,&
53 twopi,&
54 z_one,&
55 z_zero
56 USE message_passing, ONLY: mp_para_env_type
57 USE orbital_pointers, ONLY: ncoset
58 USE parallel_gemm_api, ONLY: parallel_gemm
59 USE particle_types, ONLY: particle_type
60 USE qs_energy_types, ONLY: qs_energy_type
61 USE qs_environment_types, ONLY: get_qs_env,&
62 qs_environment_type,&
63 set_qs_env
64 USE qs_force_types, ONLY: qs_force_type
65 USE qs_kind_types, ONLY: get_qs_kind,&
66 get_qs_kind_set,&
67 qs_kind_type
68 USE qs_mo_types, ONLY: get_mo_set,&
69 mo_set_type
70 USE qs_moments, ONLY: build_berry_moment_matrix
71 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
72 neighbor_list_iterate,&
73 neighbor_list_iterator_create,&
74 neighbor_list_iterator_p_type,&
75 neighbor_list_iterator_release,&
76 neighbor_list_set_p_type
77 USE qs_period_efield_types, ONLY: efield_berry_type,&
78 init_efield_matrices,&
79 set_efield_matrices
80 USE virial_methods, ONLY: virial_pair_force
81 USE virial_types, ONLY: virial_type
82#include "./base/base_uses.f90"
83
84 IMPLICIT NONE
85
86 PRIVATE
87
88 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_efield_berry'
89
90 ! *** Public subroutines ***
91
92 PUBLIC :: qs_efield_berry_phase
93
94! **************************************************************************************************
95
96CONTAINS
97
98! **************************************************************************************************
99
100! **************************************************************************************************
101!> \brief ...
102!> \param qs_env ...
103!> \param just_energy ...
104!> \param calculate_forces ...
105! **************************************************************************************************
106 SUBROUTINE qs_efield_berry_phase(qs_env, just_energy, calculate_forces)
107
108 TYPE(qs_environment_type), POINTER :: qs_env
109 LOGICAL, INTENT(IN) :: just_energy, calculate_forces
110
111 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_efield_berry_phase'
112
113 INTEGER :: handle
114 LOGICAL :: s_mstruct_changed
115 TYPE(dft_control_type), POINTER :: dft_control
116
117 CALL timeset(routinen, handle)
118
119 NULLIFY (dft_control)
120 CALL get_qs_env(qs_env, s_mstruct_changed=s_mstruct_changed, &
121 dft_control=dft_control)
122
123 IF (dft_control%apply_period_efield) THEN
124 ! check if the periodic efield should be applied in the current step
125 IF (dft_control%period_efield%start_frame <= qs_env%sim_step .AND. &
126 (dft_control%period_efield%end_frame == -1 .OR. dft_control%period_efield%end_frame >= qs_env%sim_step)) THEN
127
128 IF (s_mstruct_changed) CALL qs_efield_integrals(qs_env)
129 IF (dft_control%period_efield%displacement_field) THEN
130 CALL qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
131 ELSE
132 CALL qs_efield_derivatives(qs_env, just_energy, calculate_forces)
133 END IF
134 END IF
135 END IF
136
137 CALL timestop(handle)
138
139 END SUBROUTINE qs_efield_berry_phase
140
141! **************************************************************************************************
142!> \brief ...
143!> \param qs_env ...
144! **************************************************************************************************
145 SUBROUTINE qs_efield_integrals(qs_env)
146
147 TYPE(qs_environment_type), POINTER :: qs_env
148
149 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_efield_integrals'
150
151 INTEGER :: handle, i
152 REAL(dp), DIMENSION(3) :: kvec
153 TYPE(cell_type), POINTER :: cell
154 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: cosmat, matrix_s, sinmat
155 TYPE(dft_control_type), POINTER :: dft_control
156 TYPE(efield_berry_type), POINTER :: efield
157
158 CALL timeset(routinen, handle)
159 cpassert(ASSOCIATED(qs_env))
160
161 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
162 NULLIFY (matrix_s)
163 CALL get_qs_env(qs_env=qs_env, efield=efield, cell=cell, matrix_s=matrix_s)
164 CALL init_efield_matrices(efield)
165 ALLOCATE (cosmat(3), sinmat(3))
166 DO i = 1, 3
167 ALLOCATE (cosmat(i)%matrix, sinmat(i)%matrix)
168
169 CALL dbcsr_copy(cosmat(i)%matrix, matrix_s(1)%matrix, 'COS MAT')
170 CALL dbcsr_copy(sinmat(i)%matrix, matrix_s(1)%matrix, 'SIN MAT')
171
172 kvec(:) = twopi*cell%h_inv(i, :)
173 CALL build_berry_moment_matrix(qs_env, cosmat(i)%matrix, sinmat(i)%matrix, kvec)
174 END DO
175 CALL set_efield_matrices(efield=efield, cosmat=cosmat, sinmat=sinmat)
176 CALL set_qs_env(qs_env=qs_env, efield=efield)
177 CALL timestop(handle)
178
179 END SUBROUTINE qs_efield_integrals
180
181! **************************************************************************************************
182!> \brief ...
183!> \param qs_env ...
184!> \param just_energy ...
185!> \param calculate_forces ...
186! **************************************************************************************************
187 SUBROUTINE qs_efield_derivatives(qs_env, just_energy, calculate_forces)
188 TYPE(qs_environment_type), POINTER :: qs_env
189 LOGICAL, INTENT(IN) :: just_energy, calculate_forces
190
191 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_efield_derivatives'
192
193 COMPLEX(dp) :: zdet, zdeta, zi(3)
194 INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, j, &
195 jatom, jkind, jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, &
196 nseta, nsetb, sgfa, sgfb
197 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
198 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
199 npgfb, nsgfa, nsgfb
200 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
201 LOGICAL :: found, uniform, use_virial
202 REAL(dp) :: charge, ci(3), cqi(3), dab, dd, &
203 ener_field, f0, fab, fieldpol(3), &
204 focc, fpolvec(3), hmat(3, 3), occ, &
205 qi(3), strength, ti(3)
206 REAL(dp), DIMENSION(3) :: forcea, forceb, kvec, ra, rab, rb, ria
207 REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work
208 REAL(dp), DIMENSION(:, :, :), POINTER :: dcosab, dsinab
209 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
210 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
211 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
212 TYPE(block_p_type), DIMENSION(3, 2) :: dcost, dsint
213 TYPE(cell_type), POINTER :: cell
214 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat, inv_mat
215 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
216 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff_tmp, mo_derivs_tmp
217 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: inv_work, op_fm_set, opvec
218 TYPE(cp_fm_type), POINTER :: mo_coeff
219 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, mo_derivs
220 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: tempmat
221 TYPE(dbcsr_type), POINTER :: cosmat, mo_coeff_b, sinmat
222 TYPE(dft_control_type), POINTER :: dft_control
223 TYPE(efield_berry_type), POINTER :: efield
224 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
225 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
226 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
227 TYPE(mp_para_env_type), POINTER :: para_env
228 TYPE(neighbor_list_iterator_p_type), &
229 DIMENSION(:), POINTER :: nl_iterator
230 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
231 POINTER :: sab_orb
232 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
233 TYPE(qs_energy_type), POINTER :: energy
234 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
235 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
236 TYPE(qs_kind_type), POINTER :: qs_kind
237 TYPE(virial_type), POINTER :: virial
238
239 CALL timeset(routinen, handle)
240
241 NULLIFY (dft_control, cell, particle_set)
242 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
243 particle_set=particle_set, virial=virial)
244 NULLIFY (qs_kind_set, efield, para_env, sab_orb)
245 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
246 efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
247
248 ! calculate stress only if forces requested also
249 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
250 use_virial = use_virial .AND. calculate_forces
251 ! disable stress calculation
252 IF (use_virial) THEN
253 cpabort("Stress tensor for periodic E-field not implemented")
254 END IF
255
256 ! if an intensities list is given, select the value for the current step
257 strength = dft_control%period_efield%strength
258 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
259 strength = dft_control%period_efield%strength_list(mod(qs_env%sim_step &
260 - dft_control%period_efield%start_frame, SIZE(dft_control%period_efield%strength_list)) + 1)
261 END IF
262
263 fieldpol = dft_control%period_efield%polarisation
264 fieldpol = fieldpol/norm2(fieldpol)
265 fieldpol = -fieldpol*strength
266 hmat = cell%hmat(:, :)/twopi
267 DO idir = 1, 3
268 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
269 END DO
270
271 ! nuclear contribution
272 natom = SIZE(particle_set)
273 IF (calculate_forces) THEN
274 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
275 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
276 END IF
277 zi(:) = cmplx(1._dp, 0._dp, dp)
278 DO ia = 1, natom
279 CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
280 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
281 ria = particle_set(ia)%r
282 ria = pbc(ria, cell)
283 DO idir = 1, 3
284 kvec(:) = twopi*cell%h_inv(idir, :)
285 dd = sum(kvec(:)*ria(:))
286 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
287 zi(idir) = zi(idir)*zdeta
288 END DO
289 IF (calculate_forces) THEN
290 IF (para_env%mepos == 0) THEN
291 iatom = atom_of_kind(ia)
292 forcea(:) = fieldpol(:)*charge
293 force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + forcea(:)
294 END IF
295 END IF
296 IF (use_virial) THEN
297 IF (para_env%mepos == 0) THEN
298 CALL virial_pair_force(virial%pv_virial, 1.0_dp, forcea, ria)
299 END IF
300 END IF
301 END DO
302 qi = aimag(log(zi))
303
304 ! check uniform occupation
305 NULLIFY (mos)
306 CALL get_qs_env(qs_env=qs_env, mos=mos)
307 DO ispin = 1, dft_control%nspins
308 CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
309 IF (.NOT. uniform) THEN
310 cpabort("Berry phase moments for non uniform MOs' occupation numbers not implemented")
311 END IF
312 END DO
313
314 NULLIFY (mo_derivs)
315 CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs)
316 ! initialize all work matrices needed
317 ALLOCATE (op_fm_set(2, dft_control%nspins))
318 ALLOCATE (opvec(2, dft_control%nspins))
319 ALLOCATE (eigrmat(dft_control%nspins))
320 ALLOCATE (inv_mat(dft_control%nspins))
321 ALLOCATE (inv_work(2, dft_control%nspins))
322 ALLOCATE (mo_derivs_tmp(SIZE(mo_derivs)))
323 ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs)))
324
325 ! Allocate temp matrices for the wavefunction derivatives
326 DO ispin = 1, dft_control%nspins
327 NULLIFY (tmp_fm_struct, mo_coeff)
328 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
329 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
330 ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
331 CALL cp_fm_create(mo_derivs_tmp(ispin), mo_coeff%matrix_struct)
332 CALL cp_fm_create(mo_coeff_tmp(ispin), mo_coeff%matrix_struct)
333 CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_derivs_tmp(ispin))
334 DO i = 1, SIZE(op_fm_set, 1)
335 CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
336 CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
337 CALL cp_fm_create(inv_work(i, ispin), op_fm_set(i, ispin)%matrix_struct)
338 END DO
339 CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
340 CALL cp_cfm_create(inv_mat(ispin), op_fm_set(1, ispin)%matrix_struct)
341 CALL cp_fm_struct_release(tmp_fm_struct)
342 END DO
343 ! temp matrices for force calculation
344 IF (calculate_forces) THEN
345 NULLIFY (matrix_s)
346 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
347 ALLOCATE (tempmat(2, dft_control%nspins))
348 DO ispin = 1, dft_control%nspins
349 ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix)
350 CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
351 CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
352 CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
353 CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
354 END DO
355 ! integration
356 CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab)
357 ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab))
358 ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3))
359 lsab = max(ldab, lsab)
360 DO i = 1, 3
361 ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab))
362 ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab))
363 END DO
364 END IF
365
366 !Start the MO derivative calculation
367 !loop over all cell vectors
368 DO idir = 1, 3
369 ci(idir) = 0.0_dp
370 zi(idir) = z_zero
371 IF (abs(fpolvec(idir)) > 1.0e-12_dp) THEN
372 cosmat => efield%cosmat(idir)%matrix
373 sinmat => efield%sinmat(idir)%matrix
374 !evaluate the expression needed for the derivative (S_berry * C and [C^T S_berry C]^-1)
375 !first step S_berry * C and C^T S_berry C
376 DO ispin = 1, dft_control%nspins ! spin
377 IF (mos(ispin)%use_mo_coeff_b) THEN
378 CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo)
379 CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin))
380 ELSE
381 CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
382 mo_coeff_tmp(ispin) = mo_coeff
383 END IF
384 CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin), opvec(1, ispin), ncol=nmo)
385 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(1, ispin), 0.0_dp, &
386 op_fm_set(1, ispin))
387 CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin), opvec(2, ispin), ncol=nmo)
388 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(2, ispin), 0.0_dp, &
389 op_fm_set(2, ispin))
390 END DO
391 !second step invert C^T S_berry C
392 zdet = z_one
393 DO ispin = 1, dft_control%nspins
394 CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin), z_one, op_fm_set(1, ispin))
395 CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin), -gaussi, op_fm_set(2, ispin))
396 CALL cp_cfm_set_all(inv_mat(ispin), z_zero, z_one)
397 CALL cp_cfm_solve(eigrmat(ispin), inv_mat(ispin), zdeta)
398 zdet = zdet*zdeta
399 END DO
400 zi(idir) = zdet**occ
401 ci(idir) = aimag(log(zdet**occ))
402
403 IF (.NOT. just_energy) THEN
404 !compute the orbital derivative
405 focc = fpolvec(idir)
406 DO ispin = 1, dft_control%nspins
407 inv_work(1, ispin)%local_data(:, :) = real(inv_mat(ispin)%local_data(:, :), dp)
408 inv_work(2, ispin)%local_data(:, :) = aimag(inv_mat(ispin)%local_data(:, :))
409 CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
410 CALL parallel_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin), inv_work(2, ispin), &
411 1.0_dp, mo_derivs_tmp(ispin))
412 CALL parallel_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin), inv_work(1, ispin), &
413 1.0_dp, mo_derivs_tmp(ispin))
414 END DO
415 END IF
416
417 !compute nuclear forces
418 IF (calculate_forces) THEN
419 nkind = SIZE(qs_kind_set)
420 natom = SIZE(particle_set)
421 kvec(:) = twopi*cell%h_inv(idir, :)
422
423 ! calculate: C [C^T S_berry C]^(-1) C^T
424 ! Store this matrix in DBCSR form (only S overlap blocks)
425 DO ispin = 1, dft_control%nspins
426 CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
427 CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
428 CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
429 CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(1, ispin), 0.0_dp, &
430 opvec(1, ispin))
431 CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(2, ispin), 0.0_dp, &
432 opvec(2, ispin))
433 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, &
434 matrix_v=opvec(1, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
435 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, &
436 matrix_v=opvec(2, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
437 END DO
438
439 ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db)
440 ALLOCATE (basis_set_list(nkind))
441 DO ikind = 1, nkind
442 qs_kind => qs_kind_set(ikind)
443 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
444 IF (ASSOCIATED(basis_set_a)) THEN
445 basis_set_list(ikind)%gto_basis_set => basis_set_a
446 ELSE
447 NULLIFY (basis_set_list(ikind)%gto_basis_set)
448 END IF
449 END DO
450 !
451 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
452 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
453 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
454 iatom=iatom, jatom=jatom, r=rab)
455 basis_set_a => basis_set_list(ikind)%gto_basis_set
456 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
457 basis_set_b => basis_set_list(jkind)%gto_basis_set
458 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
459 ! basis ikind
460 first_sgfa => basis_set_a%first_sgf
461 la_max => basis_set_a%lmax
462 la_min => basis_set_a%lmin
463 npgfa => basis_set_a%npgf
464 nseta = basis_set_a%nset
465 nsgfa => basis_set_a%nsgf_set
466 rpgfa => basis_set_a%pgf_radius
467 set_radius_a => basis_set_a%set_radius
468 sphi_a => basis_set_a%sphi
469 zeta => basis_set_a%zet
470 ! basis jkind
471 first_sgfb => basis_set_b%first_sgf
472 lb_max => basis_set_b%lmax
473 lb_min => basis_set_b%lmin
474 npgfb => basis_set_b%npgf
475 nsetb = basis_set_b%nset
476 nsgfb => basis_set_b%nsgf_set
477 rpgfb => basis_set_b%pgf_radius
478 set_radius_b => basis_set_b%set_radius
479 sphi_b => basis_set_b%sphi
480 zetb => basis_set_b%zet
481
482 atom_a = atom_of_kind(iatom)
483 atom_b = atom_of_kind(jatom)
484
485 ldsa = SIZE(sphi_a, 1)
486 ldsb = SIZE(sphi_b, 1)
487 ra(:) = pbc(particle_set(iatom)%r(:), cell)
488 rb(:) = ra + rab
489 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
490
491 IF (iatom <= jatom) THEN
492 irow = iatom
493 icol = jatom
494 ELSE
495 irow = jatom
496 icol = iatom
497 END IF
498
499 IF (iatom == jatom) THEN
500 fab = 1.0_dp*occ
501 ELSE
502 fab = 2.0_dp*occ
503 END IF
504
505 DO i = 1, 3
506 dcost(i, 1)%block = 0.0_dp
507 dsint(i, 1)%block = 0.0_dp
508 dcost(i, 2)%block = 0.0_dp
509 dsint(i, 2)%block = 0.0_dp
510 END DO
511
512 DO iset = 1, nseta
513 ncoa = npgfa(iset)*ncoset(la_max(iset))
514 sgfa = first_sgfa(1, iset)
515 DO jset = 1, nsetb
516 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
517 ncob = npgfb(jset)*ncoset(lb_max(jset))
518 sgfb = first_sgfb(1, jset)
519 ! Calculate the primitive integrals (da|b)
520 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
521 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
522 ra, rb, kvec, cosab, sinab, dcosab, dsinab)
523 DO i = 1, 3
524 CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, &
525 ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
526 ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
527 dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
528 END DO
529 ! Calculate the primitive integrals (a|db)
530 CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
531 la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
532 rb, ra, kvec, cosab, sinab, dcosab, dsinab)
533 DO i = 1, 3
534 dcosab(1:ncoa, 1:ncob, i) = transpose(dcosab(1:ncob, 1:ncoa, i))
535 dsinab(1:ncoa, 1:ncob, i) = transpose(dsinab(1:ncob, 1:ncoa, i))
536 CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, &
537 ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
538 ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
539 dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
540 END DO
541 END DO
542 END DO
543 forcea = 0.0_dp
544 forceb = 0.0_dp
545 DO ispin = 1, dft_control%nspins
546 NULLIFY (rblock, iblock)
547 CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, &
548 row=irow, col=icol, block=rblock, found=found)
549 cpassert(found)
550 CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, &
551 row=irow, col=icol, block=iblock, found=found)
552 cpassert(found)
553 n1 = SIZE(rblock, 1)
554 n2 = SIZE(rblock, 2)
555 cpassert(SIZE(iblock, 1) == n1)
556 cpassert(SIZE(iblock, 2) == n2)
557 cpassert(lsab >= n1)
558 cpassert(lsab >= n2)
559 IF (iatom <= jatom) THEN
560 DO i = 1, 3
561 forcea(i) = forcea(i) + sum(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) &
562 - sum(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2))
563 forceb(i) = forceb(i) + sum(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) &
564 - sum(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2))
565 END DO
566 ELSE
567 DO i = 1, 3
568 forcea(i) = forcea(i) + sum(transpose(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) &
569 - sum(transpose(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1))
570 forceb(i) = forceb(i) + sum(transpose(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) &
571 - sum(transpose(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1))
572 END DO
573 END IF
574 END DO
575 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fab*fpolvec(idir)*forcea(1:3)
576 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fab*fpolvec(idir)*forceb(1:3)
577 IF (use_virial) THEN
578 f0 = -fab*fpolvec(idir)
579 CALL virial_pair_force(virial%pv_virial, f0, forcea, ra)
580 CALL virial_pair_force(virial%pv_virial, f0, forceb, rb)
581 END IF
582
583 END DO
584 CALL neighbor_list_iterator_release(nl_iterator)
585 DEALLOCATE (basis_set_list)
586
587 END IF
588 END IF
589 END DO
590
591 ! Energy
592 ener_field = 0.0_dp
593 ti = 0.0_dp
594 DO idir = 1, 3
595 ! make sure the total normalized polarization is within [-1:1]
596 cqi(idir) = qi(idir) + ci(idir)
597 IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
598 IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
599 ! now check for log branch
600 IF (abs(efield%polarisation(idir) - cqi(idir)) > pi) THEN
601 ti(idir) = (efield%polarisation(idir) - cqi(idir))/pi
602 DO i = 1, 10
603 cqi(idir) = cqi(idir) + sign(1.0_dp, ti(idir))*twopi
604 IF (abs(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
605 END DO
606 END IF
607 ener_field = ener_field + fpolvec(idir)*cqi(idir)
608 END DO
609
610 ! update the references
611 IF (calculate_forces) THEN
612 ! check for smoothness of energy surface
613 IF (abs(efield%field_energy - ener_field) > pi*abs(sum(fpolvec))) THEN
614 cpwarn("Large change of e-field energy detected. Correct for non-smooth energy surface")
615 END IF
616 efield%field_energy = ener_field
617 efield%polarisation(:) = cqi(:)
618 END IF
619 energy%efield = ener_field
620
621 IF (.NOT. just_energy) THEN
622 ! Add the result to mo_derivativs
623 DO ispin = 1, dft_control%nspins
624 CALL copy_fm_to_dbcsr(mo_derivs_tmp(ispin), mo_derivs(ispin)%matrix)
625 END DO
626 IF (use_virial) THEN
627 ti = 0.0_dp
628 DO i = 1, 3
629 DO j = 1, 3
630 ti(j) = ti(j) + hmat(j, i)*cqi(i)
631 END DO
632 END DO
633 DO i = 1, 3
634 DO j = 1, 3
635 virial%pv_virial(i, j) = virial%pv_virial(i, j) - fieldpol(i)*ti(j)
636 END DO
637 END DO
638 END IF
639 END IF
640
641 DO ispin = 1, dft_control%nspins
642 CALL cp_cfm_release(eigrmat(ispin))
643 CALL cp_cfm_release(inv_mat(ispin))
644 CALL cp_fm_release(mo_derivs_tmp(ispin))
645 IF (mos(ispin)%use_mo_coeff_b) CALL cp_fm_release(mo_coeff_tmp(ispin))
646 DO i = 1, SIZE(op_fm_set, 1)
647 CALL cp_fm_release(opvec(i, ispin))
648 CALL cp_fm_release(op_fm_set(i, ispin))
649 CALL cp_fm_release(inv_work(i, ispin))
650 END DO
651 END DO
652 DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat)
653 DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp)
654
655 IF (calculate_forces) THEN
656 DO ikind = 1, SIZE(atomic_kind_set)
657 CALL para_env%sum(force(ikind)%efield)
658 END DO
659 DEALLOCATE (cosab, sinab, work, dcosab, dsinab)
660 DO i = 1, 3
661 DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block)
662 DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block)
663 END DO
664 CALL dbcsr_deallocate_matrix_set(tempmat)
665 END IF
666 CALL timestop(handle)
667
668 END SUBROUTINE qs_efield_derivatives
669
670! **************************************************************************************************
671!> \brief ...
672!> \param qs_env ...
673!> \param just_energy ...
674!> \param calculate_forces ...
675! **************************************************************************************************
676 SUBROUTINE qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
677 TYPE(qs_environment_type), POINTER :: qs_env
678 LOGICAL, INTENT(IN) :: just_energy, calculate_forces
679
680 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_dispfield_derivatives'
681
682 COMPLEX(dp) :: zdet, zdeta, zi(3)
683 INTEGER :: handle, i, ia, iatom, icol, idir, ikind, iodeb, irow, iset, ispin, jatom, jkind, &
684 jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, nseta, nsetb, &
685 sgfa, sgfb
686 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
687 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
688 npgfb, nsgfa, nsgfb
689 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
690 LOGICAL :: found, uniform, use_virial
691 REAL(dp) :: charge, ci(3), cqi(3), dab, dd, di(3), ener_field, fab, fieldpol(3), focc, &
692 hmat(3, 3), occ, omega, qi(3), rlog(3), strength, zlog(3)
693 REAL(dp), DIMENSION(3) :: dfilter, forcea, forceb, kvec, ra, rab, &
694 rb, ria
695 REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work
696 REAL(dp), DIMENSION(:, :, :), POINTER :: dcosab, dsinab, force_tmp
697 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
698 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
699 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
700 TYPE(block_p_type), DIMENSION(3, 2) :: dcost, dsint
701 TYPE(cell_type), POINTER :: cell
702 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat, inv_mat
703 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
704 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff_tmp
705 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: inv_work, mo_derivs_tmp, op_fm_set, opvec
706 TYPE(cp_fm_type), POINTER :: mo_coeff
707 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, mo_derivs
708 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: tempmat
709 TYPE(dbcsr_type), POINTER :: cosmat, mo_coeff_b, sinmat
710 TYPE(dft_control_type), POINTER :: dft_control
711 TYPE(efield_berry_type), POINTER :: efield
712 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
713 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
714 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
715 TYPE(mp_para_env_type), POINTER :: para_env
716 TYPE(neighbor_list_iterator_p_type), &
717 DIMENSION(:), POINTER :: nl_iterator
718 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
719 POINTER :: sab_orb
720 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
721 TYPE(qs_energy_type), POINTER :: energy
722 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
723 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
724 TYPE(qs_kind_type), POINTER :: qs_kind
725 TYPE(virial_type), POINTER :: virial
726
727 CALL timeset(routinen, handle)
728
729 NULLIFY (dft_control, cell, particle_set)
730 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
731 particle_set=particle_set, virial=virial)
732 NULLIFY (qs_kind_set, efield, para_env, sab_orb)
733 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
734 efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
735
736 ! calculate stress only if forces requested also
737 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
738 use_virial = use_virial .AND. calculate_forces
739 ! disable stress calculation
740 IF (use_virial) THEN
741 cpabort("Stress tensor for periodic D-field not implemented")
742 END IF
743
744 dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
745
746 ! if an intensities list is given, select the value for the current step
747 strength = dft_control%period_efield%strength
748 IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
749 strength = dft_control%period_efield%strength_list(mod(qs_env%sim_step &
750 - dft_control%period_efield%start_frame, SIZE(dft_control%period_efield%strength_list)) + 1)
751 END IF
752
753 fieldpol = dft_control%period_efield%polarisation
754 fieldpol = fieldpol/norm2(fieldpol)
755 fieldpol = fieldpol*strength
756
757 omega = cell%deth
758 hmat = cell%hmat(:, :)/(twopi*omega)
759
760 ! nuclear contribution to polarization
761 natom = SIZE(particle_set)
762 IF (calculate_forces) THEN
763 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
764 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
765 ALLOCATE (force_tmp(natom, 3, 3))
766 force_tmp = 0.0_dp
767 END IF
768 zi(:) = cmplx(1._dp, 0._dp, dp)
769 DO ia = 1, natom
770 CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
771 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
772 ria = particle_set(ia)%r
773 ria = pbc(ria, cell)
774 DO idir = 1, 3
775 kvec(:) = twopi*cell%h_inv(idir, :)
776 dd = sum(kvec(:)*ria(:))
777 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
778 zi(idir) = zi(idir)*zdeta
779 END DO
780 IF (calculate_forces) THEN
781 IF (para_env%mepos == 0) THEN
782 DO i = 1, 3
783 force_tmp(ia, i, i) = force_tmp(ia, i, i) + charge/omega
784 END DO
785 END IF
786 END IF
787 END DO
788 rlog = aimag(log(zi))
789
790 ! check uniform occupation
791 NULLIFY (mos)
792 CALL get_qs_env(qs_env=qs_env, mos=mos)
793 DO ispin = 1, dft_control%nspins
794 CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
795 IF (.NOT. uniform) THEN
796 cpabort("Berry phase moments for non uniform MO occupation numbers not implemented")
797 END IF
798 END DO
799
800 ! initialize all work matrices needed
801 NULLIFY (mo_derivs)
802 CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs)
803 ALLOCATE (op_fm_set(2, dft_control%nspins))
804 ALLOCATE (opvec(2, dft_control%nspins))
805 ALLOCATE (eigrmat(dft_control%nspins))
806 ALLOCATE (inv_mat(dft_control%nspins))
807 ALLOCATE (inv_work(2, dft_control%nspins))
808 ALLOCATE (mo_derivs_tmp(3, SIZE(mo_derivs)))
809 ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs)))
810
811 ! Allocate temp matrices for the wavefunction derivatives
812 DO ispin = 1, dft_control%nspins
813 NULLIFY (tmp_fm_struct, mo_coeff)
814 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
815 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
816 ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
817 CALL cp_fm_create(mo_coeff_tmp(ispin), mo_coeff%matrix_struct)
818 DO i = 1, 3
819 CALL cp_fm_create(mo_derivs_tmp(i, ispin), mo_coeff%matrix_struct)
820 CALL cp_fm_set_all(matrix=mo_derivs_tmp(i, ispin), alpha=0.0_dp)
821 END DO
822 DO i = 1, SIZE(op_fm_set, 1)
823 CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
824 CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
825 CALL cp_fm_create(inv_work(i, ispin), op_fm_set(i, ispin)%matrix_struct)
826 END DO
827 CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
828 CALL cp_cfm_create(inv_mat(ispin), op_fm_set(1, ispin)%matrix_struct)
829 CALL cp_fm_struct_release(tmp_fm_struct)
830 END DO
831 ! temp matrices for force calculation
832 IF (calculate_forces) THEN
833 NULLIFY (matrix_s)
834 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
835 ALLOCATE (tempmat(2, dft_control%nspins))
836 DO ispin = 1, dft_control%nspins
837 ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix)
838 CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
839 CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
840 CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
841 CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
842 END DO
843 ! integration
844 CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab)
845 ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab))
846 ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3))
847 lsab = max(lsab, ldab)
848 DO i = 1, 3
849 ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab))
850 ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab))
851 END DO
852 END IF
853
854 !Start the MO derivative calculation
855 !loop over all cell vectors
856 DO idir = 1, 3
857 zi(idir) = z_zero
858 cosmat => efield%cosmat(idir)%matrix
859 sinmat => efield%sinmat(idir)%matrix
860 !evaluate the expression needed for the derivative (S_berry * C and [C^T S_berry C]^-1)
861 !first step S_berry * C and C^T S_berry C
862 DO ispin = 1, dft_control%nspins ! spin
863 IF (mos(ispin)%use_mo_coeff_b) THEN
864 CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo)
865 CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin))
866 ELSE
867 CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
868 mo_coeff_tmp(ispin) = mo_coeff
869 END IF
870 CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin), opvec(1, ispin), ncol=nmo)
871 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(1, ispin), 0.0_dp, &
872 op_fm_set(1, ispin))
873 CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin), opvec(2, ispin), ncol=nmo)
874 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(2, ispin), 0.0_dp, &
875 op_fm_set(2, ispin))
876 END DO
877 !second step invert C^T S_berry C
878 zdet = z_one
879 DO ispin = 1, dft_control%nspins
880 CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin), z_one, op_fm_set(1, ispin))
881 CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin), -gaussi, op_fm_set(2, ispin))
882 CALL cp_cfm_set_all(inv_mat(ispin), z_zero, z_one)
883 CALL cp_cfm_solve(eigrmat(ispin), inv_mat(ispin), zdeta)
884 zdet = zdet*zdeta
885 END DO
886 zi(idir) = zdet**occ
887 zlog(idir) = aimag(log(zi(idir)))
888
889 IF (.NOT. just_energy) THEN
890 !compute the orbital derivative
891 DO ispin = 1, dft_control%nspins
892 inv_work(1, ispin)%local_data(:, :) = real(inv_mat(ispin)%local_data(:, :), dp)
893 inv_work(2, ispin)%local_data(:, :) = aimag(inv_mat(ispin)%local_data(:, :))
894 CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
895 DO i = 1, 3
896 focc = hmat(idir, i)
897 CALL parallel_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin), inv_work(2, ispin), &
898 1.0_dp, mo_derivs_tmp(idir, ispin))
899 CALL parallel_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin), inv_work(1, ispin), &
900 1.0_dp, mo_derivs_tmp(idir, ispin))
901 END DO
902 END DO
903 END IF
904
905 !compute nuclear forces
906 IF (calculate_forces) THEN
907 nkind = SIZE(qs_kind_set)
908 natom = SIZE(particle_set)
909 kvec(:) = twopi*cell%h_inv(idir, :)
910
911 ! calculate: C [C^T S_berry C]^(-1) C^T
912 ! Store this matrix in DBCSR form (only S overlap blocks)
913 DO ispin = 1, dft_control%nspins
914 CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
915 CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
916 CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
917 CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(1, ispin), 0.0_dp, &
918 opvec(1, ispin))
919 CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(2, ispin), 0.0_dp, &
920 opvec(2, ispin))
921 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, &
922 matrix_v=opvec(1, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
923 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, &
924 matrix_v=opvec(2, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
925 END DO
926
927 ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db)
928 ALLOCATE (basis_set_list(nkind))
929 DO ikind = 1, nkind
930 qs_kind => qs_kind_set(ikind)
931 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
932 IF (ASSOCIATED(basis_set_a)) THEN
933 basis_set_list(ikind)%gto_basis_set => basis_set_a
934 ELSE
935 NULLIFY (basis_set_list(ikind)%gto_basis_set)
936 END IF
937 END DO
938 !
939 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
940 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
941 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
942 iatom=iatom, jatom=jatom, r=rab)
943 basis_set_a => basis_set_list(ikind)%gto_basis_set
944 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
945 basis_set_b => basis_set_list(jkind)%gto_basis_set
946 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
947 ! basis ikind
948 first_sgfa => basis_set_a%first_sgf
949 la_max => basis_set_a%lmax
950 la_min => basis_set_a%lmin
951 npgfa => basis_set_a%npgf
952 nseta = basis_set_a%nset
953 nsgfa => basis_set_a%nsgf_set
954 rpgfa => basis_set_a%pgf_radius
955 set_radius_a => basis_set_a%set_radius
956 sphi_a => basis_set_a%sphi
957 zeta => basis_set_a%zet
958 ! basis jkind
959 first_sgfb => basis_set_b%first_sgf
960 lb_max => basis_set_b%lmax
961 lb_min => basis_set_b%lmin
962 npgfb => basis_set_b%npgf
963 nsetb = basis_set_b%nset
964 nsgfb => basis_set_b%nsgf_set
965 rpgfb => basis_set_b%pgf_radius
966 set_radius_b => basis_set_b%set_radius
967 sphi_b => basis_set_b%sphi
968 zetb => basis_set_b%zet
969
970 ldsa = SIZE(sphi_a, 1)
971 ldsb = SIZE(sphi_b, 1)
972 ra(:) = pbc(particle_set(iatom)%r(:), cell)
973 rb(:) = ra + rab
974 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
975
976 IF (iatom <= jatom) THEN
977 irow = iatom
978 icol = jatom
979 ELSE
980 irow = jatom
981 icol = iatom
982 END IF
983
984 IF (iatom == jatom) THEN
985 fab = 1.0_dp*occ
986 ELSE
987 fab = 2.0_dp*occ
988 END IF
989
990 DO i = 1, 3
991 dcost(i, 1)%block = 0.0_dp
992 dsint(i, 1)%block = 0.0_dp
993 dcost(i, 2)%block = 0.0_dp
994 dsint(i, 2)%block = 0.0_dp
995 END DO
996
997 DO iset = 1, nseta
998 ncoa = npgfa(iset)*ncoset(la_max(iset))
999 sgfa = first_sgfa(1, iset)
1000 DO jset = 1, nsetb
1001 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1002 ncob = npgfb(jset)*ncoset(lb_max(jset))
1003 sgfb = first_sgfb(1, jset)
1004 ! Calculate the primitive integrals (da|b)
1005 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1006 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1007 ra, rb, kvec, cosab, sinab, dcosab, dsinab)
1008 DO i = 1, 3
1009 CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, &
1010 ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1011 ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1012 dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
1013 END DO
1014 ! Calculate the primitive integrals (a|db)
1015 CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1016 la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1017 rb, ra, kvec, cosab, sinab, dcosab, dsinab)
1018 DO i = 1, 3
1019 dcosab(1:ncoa, 1:ncob, i) = transpose(dcosab(1:ncob, 1:ncoa, i))
1020 dsinab(1:ncoa, 1:ncob, i) = transpose(dsinab(1:ncob, 1:ncoa, i))
1021 CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, &
1022 ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1023 ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1024 dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
1025 END DO
1026 END DO
1027 END DO
1028 forcea = 0.0_dp
1029 forceb = 0.0_dp
1030 DO ispin = 1, dft_control%nspins
1031 NULLIFY (rblock, iblock)
1032 CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, &
1033 row=irow, col=icol, block=rblock, found=found)
1034 cpassert(found)
1035 CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, &
1036 row=irow, col=icol, block=iblock, found=found)
1037 cpassert(found)
1038 n1 = SIZE(rblock, 1)
1039 n2 = SIZE(rblock, 2)
1040 cpassert(SIZE(iblock, 1) == n1)
1041 cpassert(SIZE(iblock, 2) == n2)
1042 cpassert(lsab >= n1)
1043 cpassert(lsab >= n2)
1044 IF (iatom <= jatom) THEN
1045 DO i = 1, 3
1046 forcea(i) = forcea(i) + sum(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) &
1047 - sum(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2))
1048 forceb(i) = forceb(i) + sum(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) &
1049 - sum(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2))
1050 END DO
1051 ELSE
1052 DO i = 1, 3
1053 forcea(i) = forcea(i) + sum(transpose(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) &
1054 - sum(transpose(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1))
1055 forceb(i) = forceb(i) + sum(transpose(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) &
1056 - sum(transpose(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1))
1057 END DO
1058 END IF
1059 END DO
1060 DO i = 1, 3
1061 force_tmp(iatom, :, i) = force_tmp(iatom, :, i) - fab*hmat(i, idir)*forcea(:)
1062 force_tmp(jatom, :, i) = force_tmp(jatom, :, i) - fab*hmat(i, idir)*forceb(:)
1063 END DO
1064 END DO
1065 CALL neighbor_list_iterator_release(nl_iterator)
1066 DEALLOCATE (basis_set_list)
1067 END IF
1068 END DO
1069
1070 ! make sure the total normalized polarization is within [-1:1]
1071 DO idir = 1, 3
1072 cqi(idir) = rlog(idir) + zlog(idir)
1073 IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
1074 IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
1075 ! now check for log branch
1076 IF (calculate_forces) THEN
1077 IF (abs(efield%polarisation(idir) - cqi(idir)) > pi) THEN
1078 di(idir) = (efield%polarisation(idir) - cqi(idir))/pi
1079 DO i = 1, 10
1080 cqi(idir) = cqi(idir) + sign(1.0_dp, di(idir))*twopi
1081 IF (abs(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
1082 END DO
1083 END IF
1084 END IF
1085 END DO
1086 DO idir = 1, 3
1087 qi(idir) = 0.0_dp
1088 ci(idir) = 0.0_dp
1089 DO i = 1, 3
1090 ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
1091 END DO
1092 END DO
1093
1094 ! update the references
1095 IF (calculate_forces) THEN
1096 ener_field = sum(ci)
1097 ! check for smoothness of energy surface
1098 IF (abs(efield%field_energy - ener_field) > pi*abs(sum(hmat))) THEN
1099 cpwarn("Large change of e-field energy detected. Correct for non-smooth energy surface")
1100 END IF
1101 efield%field_energy = ener_field
1102 efield%polarisation(:) = cqi(:)
1103 END IF
1104
1105 ! Energy
1106 ener_field = 0.0_dp
1107 DO i = 1, 3
1108 ener_field = ener_field + dfilter(i)*(fieldpol(i) - 2._dp*twopi*ci(i))**2
1109 END DO
1110 energy%efield = 0.25_dp*omega/twopi*ener_field
1111
1112 ! debugging output
1113 IF (para_env%is_source()) THEN
1114 iodeb = -1
1115 IF (iodeb > 0) THEN
1116 WRITE (iodeb, '(A,T61,F20.10)') " Polarisation Quantum: ", 2._dp*twopi*twopi*hmat(3, 3)
1117 WRITE (iodeb, '(A,T21,3F20.10)') " Polarisation: ", 2._dp*twopi*ci(1:3)
1118 WRITE (iodeb, '(A,T21,3F20.10)') " Displacement: ", fieldpol(1:3)
1119 WRITE (iodeb, '(A,T21,3F20.10)') " E-Field: ", ((fieldpol(i) - 2._dp*twopi*ci(i)), i=1, 3)
1120 WRITE (iodeb, '(A,T61,F20.10)') " Disp Free Energy:", energy%efield
1121 END IF
1122 END IF
1123
1124 IF (.NOT. just_energy) THEN
1125 DO i = 1, 3
1126 di(i) = -omega*(fieldpol(i) - 2._dp*twopi*ci(i))*dfilter(i)
1127 END DO
1128 ! Add the result to mo_derivativs
1129 DO ispin = 1, dft_control%nspins
1130 CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_coeff_tmp(ispin))
1131 DO idir = 1, 3
1132 CALL cp_fm_scale_and_add(1.0_dp, mo_coeff_tmp(ispin), di(idir), &
1133 mo_derivs_tmp(idir, ispin))
1134 END DO
1135 END DO
1136 DO ispin = 1, dft_control%nspins
1137 CALL copy_fm_to_dbcsr(mo_coeff_tmp(ispin), mo_derivs(ispin)%matrix)
1138 END DO
1139 END IF
1140
1141 IF (calculate_forces) THEN
1142 DO i = 1, 3
1143 DO ia = 1, natom
1144 CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
1145 iatom = atom_of_kind(ia)
1146 force(ikind)%efield(1:3, iatom) = force(ikind)%efield(1:3, iatom) + di(i)*force_tmp(ia, 1:3, i)
1147 END DO
1148 END DO
1149 END IF
1150
1151 DO ispin = 1, dft_control%nspins
1152 CALL cp_cfm_release(eigrmat(ispin))
1153 CALL cp_cfm_release(inv_mat(ispin))
1154 IF (mos(ispin)%use_mo_coeff_b) CALL cp_fm_release(mo_coeff_tmp(ispin))
1155 DO i = 1, 3
1156 CALL cp_fm_release(mo_derivs_tmp(i, ispin))
1157 END DO
1158 DO i = 1, SIZE(op_fm_set, 1)
1159 CALL cp_fm_release(opvec(i, ispin))
1160 CALL cp_fm_release(op_fm_set(i, ispin))
1161 CALL cp_fm_release(inv_work(i, ispin))
1162 END DO
1163 END DO
1164 DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat)
1165 DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp)
1166
1167 IF (calculate_forces) THEN
1168 DO ikind = 1, SIZE(atomic_kind_set)
1169 CALL para_env%sum(force(ikind)%efield)
1170 END DO
1171 DEALLOCATE (force_tmp)
1172 DEALLOCATE (cosab, sinab, work, dcosab, dsinab)
1173 DO i = 1, 3
1174 DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block)
1175 DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block)
1176 END DO
1177 CALL dbcsr_deallocate_matrix_set(tempmat)
1178 END IF
1179 CALL timestop(handle)
1180
1181 END SUBROUTINE qs_dispfield_derivatives
1182
1183! **************************************************************************************************
1184!> \brief ...
1185!> \param cos_block ...
1186!> \param sin_block ...
1187!> \param ncoa ...
1188!> \param nsgfa ...
1189!> \param sgfa ...
1190!> \param sphi_a ...
1191!> \param ldsa ...
1192!> \param ncob ...
1193!> \param nsgfb ...
1194!> \param sgfb ...
1195!> \param sphi_b ...
1196!> \param ldsb ...
1197!> \param cosab ...
1198!> \param sinab ...
1199!> \param ldab ...
1200!> \param work ...
1201!> \param ldwork ...
1202! **************************************************************************************************
1203 SUBROUTINE contract_all(cos_block, sin_block, &
1204 ncoa, nsgfa, sgfa, sphi_a, ldsa, &
1205 ncob, nsgfb, sgfb, sphi_b, ldsb, &
1206 cosab, sinab, ldab, work, ldwork)
1207
1208 REAL(dp), DIMENSION(:, :), POINTER :: cos_block, sin_block
1209 INTEGER, INTENT(IN) :: ncoa, nsgfa, sgfa
1210 REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_a
1211 INTEGER, INTENT(IN) :: ldsa, ncob, nsgfb, sgfb
1212 REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_b
1213 INTEGER, INTENT(IN) :: ldsb
1214 REAL(dp), DIMENSION(:, :), INTENT(IN) :: cosab, sinab
1215 INTEGER, INTENT(IN) :: ldab
1216 REAL(dp), DIMENSION(:, :) :: work
1217 INTEGER, INTENT(IN) :: ldwork
1218
1219! Calculate cosine
1220
1221 CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, cosab(1, 1), ldab, &
1222 sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
1223
1224 CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, &
1225 work(1, 1), ldwork, 1.0_dp, cos_block(sgfa, sgfb), SIZE(cos_block, 1))
1226
1227 ! Calculate sine
1228 CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, sinab(1, 1), ldab, &
1229 sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
1230
1231 CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, &
1232 work(1, 1), ldwork, 1.0_dp, sin_block(sgfa, sgfb), SIZE(sin_block, 1))
1233
1234 END SUBROUTINE contract_all
1235
1236END MODULE qs_efield_berry
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cell_types::pbc
Definition cell_types.F:104
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition ai_moments.F:17
ai_moments::cossin
subroutine, public cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
...
Definition ai_moments.F:155
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
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.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
block_p_types
collect pointers to a block of reals
Definition block_p_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_solve
subroutine, public cp_cfm_solve(matrix_a, general_a, determinant)
Solve the system of linear equations A*b=A_general using LU decomposition. Pay attention that both ma...
Definition cp_cfm_basic_linalg.F:731
cp_cfm_basic_linalg::cp_cfm_scale_and_add_fm
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
Definition cp_cfm_basic_linalg.F:351
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:202
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:548
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition cp_dbcsr_operations.F:693
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:170
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
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
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:142
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_efield_berry
Calculates the energy contribution and the mo_derivative of a static periodic electric field.
Definition qs_efield_berry.F:15
qs_efield_berry::qs_efield_berry_phase
subroutine, public qs_efield_berry_phase(qs_env, just_energy, calculate_forces)
...
Definition qs_efield_berry.F:107
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
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, matrix_vhxc, 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.
Definition qs_environment_types.F:539
qs_environment_types::set_qs_env
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.
Definition qs_environment_types.F:1128
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:401
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>.
Definition qs_moments.F:14
qs_moments::build_berry_moment_matrix
subroutine, public build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
...
Definition qs_moments.F:1045
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:550
qs_period_efield_types
type for berry phase efield matrices. At the moment only used for cosmat and sinmat
Definition qs_period_efield_types.F:16
qs_period_efield_types::set_efield_matrices
subroutine, public set_efield_matrices(efield, sinmat, cosmat, dipmat)
...
Definition qs_period_efield_types.F:81
qs_period_efield_types::init_efield_matrices
subroutine, public init_efield_matrices(efield)
...
Definition qs_period_efield_types.F:47
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:142
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
block_p_types::block_p_type
Definition block_p_types.F:23
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_control_types::dft_control_type
Definition cp_control_types.F:758
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:767
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_period_efield_types::efield_berry_type
Definition qs_period_efield_types.F:32
virial_types::virial_type
Definition virial_types.F:26