(git:b279b6b)
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-2024 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 ! **************************************************************************************************
15 MODULE 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_operations, ONLY: copy_dbcsr_to_fm,&
33  copy_fm_to_dbcsr,&
34  cp_dbcsr_plus_fm_fm_t,&
35  cp_dbcsr_sm_fm_multiply,&
36  dbcsr_deallocate_matrix_set
37  USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
38  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
39  cp_fm_struct_release,&
40  cp_fm_struct_type
41  USE cp_fm_types, ONLY: cp_fm_create,&
42  cp_fm_release,&
43  cp_fm_set_all,&
44  cp_fm_type
45  USE dbcsr_api, ONLY: dbcsr_copy,&
46  dbcsr_get_block_p,&
47  dbcsr_p_type,&
48  dbcsr_set,&
49  dbcsr_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 
96 CONTAINS
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  IF (s_mstruct_changed) CALL qs_efield_integrals(qs_env)
125  IF (dft_control%period_efield%displacement_field) THEN
126  CALL qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
127  ELSE
128  CALL qs_efield_derivatives(qs_env, just_energy, calculate_forces)
129  END IF
130  END IF
131 
132  CALL timestop(handle)
133 
134  END SUBROUTINE qs_efield_berry_phase
135 
136 ! **************************************************************************************************
137 !> \brief ...
138 !> \param qs_env ...
139 ! **************************************************************************************************
140  SUBROUTINE qs_efield_integrals(qs_env)
141 
142  TYPE(qs_environment_type), POINTER :: qs_env
143 
144  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_efield_integrals'
145 
146  INTEGER :: handle, i
147  REAL(dp), DIMENSION(3) :: kvec
148  TYPE(cell_type), POINTER :: cell
149  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: cosmat, matrix_s, sinmat
150  TYPE(dft_control_type), POINTER :: dft_control
151  TYPE(efield_berry_type), POINTER :: efield
152 
153  CALL timeset(routinen, handle)
154  cpassert(ASSOCIATED(qs_env))
155 
156  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
157  NULLIFY (matrix_s)
158  CALL get_qs_env(qs_env=qs_env, efield=efield, cell=cell, matrix_s=matrix_s)
159  CALL init_efield_matrices(efield)
160  ALLOCATE (cosmat(3), sinmat(3))
161  DO i = 1, 3
162  ALLOCATE (cosmat(i)%matrix, sinmat(i)%matrix)
163 
164  CALL dbcsr_copy(cosmat(i)%matrix, matrix_s(1)%matrix, 'COS MAT')
165  CALL dbcsr_copy(sinmat(i)%matrix, matrix_s(1)%matrix, 'SIN MAT')
166  CALL dbcsr_set(cosmat(i)%matrix, 0.0_dp)
167  CALL dbcsr_set(sinmat(i)%matrix, 0.0_dp)
168 
169  kvec(:) = twopi*cell%h_inv(i, :)
170  CALL build_berry_moment_matrix(qs_env, cosmat(i)%matrix, sinmat(i)%matrix, kvec)
171  END DO
172  CALL set_efield_matrices(efield=efield, cosmat=cosmat, sinmat=sinmat)
173  CALL set_qs_env(qs_env=qs_env, efield=efield)
174  CALL timestop(handle)
175 
176  END SUBROUTINE qs_efield_integrals
177 
178 ! **************************************************************************************************
179 !> \brief ...
180 !> \param qs_env ...
181 !> \param just_energy ...
182 !> \param calculate_forces ...
183 ! **************************************************************************************************
184  SUBROUTINE qs_efield_derivatives(qs_env, just_energy, calculate_forces)
185  TYPE(qs_environment_type), POINTER :: qs_env
186  LOGICAL, INTENT(IN) :: just_energy, calculate_forces
187 
188  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_efield_derivatives'
189 
190  COMPLEX(dp) :: zdet, zdeta, zi(3)
191  INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, j, &
192  jatom, jkind, jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, &
193  nseta, nsetb, sgfa, sgfb
194  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
195  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
196  npgfb, nsgfa, nsgfb
197  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
198  LOGICAL :: found, uniform, use_virial
199  REAL(dp) :: charge, ci(3), cqi(3), dab, dd, &
200  ener_field, f0, fab, fieldpol(3), &
201  focc, fpolvec(3), hmat(3, 3), occ, &
202  qi(3), ti(3)
203  REAL(dp), DIMENSION(3) :: forcea, forceb, kvec, ra, rab, rb, ria
204  REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work
205  REAL(dp), DIMENSION(:, :, :), POINTER :: dcosab, dsinab
206  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
207  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
208  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
209  TYPE(block_p_type), DIMENSION(3, 2) :: dcost, dsint
210  TYPE(cell_type), POINTER :: cell
211  TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat, inv_mat
212  TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
213  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff_tmp, mo_derivs_tmp
214  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: inv_work, op_fm_set, opvec
215  TYPE(cp_fm_type), POINTER :: mo_coeff
216  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, mo_derivs
217  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: tempmat
218  TYPE(dbcsr_type), POINTER :: cosmat, mo_coeff_b, sinmat
219  TYPE(dft_control_type), POINTER :: dft_control
220  TYPE(efield_berry_type), POINTER :: efield
221  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
222  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
223  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
224  TYPE(mp_para_env_type), POINTER :: para_env
225  TYPE(neighbor_list_iterator_p_type), &
226  DIMENSION(:), POINTER :: nl_iterator
227  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
228  POINTER :: sab_orb
229  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
230  TYPE(qs_energy_type), POINTER :: energy
231  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
232  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
233  TYPE(qs_kind_type), POINTER :: qs_kind
234  TYPE(virial_type), POINTER :: virial
235 
236  CALL timeset(routinen, handle)
237 
238  NULLIFY (dft_control, cell, particle_set)
239  CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
240  particle_set=particle_set, virial=virial)
241  NULLIFY (qs_kind_set, efield, para_env, sab_orb)
242  CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
243  efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
244 
245  ! calculate stress only if forces requested also
246  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
247  use_virial = use_virial .AND. calculate_forces
248  ! disable stress calculation
249  IF (use_virial) THEN
250  cpabort("Stress tensor for periodic E-field not implemented")
251  END IF
252 
253  fieldpol = dft_control%period_efield%polarisation
254  fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
255  fieldpol = -fieldpol*dft_control%period_efield%strength
256  hmat = cell%hmat(:, :)/twopi
257  DO idir = 1, 3
258  fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir)
259  END DO
260 
261  ! nuclear contribution
262  natom = SIZE(particle_set)
263  IF (calculate_forces) THEN
264  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
265  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
266  END IF
267  zi(:) = cmplx(1._dp, 0._dp, dp)
268  DO ia = 1, natom
269  CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
270  CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
271  ria = particle_set(ia)%r
272  ria = pbc(ria, cell)
273  DO idir = 1, 3
274  kvec(:) = twopi*cell%h_inv(idir, :)
275  dd = sum(kvec(:)*ria(:))
276  zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
277  zi(idir) = zi(idir)*zdeta
278  END DO
279  IF (calculate_forces) THEN
280  IF (para_env%mepos == 0) THEN
281  iatom = atom_of_kind(ia)
282  forcea(:) = fieldpol(:)*charge
283  force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + forcea(:)
284  END IF
285  END IF
286  IF (use_virial) THEN
287  IF (para_env%mepos == 0) &
288  CALL virial_pair_force(virial%pv_virial, 1.0_dp, forcea, ria)
289  END IF
290  END DO
291  qi = aimag(log(zi))
292 
293  ! check uniform occupation
294  NULLIFY (mos)
295  CALL get_qs_env(qs_env=qs_env, mos=mos)
296  DO ispin = 1, dft_control%nspins
297  CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
298  IF (.NOT. uniform) THEN
299  cpabort("Berry phase moments for non uniform MOs' occupation numbers not implemented")
300  END IF
301  END DO
302 
303  NULLIFY (mo_derivs)
304  CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs)
305  ! initialize all work matrices needed
306  ALLOCATE (op_fm_set(2, dft_control%nspins))
307  ALLOCATE (opvec(2, dft_control%nspins))
308  ALLOCATE (eigrmat(dft_control%nspins))
309  ALLOCATE (inv_mat(dft_control%nspins))
310  ALLOCATE (inv_work(2, dft_control%nspins))
311  ALLOCATE (mo_derivs_tmp(SIZE(mo_derivs)))
312  ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs)))
313 
314  ! Allocate temp matrices for the wavefunction derivatives
315  DO ispin = 1, dft_control%nspins
316  NULLIFY (tmp_fm_struct, mo_coeff)
317  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
318  CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
319  ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
320  CALL cp_fm_create(mo_derivs_tmp(ispin), mo_coeff%matrix_struct)
321  CALL cp_fm_create(mo_coeff_tmp(ispin), mo_coeff%matrix_struct)
322  CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_derivs_tmp(ispin))
323  DO i = 1, SIZE(op_fm_set, 1)
324  CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
325  CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
326  CALL cp_fm_create(inv_work(i, ispin), op_fm_set(i, ispin)%matrix_struct)
327  END DO
328  CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
329  CALL cp_cfm_create(inv_mat(ispin), op_fm_set(1, ispin)%matrix_struct)
330  CALL cp_fm_struct_release(tmp_fm_struct)
331  END DO
332  ! temp matrices for force calculation
333  IF (calculate_forces) THEN
334  NULLIFY (matrix_s)
335  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
336  ALLOCATE (tempmat(2, dft_control%nspins))
337  DO ispin = 1, dft_control%nspins
338  ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix)
339  CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
340  CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
341  CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
342  CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
343  END DO
344  ! integration
345  CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab)
346  ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab))
347  ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3))
348  lsab = max(ldab, lsab)
349  DO i = 1, 3
350  ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab))
351  ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab))
352  END DO
353  END IF
354 
355  !Start the MO derivative calculation
356  !loop over all cell vectors
357  DO idir = 1, 3
358  ci(idir) = 0.0_dp
359  zi(idir) = z_zero
360  IF (abs(fpolvec(idir)) .GT. 1.0e-12_dp) THEN
361  cosmat => efield%cosmat(idir)%matrix
362  sinmat => efield%sinmat(idir)%matrix
363  !evaluate the expression needed for the derivative (S_berry * C and [C^T S_berry C]^-1)
364  !first step S_berry * C and C^T S_berry C
365  DO ispin = 1, dft_control%nspins ! spin
366  IF (mos(ispin)%use_mo_coeff_b) THEN
367  CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo)
368  CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin))
369  ELSE
370  CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
371  mo_coeff_tmp(ispin) = mo_coeff
372  END IF
373  CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin), opvec(1, ispin), ncol=nmo)
374  CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(1, ispin), 0.0_dp, &
375  op_fm_set(1, ispin))
376  CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin), opvec(2, ispin), ncol=nmo)
377  CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(2, ispin), 0.0_dp, &
378  op_fm_set(2, ispin))
379  END DO
380  !second step invert C^T S_berry C
381  zdet = z_one
382  DO ispin = 1, dft_control%nspins
383  CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin), z_one, op_fm_set(1, ispin))
384  CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin), -gaussi, op_fm_set(2, ispin))
385  CALL cp_cfm_set_all(inv_mat(ispin), z_zero, z_one)
386  CALL cp_cfm_solve(eigrmat(ispin), inv_mat(ispin), zdeta)
387  zdet = zdet*zdeta
388  END DO
389  zi(idir) = zdet**occ
390  ci(idir) = aimag(log(zdet**occ))
391 
392  IF (.NOT. just_energy) THEN
393  !compute the orbital derivative
394  focc = fpolvec(idir)
395  DO ispin = 1, dft_control%nspins
396  inv_work(1, ispin)%local_data(:, :) = real(inv_mat(ispin)%local_data(:, :), dp)
397  inv_work(2, ispin)%local_data(:, :) = aimag(inv_mat(ispin)%local_data(:, :))
398  CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
399  CALL parallel_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin), inv_work(2, ispin), &
400  1.0_dp, mo_derivs_tmp(ispin))
401  CALL parallel_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin), inv_work(1, ispin), &
402  1.0_dp, mo_derivs_tmp(ispin))
403  END DO
404  END IF
405 
406  !compute nuclear forces
407  IF (calculate_forces) THEN
408  nkind = SIZE(qs_kind_set)
409  natom = SIZE(particle_set)
410  kvec(:) = twopi*cell%h_inv(idir, :)
411 
412  ! calculate: C [C^T S_berry C]^(-1) C^T
413  ! Store this matrix in DBCSR form (only S overlap blocks)
414  DO ispin = 1, dft_control%nspins
415  CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
416  CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
417  CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
418  CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(1, ispin), 0.0_dp, &
419  opvec(1, ispin))
420  CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(2, ispin), 0.0_dp, &
421  opvec(2, ispin))
422  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, &
423  matrix_v=opvec(1, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
424  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, &
425  matrix_v=opvec(2, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
426  END DO
427 
428  ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db)
429  ALLOCATE (basis_set_list(nkind))
430  DO ikind = 1, nkind
431  qs_kind => qs_kind_set(ikind)
432  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
433  IF (ASSOCIATED(basis_set_a)) THEN
434  basis_set_list(ikind)%gto_basis_set => basis_set_a
435  ELSE
436  NULLIFY (basis_set_list(ikind)%gto_basis_set)
437  END IF
438  END DO
439  !
440  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
441  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
442  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
443  iatom=iatom, jatom=jatom, r=rab)
444  basis_set_a => basis_set_list(ikind)%gto_basis_set
445  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
446  basis_set_b => basis_set_list(jkind)%gto_basis_set
447  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
448  ! basis ikind
449  first_sgfa => basis_set_a%first_sgf
450  la_max => basis_set_a%lmax
451  la_min => basis_set_a%lmin
452  npgfa => basis_set_a%npgf
453  nseta = basis_set_a%nset
454  nsgfa => basis_set_a%nsgf_set
455  rpgfa => basis_set_a%pgf_radius
456  set_radius_a => basis_set_a%set_radius
457  sphi_a => basis_set_a%sphi
458  zeta => basis_set_a%zet
459  ! basis jkind
460  first_sgfb => basis_set_b%first_sgf
461  lb_max => basis_set_b%lmax
462  lb_min => basis_set_b%lmin
463  npgfb => basis_set_b%npgf
464  nsetb = basis_set_b%nset
465  nsgfb => basis_set_b%nsgf_set
466  rpgfb => basis_set_b%pgf_radius
467  set_radius_b => basis_set_b%set_radius
468  sphi_b => basis_set_b%sphi
469  zetb => basis_set_b%zet
470 
471  atom_a = atom_of_kind(iatom)
472  atom_b = atom_of_kind(jatom)
473 
474  ldsa = SIZE(sphi_a, 1)
475  ldsb = SIZE(sphi_b, 1)
476  ra(:) = pbc(particle_set(iatom)%r(:), cell)
477  rb(:) = ra + rab
478  dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
479 
480  IF (iatom <= jatom) THEN
481  irow = iatom
482  icol = jatom
483  ELSE
484  irow = jatom
485  icol = iatom
486  END IF
487 
488  IF (iatom == jatom .AND. dab < 1.e-10_dp) THEN
489  fab = 1.0_dp*occ
490  ELSE
491  fab = 2.0_dp*occ
492  END IF
493 
494  DO i = 1, 3
495  dcost(i, 1)%block = 0.0_dp
496  dsint(i, 1)%block = 0.0_dp
497  dcost(i, 2)%block = 0.0_dp
498  dsint(i, 2)%block = 0.0_dp
499  END DO
500 
501  DO iset = 1, nseta
502  ncoa = npgfa(iset)*ncoset(la_max(iset))
503  sgfa = first_sgfa(1, iset)
504  DO jset = 1, nsetb
505  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
506  ncob = npgfb(jset)*ncoset(lb_max(jset))
507  sgfb = first_sgfb(1, jset)
508  ! Calculate the primitive integrals (da|b)
509  CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
510  lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
511  ra, rb, kvec, cosab, sinab, dcosab, dsinab)
512  DO i = 1, 3
513  CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, &
514  ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
515  ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
516  dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
517  END DO
518  ! Calculate the primitive integrals (a|db)
519  CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
520  la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
521  rb, ra, kvec, cosab, sinab, dcosab, dsinab)
522  DO i = 1, 3
523  dcosab(1:ncoa, 1:ncob, i) = transpose(dcosab(1:ncob, 1:ncoa, i))
524  dsinab(1:ncoa, 1:ncob, i) = transpose(dsinab(1:ncob, 1:ncoa, i))
525  CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, &
526  ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
527  ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
528  dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
529  END DO
530  END DO
531  END DO
532  forcea = 0.0_dp
533  forceb = 0.0_dp
534  DO ispin = 1, dft_control%nspins
535  NULLIFY (rblock, iblock)
536  CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, &
537  row=irow, col=icol, block=rblock, found=found)
538  cpassert(found)
539  CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, &
540  row=irow, col=icol, block=iblock, found=found)
541  cpassert(found)
542  n1 = SIZE(rblock, 1)
543  n2 = SIZE(rblock, 2)
544  cpassert(SIZE(iblock, 1) == n1)
545  cpassert(SIZE(iblock, 2) == n2)
546  cpassert(lsab >= n1)
547  cpassert(lsab >= n2)
548  IF (iatom <= jatom) THEN
549  DO i = 1, 3
550  forcea(i) = forcea(i) + sum(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) &
551  - sum(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2))
552  forceb(i) = forceb(i) + sum(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) &
553  - sum(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2))
554  END DO
555  ELSE
556  DO i = 1, 3
557  forcea(i) = forcea(i) + sum(transpose(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) &
558  - sum(transpose(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1))
559  forceb(i) = forceb(i) + sum(transpose(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) &
560  - sum(transpose(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1))
561  END DO
562  END IF
563  END DO
564  force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fab*fpolvec(idir)*forcea(1:3)
565  force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fab*fpolvec(idir)*forceb(1:3)
566  IF (use_virial) THEN
567  f0 = -fab*fpolvec(idir)
568  CALL virial_pair_force(virial%pv_virial, f0, forcea, ra)
569  CALL virial_pair_force(virial%pv_virial, f0, forceb, rb)
570  END IF
571 
572  END DO
573  CALL neighbor_list_iterator_release(nl_iterator)
574  DEALLOCATE (basis_set_list)
575 
576  END IF
577  END IF
578  END DO
579 
580  ! Energy
581  ener_field = 0.0_dp
582  ti = 0.0_dp
583  DO idir = 1, 3
584  ! make sure the total normalized polarization is within [-1:1]
585  cqi(idir) = qi(idir) + ci(idir)
586  IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
587  IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
588  ! now check for log branch
589  IF (abs(efield%polarisation(idir) - cqi(idir)) > pi) THEN
590  ti(idir) = (efield%polarisation(idir) - cqi(idir))/pi
591  DO i = 1, 10
592  cqi(idir) = cqi(idir) + sign(1.0_dp, ti(idir))*twopi
593  IF (abs(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
594  END DO
595  END IF
596  ener_field = ener_field + fpolvec(idir)*cqi(idir)
597  END DO
598 
599  ! update the references
600  IF (calculate_forces) THEN
601  ! check for smoothness of energy surface
602  IF (abs(efield%field_energy - ener_field) > pi*abs(sum(fpolvec))) THEN
603  cpwarn("Large change of e-field energy detected. Correct for non-smooth energy surface")
604  END IF
605  efield%field_energy = ener_field
606  efield%polarisation(:) = cqi(:)
607  END IF
608  energy%efield = ener_field
609 
610  IF (.NOT. just_energy) THEN
611  ! Add the result to mo_derivativs
612  DO ispin = 1, dft_control%nspins
613  CALL copy_fm_to_dbcsr(mo_derivs_tmp(ispin), mo_derivs(ispin)%matrix)
614  END DO
615  IF (use_virial) THEN
616  ti = 0.0_dp
617  DO i = 1, 3
618  DO j = 1, 3
619  ti(j) = ti(j) + hmat(j, i)*cqi(i)
620  END DO
621  END DO
622  DO i = 1, 3
623  DO j = 1, 3
624  virial%pv_virial(i, j) = virial%pv_virial(i, j) - fieldpol(i)*ti(j)
625  END DO
626  END DO
627  END IF
628  END IF
629 
630  DO ispin = 1, dft_control%nspins
631  CALL cp_cfm_release(eigrmat(ispin))
632  CALL cp_cfm_release(inv_mat(ispin))
633  CALL cp_fm_release(mo_derivs_tmp(ispin))
634  IF (mos(ispin)%use_mo_coeff_b) CALL cp_fm_release(mo_coeff_tmp(ispin))
635  DO i = 1, SIZE(op_fm_set, 1)
636  CALL cp_fm_release(opvec(i, ispin))
637  CALL cp_fm_release(op_fm_set(i, ispin))
638  CALL cp_fm_release(inv_work(i, ispin))
639  END DO
640  END DO
641  DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat)
642  DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp)
643 
644  IF (calculate_forces) THEN
645  DO ikind = 1, SIZE(atomic_kind_set)
646  CALL para_env%sum(force(ikind)%efield)
647  END DO
648  DEALLOCATE (cosab, sinab, work, dcosab, dsinab)
649  DO i = 1, 3
650  DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block)
651  DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block)
652  END DO
653  CALL dbcsr_deallocate_matrix_set(tempmat)
654  END IF
655  CALL timestop(handle)
656 
657  END SUBROUTINE qs_efield_derivatives
658 
659 ! **************************************************************************************************
660 !> \brief ...
661 !> \param qs_env ...
662 !> \param just_energy ...
663 !> \param calculate_forces ...
664 ! **************************************************************************************************
665  SUBROUTINE qs_dispfield_derivatives(qs_env, just_energy, calculate_forces)
666  TYPE(qs_environment_type), POINTER :: qs_env
667  LOGICAL, INTENT(IN) :: just_energy, calculate_forces
668 
669  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_dispfield_derivatives'
670 
671  COMPLEX(dp) :: zdet, zdeta, zi(3)
672  INTEGER :: handle, i, ia, iatom, icol, idir, ikind, iodeb, irow, iset, ispin, jatom, jkind, &
673  jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, nseta, nsetb, &
674  sgfa, sgfb
675  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
676  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
677  npgfb, nsgfa, nsgfb
678  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
679  LOGICAL :: found, uniform, use_virial
680  REAL(dp) :: charge, ci(3), cqi(3), dab, dd, di(3), &
681  ener_field, fab, fieldpol(3), focc, &
682  hmat(3, 3), occ, omega, qi(3), &
683  rlog(3), zlog(3)
684  REAL(dp), DIMENSION(3) :: dfilter, forcea, forceb, kvec, ra, rab, &
685  rb, ria
686  REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work
687  REAL(dp), DIMENSION(:, :, :), POINTER :: dcosab, dsinab, force_tmp
688  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
689  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
690  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
691  TYPE(block_p_type), DIMENSION(3, 2) :: dcost, dsint
692  TYPE(cell_type), POINTER :: cell
693  TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat, inv_mat
694  TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
695  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff_tmp
696  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: inv_work, mo_derivs_tmp, op_fm_set, opvec
697  TYPE(cp_fm_type), POINTER :: mo_coeff
698  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, mo_derivs
699  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: tempmat
700  TYPE(dbcsr_type), POINTER :: cosmat, mo_coeff_b, sinmat
701  TYPE(dft_control_type), POINTER :: dft_control
702  TYPE(efield_berry_type), POINTER :: efield
703  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
704  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
705  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
706  TYPE(mp_para_env_type), POINTER :: para_env
707  TYPE(neighbor_list_iterator_p_type), &
708  DIMENSION(:), POINTER :: nl_iterator
709  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
710  POINTER :: sab_orb
711  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
712  TYPE(qs_energy_type), POINTER :: energy
713  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
714  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
715  TYPE(qs_kind_type), POINTER :: qs_kind
716  TYPE(virial_type), POINTER :: virial
717 
718  CALL timeset(routinen, handle)
719 
720  NULLIFY (dft_control, cell, particle_set)
721  CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, &
722  particle_set=particle_set, virial=virial)
723  NULLIFY (qs_kind_set, efield, para_env, sab_orb)
724  CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
725  efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb)
726 
727  ! calculate stress only if forces requested also
728  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
729  use_virial = use_virial .AND. calculate_forces
730  ! disable stress calculation
731  IF (use_virial) THEN
732  cpabort("Stress tensor for periodic D-field not implemented")
733  END IF
734 
735  dfilter(1:3) = dft_control%period_efield%d_filter(1:3)
736 
737  fieldpol = dft_control%period_efield%polarisation
738  fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
739  fieldpol = fieldpol*dft_control%period_efield%strength
740 
741  omega = cell%deth
742  hmat = cell%hmat(:, :)/(twopi*omega)
743 
744  ! nuclear contribution to polarization
745  natom = SIZE(particle_set)
746  IF (calculate_forces) THEN
747  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
748  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
749  ALLOCATE (force_tmp(natom, 3, 3))
750  force_tmp = 0.0_dp
751  END IF
752  zi(:) = cmplx(1._dp, 0._dp, dp)
753  DO ia = 1, natom
754  CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
755  CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
756  ria = particle_set(ia)%r
757  ria = pbc(ria, cell)
758  DO idir = 1, 3
759  kvec(:) = twopi*cell%h_inv(idir, :)
760  dd = sum(kvec(:)*ria(:))
761  zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
762  zi(idir) = zi(idir)*zdeta
763  END DO
764  IF (calculate_forces) THEN
765  IF (para_env%mepos == 0) THEN
766  DO i = 1, 3
767  force_tmp(ia, i, i) = force_tmp(ia, i, i) + charge/omega
768  END DO
769  END IF
770  END IF
771  END DO
772  rlog = aimag(log(zi))
773 
774  ! check uniform occupation
775  NULLIFY (mos)
776  CALL get_qs_env(qs_env=qs_env, mos=mos)
777  DO ispin = 1, dft_control%nspins
778  CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
779  IF (.NOT. uniform) THEN
780  cpabort("Berry phase moments for non uniform MO occupation numbers not implemented")
781  END IF
782  END DO
783 
784  ! initialize all work matrices needed
785  NULLIFY (mo_derivs)
786  CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs)
787  ALLOCATE (op_fm_set(2, dft_control%nspins))
788  ALLOCATE (opvec(2, dft_control%nspins))
789  ALLOCATE (eigrmat(dft_control%nspins))
790  ALLOCATE (inv_mat(dft_control%nspins))
791  ALLOCATE (inv_work(2, dft_control%nspins))
792  ALLOCATE (mo_derivs_tmp(3, SIZE(mo_derivs)))
793  ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs)))
794 
795  ! Allocate temp matrices for the wavefunction derivatives
796  DO ispin = 1, dft_control%nspins
797  NULLIFY (tmp_fm_struct, mo_coeff)
798  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
799  CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
800  ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
801  CALL cp_fm_create(mo_coeff_tmp(ispin), mo_coeff%matrix_struct)
802  DO i = 1, 3
803  CALL cp_fm_create(mo_derivs_tmp(i, ispin), mo_coeff%matrix_struct)
804  CALL cp_fm_set_all(matrix=mo_derivs_tmp(i, ispin), alpha=0.0_dp)
805  END DO
806  DO i = 1, SIZE(op_fm_set, 1)
807  CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
808  CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
809  CALL cp_fm_create(inv_work(i, ispin), op_fm_set(i, ispin)%matrix_struct)
810  END DO
811  CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
812  CALL cp_cfm_create(inv_mat(ispin), op_fm_set(1, ispin)%matrix_struct)
813  CALL cp_fm_struct_release(tmp_fm_struct)
814  END DO
815  ! temp matrices for force calculation
816  IF (calculate_forces) THEN
817  NULLIFY (matrix_s)
818  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
819  ALLOCATE (tempmat(2, dft_control%nspins))
820  DO ispin = 1, dft_control%nspins
821  ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix)
822  CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
823  CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT')
824  CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
825  CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
826  END DO
827  ! integration
828  CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab)
829  ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab))
830  ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3))
831  lsab = max(lsab, ldab)
832  DO i = 1, 3
833  ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab))
834  ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab))
835  END DO
836  END IF
837 
838  !Start the MO derivative calculation
839  !loop over all cell vectors
840  DO idir = 1, 3
841  zi(idir) = z_zero
842  cosmat => efield%cosmat(idir)%matrix
843  sinmat => efield%sinmat(idir)%matrix
844  !evaluate the expression needed for the derivative (S_berry * C and [C^T S_berry C]^-1)
845  !first step S_berry * C and C^T S_berry C
846  DO ispin = 1, dft_control%nspins ! spin
847  IF (mos(ispin)%use_mo_coeff_b) THEN
848  CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo)
849  CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin))
850  ELSE
851  CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
852  mo_coeff_tmp(ispin) = mo_coeff
853  END IF
854  CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin), opvec(1, ispin), ncol=nmo)
855  CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(1, ispin), 0.0_dp, &
856  op_fm_set(1, ispin))
857  CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin), opvec(2, ispin), ncol=nmo)
858  CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin), opvec(2, ispin), 0.0_dp, &
859  op_fm_set(2, ispin))
860  END DO
861  !second step invert C^T S_berry C
862  zdet = z_one
863  DO ispin = 1, dft_control%nspins
864  CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin), z_one, op_fm_set(1, ispin))
865  CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin), -gaussi, op_fm_set(2, ispin))
866  CALL cp_cfm_set_all(inv_mat(ispin), z_zero, z_one)
867  CALL cp_cfm_solve(eigrmat(ispin), inv_mat(ispin), zdeta)
868  zdet = zdet*zdeta
869  END DO
870  zi(idir) = zdet**occ
871  zlog(idir) = aimag(log(zi(idir)))
872 
873  IF (.NOT. just_energy) THEN
874  !compute the orbital derivative
875  DO ispin = 1, dft_control%nspins
876  inv_work(1, ispin)%local_data(:, :) = real(inv_mat(ispin)%local_data(:, :), dp)
877  inv_work(2, ispin)%local_data(:, :) = aimag(inv_mat(ispin)%local_data(:, :))
878  CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
879  DO i = 1, 3
880  focc = hmat(idir, i)
881  CALL parallel_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin), inv_work(2, ispin), &
882  1.0_dp, mo_derivs_tmp(idir, ispin))
883  CALL parallel_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin), inv_work(1, ispin), &
884  1.0_dp, mo_derivs_tmp(idir, ispin))
885  END DO
886  END DO
887  END IF
888 
889  !compute nuclear forces
890  IF (calculate_forces) THEN
891  nkind = SIZE(qs_kind_set)
892  natom = SIZE(particle_set)
893  kvec(:) = twopi*cell%h_inv(idir, :)
894 
895  ! calculate: C [C^T S_berry C]^(-1) C^T
896  ! Store this matrix in DBCSR form (only S overlap blocks)
897  DO ispin = 1, dft_control%nspins
898  CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp)
899  CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp)
900  CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
901  CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(1, ispin), 0.0_dp, &
902  opvec(1, ispin))
903  CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin), inv_work(2, ispin), 0.0_dp, &
904  opvec(2, ispin))
905  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, &
906  matrix_v=opvec(1, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
907  CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, &
908  matrix_v=opvec(2, ispin), matrix_g=mo_coeff_tmp(ispin), ncol=nmo)
909  END DO
910 
911  ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db)
912  ALLOCATE (basis_set_list(nkind))
913  DO ikind = 1, nkind
914  qs_kind => qs_kind_set(ikind)
915  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
916  IF (ASSOCIATED(basis_set_a)) THEN
917  basis_set_list(ikind)%gto_basis_set => basis_set_a
918  ELSE
919  NULLIFY (basis_set_list(ikind)%gto_basis_set)
920  END IF
921  END DO
922  !
923  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
924  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
925  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
926  iatom=iatom, jatom=jatom, r=rab)
927  basis_set_a => basis_set_list(ikind)%gto_basis_set
928  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
929  basis_set_b => basis_set_list(jkind)%gto_basis_set
930  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
931  ! basis ikind
932  first_sgfa => basis_set_a%first_sgf
933  la_max => basis_set_a%lmax
934  la_min => basis_set_a%lmin
935  npgfa => basis_set_a%npgf
936  nseta = basis_set_a%nset
937  nsgfa => basis_set_a%nsgf_set
938  rpgfa => basis_set_a%pgf_radius
939  set_radius_a => basis_set_a%set_radius
940  sphi_a => basis_set_a%sphi
941  zeta => basis_set_a%zet
942  ! basis jkind
943  first_sgfb => basis_set_b%first_sgf
944  lb_max => basis_set_b%lmax
945  lb_min => basis_set_b%lmin
946  npgfb => basis_set_b%npgf
947  nsetb = basis_set_b%nset
948  nsgfb => basis_set_b%nsgf_set
949  rpgfb => basis_set_b%pgf_radius
950  set_radius_b => basis_set_b%set_radius
951  sphi_b => basis_set_b%sphi
952  zetb => basis_set_b%zet
953 
954  ldsa = SIZE(sphi_a, 1)
955  ldsb = SIZE(sphi_b, 1)
956  ra(:) = pbc(particle_set(iatom)%r(:), cell)
957  rb(:) = ra + rab
958  dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
959 
960  IF (iatom <= jatom) THEN
961  irow = iatom
962  icol = jatom
963  ELSE
964  irow = jatom
965  icol = iatom
966  END IF
967 
968  IF (iatom == jatom .AND. dab < 1.e-10_dp) THEN
969  fab = 1.0_dp*occ
970  ELSE
971  fab = 2.0_dp*occ
972  END IF
973 
974  DO i = 1, 3
975  dcost(i, 1)%block = 0.0_dp
976  dsint(i, 1)%block = 0.0_dp
977  dcost(i, 2)%block = 0.0_dp
978  dsint(i, 2)%block = 0.0_dp
979  END DO
980 
981  DO iset = 1, nseta
982  ncoa = npgfa(iset)*ncoset(la_max(iset))
983  sgfa = first_sgfa(1, iset)
984  DO jset = 1, nsetb
985  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
986  ncob = npgfb(jset)*ncoset(lb_max(jset))
987  sgfb = first_sgfb(1, jset)
988  ! Calculate the primitive integrals (da|b)
989  CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
990  lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
991  ra, rb, kvec, cosab, sinab, dcosab, dsinab)
992  DO i = 1, 3
993  CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, &
994  ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
995  ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
996  dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
997  END DO
998  ! Calculate the primitive integrals (a|db)
999  CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1000  la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1001  rb, ra, kvec, cosab, sinab, dcosab, dsinab)
1002  DO i = 1, 3
1003  dcosab(1:ncoa, 1:ncob, i) = transpose(dcosab(1:ncob, 1:ncoa, i))
1004  dsinab(1:ncoa, 1:ncob, i) = transpose(dsinab(1:ncob, 1:ncoa, i))
1005  CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, &
1006  ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1007  ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1008  dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab)
1009  END DO
1010  END DO
1011  END DO
1012  forcea = 0.0_dp
1013  forceb = 0.0_dp
1014  DO ispin = 1, dft_control%nspins
1015  NULLIFY (rblock, iblock)
1016  CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, &
1017  row=irow, col=icol, block=rblock, found=found)
1018  cpassert(found)
1019  CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, &
1020  row=irow, col=icol, block=iblock, found=found)
1021  cpassert(found)
1022  n1 = SIZE(rblock, 1)
1023  n2 = SIZE(rblock, 2)
1024  cpassert(SIZE(iblock, 1) == n1)
1025  cpassert(SIZE(iblock, 2) == n2)
1026  cpassert(lsab >= n1)
1027  cpassert(lsab >= n2)
1028  IF (iatom <= jatom) THEN
1029  DO i = 1, 3
1030  forcea(i) = forcea(i) + sum(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) &
1031  - sum(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2))
1032  forceb(i) = forceb(i) + sum(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) &
1033  - sum(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2))
1034  END DO
1035  ELSE
1036  DO i = 1, 3
1037  forcea(i) = forcea(i) + sum(transpose(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) &
1038  - sum(transpose(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1))
1039  forceb(i) = forceb(i) + sum(transpose(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) &
1040  - sum(transpose(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1))
1041  END DO
1042  END IF
1043  END DO
1044  DO i = 1, 3
1045  force_tmp(iatom, :, i) = force_tmp(iatom, :, i) - fab*hmat(i, idir)*forcea(:)
1046  force_tmp(jatom, :, i) = force_tmp(jatom, :, i) - fab*hmat(i, idir)*forceb(:)
1047  END DO
1048  END DO
1049  CALL neighbor_list_iterator_release(nl_iterator)
1050  DEALLOCATE (basis_set_list)
1051  END IF
1052  END DO
1053 
1054  ! make sure the total normalized polarization is within [-1:1]
1055  DO idir = 1, 3
1056  cqi(idir) = rlog(idir) + zlog(idir)
1057  IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi
1058  IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi
1059  ! now check for log branch
1060  IF (calculate_forces) THEN
1061  IF (abs(efield%polarisation(idir) - cqi(idir)) > pi) THEN
1062  di(idir) = (efield%polarisation(idir) - cqi(idir))/pi
1063  DO i = 1, 10
1064  cqi(idir) = cqi(idir) + sign(1.0_dp, di(idir))*twopi
1065  IF (abs(efield%polarisation(idir) - cqi(idir)) < pi) EXIT
1066  END DO
1067  END IF
1068  END IF
1069  END DO
1070  DO idir = 1, 3
1071  qi(idir) = 0.0_dp
1072  ci(idir) = 0.0_dp
1073  DO i = 1, 3
1074  ci(idir) = ci(idir) + hmat(idir, i)*cqi(i)
1075  END DO
1076  END DO
1077 
1078  ! update the references
1079  IF (calculate_forces) THEN
1080  ener_field = sum(ci)
1081  ! check for smoothness of energy surface
1082  IF (abs(efield%field_energy - ener_field) > pi*abs(sum(hmat))) THEN
1083  cpwarn("Large change of e-field energy detected. Correct for non-smooth energy surface")
1084  END IF
1085  efield%field_energy = ener_field
1086  efield%polarisation(:) = cqi(:)
1087  END IF
1088 
1089  ! Energy
1090  ener_field = 0.0_dp
1091  DO i = 1, 3
1092  ener_field = ener_field + dfilter(i)*(fieldpol(i) - 2._dp*twopi*ci(i))**2
1093  END DO
1094  energy%efield = 0.25_dp*omega/twopi*ener_field
1095 
1096  ! debugging output
1097  IF (para_env%is_source()) THEN
1098  iodeb = -1
1099  IF (iodeb > 0) THEN
1100  WRITE (iodeb, '(A,T61,F20.10)') " Polarisation Quantum: ", 2._dp*twopi*twopi*hmat(3, 3)
1101  WRITE (iodeb, '(A,T21,3F20.10)') " Polarisation: ", 2._dp*twopi*ci(1:3)
1102  WRITE (iodeb, '(A,T21,3F20.10)') " Displacement: ", fieldpol(1:3)
1103  WRITE (iodeb, '(A,T21,3F20.10)') " E-Field: ", ((fieldpol(i) - 2._dp*twopi*ci(i)), i=1, 3)
1104  WRITE (iodeb, '(A,T61,F20.10)') " Disp Free Energy:", energy%efield
1105  END IF
1106  END IF
1107 
1108  IF (.NOT. just_energy) THEN
1109  DO i = 1, 3
1110  di(i) = -omega*(fieldpol(i) - 2._dp*twopi*ci(i))*dfilter(i)
1111  END DO
1112  ! Add the result to mo_derivativs
1113  DO ispin = 1, dft_control%nspins
1114  CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_coeff_tmp(ispin))
1115  DO idir = 1, 3
1116  CALL cp_fm_scale_and_add(1.0_dp, mo_coeff_tmp(ispin), di(idir), &
1117  mo_derivs_tmp(idir, ispin))
1118  END DO
1119  END DO
1120  DO ispin = 1, dft_control%nspins
1121  CALL copy_fm_to_dbcsr(mo_coeff_tmp(ispin), mo_derivs(ispin)%matrix)
1122  END DO
1123  END IF
1124 
1125  IF (calculate_forces) THEN
1126  DO i = 1, 3
1127  DO ia = 1, natom
1128  CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
1129  iatom = atom_of_kind(ia)
1130  force(ikind)%efield(1:3, iatom) = force(ikind)%efield(1:3, iatom) + di(i)*force_tmp(ia, 1:3, i)
1131  END DO
1132  END DO
1133  END IF
1134 
1135  DO ispin = 1, dft_control%nspins
1136  CALL cp_cfm_release(eigrmat(ispin))
1137  CALL cp_cfm_release(inv_mat(ispin))
1138  IF (mos(ispin)%use_mo_coeff_b) CALL cp_fm_release(mo_coeff_tmp(ispin))
1139  DO i = 1, 3
1140  CALL cp_fm_release(mo_derivs_tmp(i, ispin))
1141  END DO
1142  DO i = 1, SIZE(op_fm_set, 1)
1143  CALL cp_fm_release(opvec(i, ispin))
1144  CALL cp_fm_release(op_fm_set(i, ispin))
1145  CALL cp_fm_release(inv_work(i, ispin))
1146  END DO
1147  END DO
1148  DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat)
1149  DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp)
1150 
1151  IF (calculate_forces) THEN
1152  DO ikind = 1, SIZE(atomic_kind_set)
1153  CALL para_env%sum(force(ikind)%efield)
1154  END DO
1155  DEALLOCATE (force_tmp)
1156  DEALLOCATE (cosab, sinab, work, dcosab, dsinab)
1157  DO i = 1, 3
1158  DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block)
1159  DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block)
1160  END DO
1161  CALL dbcsr_deallocate_matrix_set(tempmat)
1162  END IF
1163  CALL timestop(handle)
1164 
1165  END SUBROUTINE qs_dispfield_derivatives
1166 
1167 ! **************************************************************************************************
1168 !> \brief ...
1169 !> \param cos_block ...
1170 !> \param sin_block ...
1171 !> \param ncoa ...
1172 !> \param nsgfa ...
1173 !> \param sgfa ...
1174 !> \param sphi_a ...
1175 !> \param ldsa ...
1176 !> \param ncob ...
1177 !> \param nsgfb ...
1178 !> \param sgfb ...
1179 !> \param sphi_b ...
1180 !> \param ldsb ...
1181 !> \param cosab ...
1182 !> \param sinab ...
1183 !> \param ldab ...
1184 !> \param work ...
1185 !> \param ldwork ...
1186 ! **************************************************************************************************
1187  SUBROUTINE contract_all(cos_block, sin_block, &
1188  ncoa, nsgfa, sgfa, sphi_a, ldsa, &
1189  ncob, nsgfb, sgfb, sphi_b, ldsb, &
1190  cosab, sinab, ldab, work, ldwork)
1191 
1192  REAL(dp), DIMENSION(:, :), POINTER :: cos_block, sin_block
1193  INTEGER, INTENT(IN) :: ncoa, nsgfa, sgfa
1194  REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_a
1195  INTEGER, INTENT(IN) :: ldsa, ncob, nsgfb, sgfb
1196  REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_b
1197  INTEGER, INTENT(IN) :: ldsb
1198  REAL(dp), DIMENSION(:, :), INTENT(IN) :: cosab, sinab
1199  INTEGER, INTENT(IN) :: ldab
1200  REAL(dp), DIMENSION(:, :) :: work
1201  INTEGER, INTENT(IN) :: ldwork
1202 
1203 ! Calculate cosine
1204 
1205  CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, cosab(1, 1), ldab, &
1206  sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
1207 
1208  CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, &
1209  work(1, 1), ldwork, 1.0_dp, cos_block(sgfa, sgfb), SIZE(cos_block, 1))
1210 
1211  ! Calculate sine
1212  CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, sinab(1, 1), ldab, &
1213  sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
1214 
1215  CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, &
1216  work(1, 1), ldwork, 1.0_dp, sin_block(sgfa, sgfb), SIZE(sin_block, 1))
1217 
1218  END SUBROUTINE contract_all
1219 
1220 END MODULE qs_efield_berry
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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:195
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
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:339
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:748
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:355
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition: cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition: cp_cfm_types.F:159
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:179
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_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:744
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:208
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:887
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:118
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:167
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:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
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:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
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:23
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, 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, 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, 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, rhs)
Set the QUICKSTEP environment.
Definition: qs_environment_types.F:1035
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, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, 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_r3d_rs_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, 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, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
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, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
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:397
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:1336
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:549
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:139
virial_types
Definition: virial_types.F:13