Octopus
casida.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!! Copyright (C) 2012-2013 D. Strubbe
3!! Copyright (C) 2017-2018 J. Flick, S. Ohlmann
4!!
5!! This program is free software; you can redistribute it and/or modify
6!! it under the terms of the GNU General Public License as published by
7!! the Free Software Foundation; either version 2, or (at your option)
8!! any later version.
9!!
10!! This program is distributed in the hope that it will be useful,
11!! but WITHOUT ANY WARRANTY; without even the implied warranty of
12!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!! GNU General Public License for more details.
14!!
15!! You should have received a copy of the GNU General Public License
16!! along with this program; if not, write to the Free Software
17!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18!! 02110-1301, USA.
19!!
20
21#include "global.h"
22
24
25module casida_oct_m
26 use batch_oct_m
29 use comm_oct_m
30 use debug_oct_m
35#ifdef HAVE_ELPA
36 use elpa
37#endif
40 use forces_oct_m
42 use global_oct_m
43 use grid_oct_m
45 use io_oct_m
47 use iso_c_binding
48 use, intrinsic :: iso_fortran_env
53 use lda_u_oct_m
54 use loct_oct_m
56 use mesh_oct_m
59 use mpi_oct_m
63 use parser_oct_m
64 use pblas_oct_m
65 use pcm_oct_m
74 use sort_oct_m
75 use space_oct_m
81 use unit_oct_m
83 use utils_oct_m
85 use v_ks_oct_m
86 use xc_oct_m
87 use xc_sic_oct_m
88 use xc_f03_lib_m
89
90 implicit none
91
92 private
93 public :: &
94 casida_run, &
96
97 integer, parameter :: &
98 CASIDA_EPS_DIFF = 1, &
101 casida_variational = 8, &
102 casida_casida = 16
103
104 integer, parameter :: &
105 SOLVER_ELPA = 1, &
107
109 type casida_t
110 private
111 integer :: type
112 ! !< CASIDA_VARIATIONAL | CASIDA_CASIDA
113
114 logical :: states_are_real
115 integer, allocatable :: n_occ(:)
116 integer, allocatable :: n_unocc(:)
117 integer :: nst
118 integer :: nik
119 integer :: space_dim
120 integer :: el_per_state
121 character(len=80) :: trandens
122 character(len=80) :: print_exst
123 real(real64) :: weight_thresh
124 logical :: triplet
125 logical :: calc_forces
126 logical :: calc_forces_kernel
127 logical :: calc_forces_scf
128 logical :: herm_conj
129 type(restart_t) :: restart_load
130 type(restart_t) :: restart_dump
131
132 logical, allocatable :: is_included(:,:,:)
133 integer :: n_pairs
134 type(states_pair_t), allocatable :: pair(:)
135 integer, allocatable :: index(:,:,:)
136 integer, allocatable :: ind(:)
137
138 real(real64), allocatable :: dmat(:,:)
139 real(real64), allocatable :: dmat_save(:,:)
140 complex(real64), allocatable :: zmat(:,:)
141 complex(real64), allocatable :: zmat_save(:,:)
142 real(real64), allocatable :: w(:)
143 real(real64), allocatable :: dtm(:, :)
144 complex(real64), allocatable :: ztm(:, :)
145 real(real64), allocatable :: f(:)
146 real(real64), allocatable :: s(:)
147
148 real(real64), allocatable :: rho(:,:)
149 real(real64), allocatable :: fxc(:,:,:)
150 real(real64) :: kernel_lrc_alpha
151 logical :: gga
152 real(real64), allocatable :: fxc_grad(:,:,:,:,:)
153 real(real64), allocatable :: fxc_grad_spin(:,:,:,:)
154
155 real(real64), allocatable :: dmat2(:,:)
156 complex(real64), allocatable :: zmat2(:,:)
157 real(real64), allocatable :: dlr_hmat2(:,:)
158 complex(real64), allocatable :: zlr_hmat2(:,:)
159 real(real64), allocatable :: forces(:,:,:)
160 real(real64), allocatable :: dw2(:)
161 real(real64), allocatable :: zw2(:)
162
163 ! variables for momentum-transfer-dependent calculation
164 logical :: qcalc
165 real(real64), allocatable :: qvector(:)
166 real(real64), allocatable :: qf(:)
167 real(real64), allocatable :: qf_avg(:)
168 integer :: avg_order
169
170 logical :: parallel_in_eh_pairs
171 logical :: parallel_in_domains
172 logical :: distributed_matrix
173 logical :: write_matrix
174 integer :: parallel_solver
175 type(mpi_grp_t) :: mpi_grp
176 logical :: fromScratch
177 logical :: has_photons
178 integer :: pt_nmodes
179 type(photon_mode_t), pointer :: photon_modes => null()
180
181 integer :: n, nb_rows, nb_cols, block_size
182 type(blacs_proc_grid_t) :: proc_grid
183 integer :: desc(BLACS_DLEN)
184 type(MPI_Datatype) :: darray
185 end type casida_t
186
188 private
189 integer :: qi
190 integer :: qa
191 integer :: qk
192 real(real64), allocatable :: dpot(:)
193 complex(real64), allocatable :: zpot(:)
194 end type casida_save_pot_t
195
196contains
198 subroutine casida_run_init()
199
200 push_sub(casida_run_init)
201
202 ! Pure 'other' parallelization is a bad idea. Trying to solve the Poisson equation separately on each node
203 ! consumes excessive memory and time (easily more than is available). In principle, the line below would setup
204 ! joint domain/other parallelization, but 'other' parallelization takes precedence, especially since
205 ! multicomm_init does not know the actual problem size and uses a fictitious value of 10000, making it
206 ! impossible to choose joint parallelization wisely, and generally resulting in a choice of only one domain
207 ! group. FIXME! --DAS
208 ! With the recent improvements, the 'other' parallelization over electron-hole
209 ! pairs works quite well now. For smaller matrices, a combination
210 ! seems to give the fastest run times. For larger matrices (more than a few
211 ! thousand entries per dimension), CasidaDistributedMatrix is needed for
212 ! the Casida matrix to fit into memory; this takes the cores from the
213 ! 'other' strategy. For very large matrices (more than 100000), it is
214 ! advisable to use only the 'other' strategy because the diagonalization
215 ! uses most of the computation time.
216 ! Thus you may want to enable this or a combination of other and domain to get
217 ! better performance - STO
219 call calc_mode_par%set_parallelization(p_strategy_other, default = .false.) ! enabled, but not default
221 call calc_mode_par%unset_parallelization(p_strategy_kpoints) ! disabled. FIXME: could be implemented.
224 end subroutine casida_run_init
226 ! ---------------------------------------------------------
227 subroutine casida_run(system, from_scratch)
228 class(*), intent(inout) :: system
229 logical, intent(in) :: from_scratch
230
231 push_sub(casida_run)
233 select type (system)
235 message(1) = "CalculationMode = casida not implemented for multi-system calculations"
236 call messages_fatal(1, namespace=system%namespace)
237 type is (electrons_t)
238 call casida_run_legacy(system, from_scratch)
239 end select
240
241 pop_sub(casida_run)
242 end subroutine casida_run
244 ! ---------------------------------------------------------
245 subroutine casida_run_legacy(sys, fromScratch)
246 type(electrons_t), intent(inout) :: sys
247 logical, intent(in) :: fromscratch
249 type(casida_t) :: cas
250 type(block_t) :: blk
251 integer :: idir, theorylevel, iatom, ierr, default_int
252 character(len=100) :: restart_filename
253 logical :: is_frac_occ
254 type(restart_t) :: gs_restart
255
256 push_sub(casida_run_legacy)
257 call profiling_in('CASIDA')
259 if (sys%hm%pcm%run_pcm) then
260 call messages_not_implemented("PCM for CalculationMode /= gs or td", namespace=sys%namespace)
261 end if
262
263 if (sys%space%is_periodic()) then
264 message(1) = "Casida oscillator strengths will be incorrect in periodic systems."
265 call messages_warning(1, namespace=sys%namespace)
266 end if
268 if (kpoints_number(sys%kpoints) > 1) then
269 ! Hartree matrix elements may not be correct, not tested anyway. --DAS
270 call messages_not_implemented("Casida with k-points", namespace=sys%namespace)
271 end if
272 if (family_is_mgga_with_exc(sys%hm%xc)) then
273 call messages_not_implemented("Casida with MGGA and non-local terms", namespace=sys%namespace)
274 end if
275 if (sys%hm%lda_u_level /= dft_u_none) then
276 call messages_not_implemented("Casida with DFT+U", namespace=sys%namespace)
277 end if
278 if (sys%hm%theory_level == hartree_fock) then
279 call messages_not_implemented("Casida for Hartree-Fock", namespace=sys%namespace)
280 end if
281 if (sys%hm%theory_level == generalized_kohn_sham_dft) then
282 call messages_not_implemented("Casida for generalized Kohn-Sham", namespace=sys%namespace)
283 end if
285 message(1) = 'Info: Starting Casida linear-response calculation.'
286 call messages_info(1, namespace=sys%namespace)
287
288 call restart_init(gs_restart, sys%namespace, restart_gs, restart_type_load, sys%mc, ierr, mesh=sys%gr, exact=.true.)
289 if (ierr == 0) then
290 call states_elec_look_and_load(gs_restart, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints)
291 call restart_end(gs_restart)
292 else
293 message(1) = "Previous gs calculation is required."
294 call messages_fatal(1, namespace=sys%namespace)
295 end if
296
297 cas%el_per_state = sys%st%smear%el_per_state
298 cas%nst = sys%st%nst
299 cas%nik = sys%st%nik
300 cas%space_dim = sys%space%dim
301 safe_allocate(cas%n_occ(1:sys%st%nik))
302 safe_allocate(cas%n_unocc(1:sys%st%nik))
303
304 call states_elec_count_pairs(sys%st, sys%namespace, cas%n_pairs, cas%n_occ, cas%n_unocc, cas%is_included, is_frac_occ)
305 if (is_frac_occ) then
306 call messages_not_implemented("Casida with partial occupations", namespace=sys%namespace)
307 ! Formulas are in Casida 1995 reference. The occupations are not used at all here currently.
308 end if
309
310 select case (sys%st%d%ispin)
311 case (unpolarized, spinors)
312 write(message(1),'(a,i4,a)') "Info: Found", cas%n_occ(1), " occupied states."
313 write(message(2),'(a,i4,a)') "Info: Found", cas%n_unocc(1), " unoccupied states."
314 call messages_info(2, namespace=sys%namespace)
315 case (spin_polarized)
316 write(message(1),'(a,i4,a)') "Info: Found", cas%n_occ(1), " occupied states with spin up."
317 write(message(2),'(a,i4,a)') "Info: Found", cas%n_unocc(1), " unoccupied states with spin up."
318 write(message(3),'(a,i4,a)') "Info: Found", cas%n_occ(2), " occupied states with spin down."
319 write(message(4),'(a,i4,a)') "Info: Found", cas%n_unocc(2), " unoccupied states with spin down."
320 call messages_info(4, namespace=sys%namespace)
321 end select
322
323
324 ! setup Hamiltonian, without recalculating eigenvalues (use the ones from the restart information)
325 message(1) = 'Info: Setting up Hamiltonian.'
326 call messages_info(1, namespace=sys%namespace)
327 call v_ks_h_setup(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, &
328 sys%hm, calc_eigenval=.false.)
329
330 !%Variable CasidaTheoryLevel
331 !%Type flag
332 !%Section Linear Response::Casida
333 !%Default <tt>eps_diff + petersilka + lrtddft_casida</tt>
334 !%Description
335 !% Choose which electron-hole matrix-based theory levels to use in calculating excitation energies.
336 !% More than one may be used to take advantage of the significant commonality between the calculations.
337 !% <tt>variational</tt> and <tt>lrttdft_casida</tt> are not usable with complex wavefunctions.
338 !% Note the restart data saved by each theory level is compatible with all the others.
339 !%Option eps_diff 1
340 !% Difference of eigenvalues, <i>i.e.</i> independent-particle approximation.
341 !%Option petersilka 2
342 !% The Petersilka approximation uses only elements of the Tamm-Dancoff matrix between degenerate
343 !% transitions (if no degeneracy, this is just the diagonal elements). Also called the "single-pole" approximation.
344 !% This is acceptable if there is little mixing between single-particle transitions.
345 !% Ref: M Petersilka, UJ Gossmann, and EKU Gross, <i>Phys. Rev. Lett.</i> <b>76</b>, 1212 (1996);
346 !% T Grabo, M Petersilka,and EKU Gross, <i>Theochem</i> <b>501-502</b> 353 (2000).
347 !%Option tamm_dancoff 4
348 !% The Tamm-Dancoff approximation uses only occupied-unoccupied transitions and not
349 !% unoccupied-occupied transitions.
350 !% Ref: S Hirata and M Head-Gordon, <i>Chem. Phys. Lett.</i> <b>314</b>, 291 (1999).
351 !%Option variational 8
352 !% Second-order constrained variational theory CV(2)-DFT. Only applies to real wavefunctions.
353 !% Ref: T Ziegler, M Seth, M Krykunov, J Autschbach, and F Wang,
354 !% <i>J. Chem. Phys.</i> <b>130</b>, 154102 (2009).
355 !%Option lrtddft_casida 16
356 !% The full Casida method. Only applies to real wavefunctions.
357 !% Ref: C Jamorski, ME Casida, and DR Salahub, <i>J. Chem. Phys.</i> <b>104</b>, 5134 (1996)
358 !% and ME Casida, "Time-dependent density functional response theory for molecules,"
359 !% in <i>Recent Advances in Density Functional Methods</i>, edited by DE Chong, vol. 1
360 !% of <i>Recent Advances in Computational Chemistry</i>, pp. 155-192 (World Scientific,
361 !% Singapore, 1995).
362 !%End
363
364 call parse_variable(sys%namespace, 'CasidaTheoryLevel', casida_eps_diff + casida_petersilka + casida_casida, theorylevel)
365
366 if (states_are_complex(sys%st)) then
367 if ((bitand(theorylevel, casida_variational) /= 0 &
368 .or. bitand(theorylevel, casida_casida) /= 0)) then
369 message(1) = "Variational and full Casida theory levels do not apply to complex wavefunctions."
370 call messages_fatal(1, only_root_writes = .true., namespace=sys%namespace)
371 ! see section II.D of CV(2) paper regarding this assumption. Would be Eq. 30 with complex wfns.
372 end if
373 end if
374
375 ! This variable is documented in xc_oep_init.
376 call parse_variable(sys%namespace, 'EnablePhotons', .false., cas%has_photons)
377 cas%pt_nmodes = 0
378 if (cas%has_photons) then
379 call messages_experimental('EnablePhotons = yes', namespace=sys%namespace)
380 cas%photon_modes => sys%photons%modes
381 call photon_mode_set_n_electrons(cas%photon_modes, sys%st%qtot)
382 write(message(1), '(a,i7,a)') 'INFO: Solving Casida equation with ', &
383 cas%photon_modes%nmodes, ' photon modes.'
384 write(message(2), '(a)') 'as described in ACS Photonics 2019, 6, 11, 2757-2778.'
385 call messages_info(2, namespace=sys%namespace)
386 cas%pt_nmodes = cas%photon_modes%nmodes
387 end if
388
389 !%Variable CasidaTransitionDensities
390 !%Type string
391 !%Section Linear Response::Casida
392 !%Default write none
393 !%Description
394 !% Specifies which transition densities are to be calculated and written down. The
395 !% transition density for the many-body state <i>n</i> will be written to a file called
396 !% <tt>rho_0n</tt> prefixed by the theory level. Format is set by <tt>OutputFormat</tt>.
397 !%
398 !% This variable is a string in list form, <i>i.e.</i> expressions such as "1,2-5,8-15" are
399 !% valid.
400 !%End
401 call parse_variable(sys%namespace, 'CasidaTransitionDensities', "0", cas%trandens)
402
403 if (cas%trandens /= "0") then
404 call io_function_read_what_how_when(sys%namespace, sys%space, sys%outp%what,&
405 sys%outp%how, sys%outp%output_interval)
406 end if
407
408 !%Variable CasidaMomentumTransfer
409 !%Type block
410 !%Section Linear Response::Casida
411 !%Default 0
412 !%Description
413 !% Momentum-transfer vector for the calculation of the dynamic structure
414 !% factor. When this variable is set, the transition rates are determined
415 !% using an exponential operator instead of the normal dipole one.
416 !%End
417
418 safe_allocate(cas%qvector(1:cas%space_dim))
419 if (parse_block(sys%namespace, 'CasidaMomentumTransfer', blk) == 0) then
420 do idir = 1, cas%space_dim
421 call parse_block_float(blk, 0, idir - 1, cas%qvector(idir))
422 end do
423 call parse_block_end(blk)
424 call messages_experimental("IXS/EELS transition rate calculation", namespace=sys%namespace)
425 message(1) = "Info: Calculating IXS/EELS transition rates."
426 call messages_info(1, namespace=sys%namespace)
427 cas%qcalc = .true.
428
429 !%Variable CasidaQuadratureOrder
430 !%Type integer
431 !%Section Linear Response::Casida
432 !%Default 5
433 !%Description
434 !% Only applies if <tt>CasidaMomentumTransfer</tt> is nonzero.
435 !% Directionally averaged dynamic structure factor is calculated by
436 !% averaging over the results from a set of <math>\vec{q}</math>-vectors. The vectors
437 !% are generated using Gauss-Legendre quadrature scheme [see <i>e.g.</i>
438 !% K. Atkinson, <i>J. Austral. Math. Soc.</i> <b>23</b>, 332 (1982)], and this
439 !% variable determines the order of the scheme.
440 !%End
441 call parse_variable(sys%namespace, 'CasidaQuadratureOrder', 5, cas%avg_order)
442 else
443 cas%qvector(:) = m_zero
444 cas%qcalc = .false.
445 end if
446
447 !%Variable CasidaCalcTriplet
448 !%Type logical
449 !%Section Linear Response::Casida
450 !%Default false
451 !%Description
452 !% For a non-spin-polarized ground state, singlet or triplet excitations can be calculated
453 !% using different matrix elements. Default is to calculate singlets. This variable has no
454 !% effect for a spin-polarized calculation.
455 !%End
456 if (sys%st%d%ispin == unpolarized) then
457 call parse_variable(sys%namespace, 'CasidaCalcTriplet', .false., cas%triplet)
458 else
459 cas%triplet = .false.
460 end if
461
462 if (cas%triplet) then
463 message(1) = "Info: Using triplet kernel. Oscillator strengths will be for spin magnetic-dipole field."
464 call messages_info(1, namespace=sys%namespace)
465 end if
466
467 !%Variable CasidaHermitianConjugate
468 !%Type logical
469 !%Section Linear Response::Casida
470 !%Default false
471 !%Description
472 !% The Casida matrix is Hermitian, so it should not matter whether we calculate the upper or
473 !% lower diagonal. Numerical issues may cause small differences however. Use this variable to
474 !% calculate the Hermitian conjugate of the usual matrix, for testing.
475 !%End
476 call parse_variable(sys%namespace, 'CasidaHermitianConjugate', .false., cas%herm_conj)
477
478 !%Variable CasidaDistributedMatrix
479 !%Type logical
480 !%Section Linear Response::Casida
481 !%Default false
482 !%Description
483 !% Large matrices with more than a few thousand rows and columns usually do
484 !% not fit into the memory of one processor anymore. With this option, the
485 !% Casida matrix is distributed in block-cyclic fashion over all cores in the
486 !% ParOther group. The diagonalization is done in parallel using ScaLAPACK
487 !% or ELPA, if available. For very large matrices (>100000), only the
488 !% ParOther strategy should be used because the diagonalization dominates
489 !% the run time of the computation.
490 !%End
491 call parse_variable(sys%namespace, 'CasidaDistributedMatrix', .false., cas%distributed_matrix)
492#ifndef HAVE_SCALAPACK
493 if (cas%distributed_matrix) then
494 message(1) = "ScaLAPACK layout requested, but code not compiled with ScaLAPACK"
495 call messages_fatal(1, namespace=sys%namespace)
496 end if
497#endif
498 call messages_obsolete_variable(sys%namespace, 'CasidaUseScalapackLayout', 'CasidaDistributedMatrix')
499
500 !%Variable CasidaWriteDistributedMatrix
501 !%Type logical
502 !%Section Linear Response::Casida
503 !%Default false
504 !%Description
505 !% Set to true to write out the full distributed Casida matrix to a file
506 !% using MPI-IO.
507 !%End
508 call parse_variable(sys%namespace, 'CasidaWriteDistributedMatrix', .false., cas%write_matrix)
509 if (.not. cas%distributed_matrix .and. cas%write_matrix) then
510 message(1) = "CasidaWriteDistributedMatrix con only be used with CasidaDistributedMatrix"
511 call messages_fatal(1, namespace=sys%namespace)
512 end if
513
514 !%Variable CasidaParallelEigensolver
515 !%Type integer
516 !%Section Linear Response::Casida
517 !%Description
518 !% Choose library to use for solving the parallel eigenproblem
519 !% of the Casida problem. This options is only relevant if a
520 !% distributed matrix is used (CasidaDistributedMatrix=true).
521 !% By default, elpa is chosen if available.
522 !%Option casida_elpa 1
523 !% Use ELPA library as parallel eigensolver
524 !%Option casida_scalapack 2
525 !% Use Scalapack as parallel eigensolver
526 !%End
527#ifdef HAVE_ELPA
528 default_int = solver_elpa
529#else
530 default_int = solver_scalapack
531#endif
532 call parse_variable(sys%namespace, 'CasidaParallelEigensolver', default_int, cas%parallel_solver)
533 if (.not. varinfo_valid_option('CasidaParallelEigensolver', cas%parallel_solver)) then
534 call messages_input_error(sys%namespace, 'CasidaParallelEigensolver')
535 end if
536#ifndef HAVE_ELPA
537 if (cas%distributed_matrix .and. cas%parallel_solver == solver_elpa) then
538 message(1) = "ELPA solver requested, but code not compiled with ELPA"
539 call messages_fatal(1, namespace=sys%namespace)
540 end if
541#endif
542
543 !%Variable CasidaPrintExcitations
544 !%Type string
545 !%Section Linear Response::Casida
546 !%Default write all
547 !%Description
548 !% Specifies which excitations are written at the end of the calculation.
549 !%
550 !% This variable is a string in list form, <i>i.e.</i> expressions such as "1,2-5,8-15" are
551 !% valid.
552 !%End
553 call parse_variable(sys%namespace, 'CasidaPrintExcitations', "all", cas%print_exst)
554 if (cas%distributed_matrix) then
555 ! do not print excited states -> too many files generated!
556 cas%print_exst = "none"
557 message(1) = "Using ScaLAPACK layout, thus disabling output of excited states."
558 message(2) = "This options creates too many files for large Casida matrices."
559 call messages_info(2, namespace=sys%namespace)
560 end if
561
562 !%Variable CasidaWeightThreshold
563 !%Type float
564 !%Section Linear Response::Casida
565 !%Default -1.
566 !%Description
567 !% Specifies the threshold value for which the individual excitations are printed.
568 !% i.e. juste-h pairs with weight larger than this threshold will be printed.
569 !%
570 !% If a negative value (default) is set, all coefficients will be printed.
571 !% For many case, a 0.01 value is a valid option.
572 !%End
573 call parse_variable(sys%namespace, 'CasidaWeightThreshold', -m_one, cas%weight_thresh)
574 if (cas%weight_thresh > m_one) then
575 message(1) = 'Casida coefficients have values between 0 and 1'
576 message(2) = 'Threshold values reset to default value'
577 call messages_warning(2, namespace=sys%namespace)
578 cas%weight_thresh = -m_one
579 end if
580
581 !%Variable CasidaCalcForces
582 !%Type logical
583 !%Section Linear Response::Casida
584 !%Default false
585 !%Description
586 !% (Experimental) Enable calculation of excited-state forces. Requires previous <tt>vib_modes</tt> calculation.
587 !%End
588 call parse_variable(sys%namespace, 'CasidaCalcForces', .false., cas%calc_forces)
589 if (cas%calc_forces) then
590 call messages_experimental("Excited-state forces calculation", namespace=sys%namespace)
591
592 !%Variable CasidaCalcForcesKernel
593 !%Type logical
594 !%Section Linear Response::Casida
595 !%Default true
596 !%Description
597 !% If false, the derivative of the kernel will not be included in the excited-state force calculation.
598 !%End
599 call parse_variable(sys%namespace, 'CasidaCalcForcesKernel', .true., cas%calc_forces_kernel)
600
601 !%Variable CasidaCalcForcesSCF
602 !%Type logical
603 !%Section Linear Response::Casida
604 !%Default false
605 !%Description
606 !% If true, the ground-state forces will be included in the excited-state forces, so they are total forces.
607 !% If false, the excited-state forces that are produced are only the gradients of the excitation energy.
608 !%End
609 call parse_variable(sys%namespace, 'CasidaCalcForcesSCF', .false., cas%calc_forces_scf)
610
611 if (cas%distributed_matrix) then
612 message(1) = "Info: Forces calculation not compatible with ScaLAPACK layout."
613 message(2) = "Using normal layout."
614 call messages_info(2, namespace=sys%namespace)
615 cas%distributed_matrix = .false.
616 end if
617 end if
618
619 ! Initialize structure
620 call casida_type_init(cas, sys)
621
622 cas%fromScratch = fromscratch
623
624 if (cas%fromScratch) then ! remove old restart files
625 if (cas%triplet) then
626 call restart_rm(cas%restart_dump, 'kernel_triplet')
627 else
628 call restart_rm(cas%restart_dump, 'kernel')
629 end if
630
631 if (cas%calc_forces) then
632 do iatom = 1, sys%ions%natoms
633 do idir = 1, cas%space_dim
634 write(restart_filename,'(a,i6.6,a,i1)') 'lr_kernel_', iatom, '_', idir
635 if (cas%triplet) restart_filename = trim(restart_filename)//'_triplet'
636 call restart_rm(cas%restart_dump, restart_filename)
637
638 write(restart_filename,'(a,i6.6,a,i1)') 'lr_hmat1_', iatom, '_', idir
639 call restart_rm(cas%restart_dump, restart_filename)
640 end do
641 end do
642 end if
643 end if
644
645 ! First, print the differences between KS eigenvalues (first approximation to the excitation energies).
646 if (bitand(theorylevel, casida_eps_diff) /= 0) then
647 message(1) = "Info: Approximating resonance energies through KS eigenvalue differences"
648 call messages_info(1, namespace=sys%namespace)
649 cas%type = casida_eps_diff
650 call casida_work(sys, cas)
651 end if
652
653 if (sys%st%d%ispin /= spinors) then
654
655 if (bitand(theorylevel, casida_tamm_dancoff) /= 0) then
656 call messages_experimental("Tamm-Dancoff calculation", namespace=sys%namespace)
657 message(1) = "Info: Calculating matrix elements in the Tamm-Dancoff approximation"
658 call messages_info(1, namespace=sys%namespace)
659 cas%type = casida_tamm_dancoff
660 call casida_work(sys, cas)
661 end if
662
663 if (bitand(theorylevel, casida_variational) /= 0) then
664 call messages_experimental("CV(2)-DFT calculation", namespace=sys%namespace)
665 message(1) = "Info: Calculating matrix elements with the CV(2)-DFT theory"
666 call messages_info(1, namespace=sys%namespace)
667 cas%type = casida_variational
668 call casida_work(sys, cas)
669 end if
670
671 if (bitand(theorylevel, casida_casida) /= 0) then
672 message(1) = "Info: Calculating matrix elements with the full Casida method"
673 call messages_info(1, namespace=sys%namespace)
674 cas%type = casida_casida
675 call casida_work(sys, cas)
676 end if
677
678 ! Doing this first, if doing the others later, takes longer, because we would use
679 ! each Poisson solution for only one matrix element instead of a whole column.
680 if (bitand(theorylevel, casida_petersilka) /= 0) then
681 message(1) = "Info: Calculating resonance energies via the Petersilka approximation"
682 call messages_info(1, namespace=sys%namespace)
683 cas%type = casida_petersilka
684 call casida_work(sys, cas)
685 end if
686
687 end if
688
689 call casida_type_end(cas)
690
691 call profiling_out('CASIDA')
692 pop_sub(casida_run_legacy)
693 end subroutine casida_run_legacy
694
695 ! ---------------------------------------------------------
697 subroutine casida_type_init(cas, sys)
698 type(casida_t), intent(inout) :: cas
699 type(electrons_t), intent(in) :: sys
700
701 integer :: ist, ast, jpair, ik, ierr
702#ifdef HAVE_SCALAPACK
703 integer :: np, np_rows, np_cols, ii, info
704#endif
705
706 push_sub(casida_type_init)
707
708 cas%kernel_lrc_alpha = sys%ks%xc%kernel_lrc_alpha
709 cas%states_are_real = states_are_real(sys%st)
710 if (cas%distributed_matrix .and. .not. cas%states_are_real) then
711 call messages_not_implemented("Complex wavefunctions with ScaLAPACK layout", namespace=sys%namespace)
712 end if
713
714 write(message(1), '(a,i9)') "Number of occupied-unoccupied pairs: ", cas%n_pairs
715 call messages_info(1, namespace=sys%namespace)
716
717 if (cas%n_pairs < 1) then
718 message(1) = "No Casida pairs -- maybe there are no unoccupied states?"
719 call messages_fatal(1, only_root_writes = .true., namespace=sys%namespace)
720 end if
721
722 if (mpi_grp_is_root(mpi_world)) write(*, "(1x)")
723
724 ! now let us take care of initializing the parallel stuff
725 cas%parallel_in_eh_pairs = multicomm_strategy_is_parallel(sys%mc, p_strategy_other)
726 if (cas%parallel_in_eh_pairs) then
727 call mpi_grp_init(cas%mpi_grp, sys%mc%group_comm(p_strategy_other))
728 else
729 call mpi_grp_init(cas%mpi_grp, mpi_comm_undefined)
730 end if
731 cas%parallel_in_domains = multicomm_strategy_is_parallel(sys%mc, p_strategy_domains)
732
733 if (cas%distributed_matrix .and. .not. cas%parallel_in_eh_pairs) then
734 message(1) = "ScaLAPACK layout requested, but 'Other' parallelization strategy not available."
735 message(2) = "Please set ParOther to use the ScaLAPACK layout."
736 message(3) = "Continuing without ScaLAPACK layout."
737 call messages_info(3, namespace=sys%namespace)
738 cas%distributed_matrix = .false.
739 end if
740
741 ! dimension of matrix
742 cas%n = cas%n_pairs + cas%pt_nmodes
743
744 ! initialize block-cyclic matrix
745 if (cas%distributed_matrix) then
746#ifdef HAVE_SCALAPACK
747 ! processor layout: always use more processors for rows, this leads to
748 ! better load balancing when computing the matrix elements
749 np = cas%mpi_grp%size
750 np_cols = 1
751 if (np > 3) then
752 do ii = floor(sqrt(real(np))), 2, -1
753 if (mod(np, ii) == 0) then
754 np_cols = ii
755 exit
756 end if
757 end do
758 end if
759 np_rows = np / np_cols
760
761 ! recommended block size: 64, take smaller value for smaller matrices for
762 ! better load balancing
763 cas%block_size = min(64, cas%n / np_rows)
764 ! limit to a minimum block size of 5 for diagonalization efficiency
765 cas%block_size = max(5, cas%block_size)
766 write(message(1), '(A,I5,A,I5,A,I5,A)') 'Parallel layout: using block size of ',&
767 cas%block_size, ' and a processor grid with ', np_rows, 'x', np_cols, &
768 ' processors (rows x cols)'
769 call messages_info(1, namespace=sys%namespace)
770
771 call blacs_proc_grid_init(cas%proc_grid, cas%mpi_grp, procdim = (/np_rows, np_cols/))
772
773 ! get size of local matrices
774 cas%nb_rows = numroc(cas%n, cas%block_size, cas%proc_grid%myrow, 0, cas%proc_grid%nprow)
775 cas%nb_cols = numroc(cas%n, cas%block_size, cas%proc_grid%mycol, 0, cas%proc_grid%npcol)
776
777 ! get ScaLAPACK descriptor
778 call descinit(cas%desc(1), cas%n, cas%n, cas%block_size, cas%block_size, 0, 0, &
779 cas%proc_grid%context, cas%nb_rows, info)
780#endif
781 else
782 ! set to full size
783 cas%nb_rows = cas%n
784 cas%nb_cols = cas%n
785 end if
786
787
788 ! allocate stuff
789 safe_allocate(cas%pair(1:cas%n))
790 if (cas%states_are_real) then
791 safe_allocate( cas%dmat(1:cas%nb_rows, 1:cas%nb_cols))
792 safe_allocate( cas%dtm(1:cas%n, 1:cas%space_dim))
793 else
794 ! caution: ScaLAPACK layout not yet tested for complex wavefunctions!
795 safe_allocate( cas%zmat(1:cas%nb_rows, 1:cas%nb_cols))
796 safe_allocate( cas%ztm(1:cas%n, 1:cas%space_dim))
797 end if
798 safe_allocate( cas%f(1:cas%n))
799 safe_allocate( cas%s(1:cas%n_pairs))
800 safe_allocate( cas%w(1:cas%n))
801 safe_allocate( cas%index(1:maxval(cas%n_occ), cas%nst - maxval(cas%n_unocc) + 1:cas%nst, 1:cas%nik))
802 safe_allocate( cas%ind(1:cas%n))
803
804 if (cas%calc_forces) then
805 if (cas%states_are_real) then
806 safe_allocate(cas%dmat_save(1:cas%n_pairs, 1:cas%n_pairs))
807 else
808 safe_allocate(cas%zmat_save(1:cas%n_pairs, 1:cas%n_pairs))
809 end if
810 safe_allocate(cas%forces(1:cas%space_dim, 1:sys%ions%natoms, 1:cas%n_pairs))
811 end if
812
813 if (cas%qcalc) then
814 safe_allocate( cas%qf (1:cas%n_pairs))
815 safe_allocate( cas%qf_avg(1:cas%n_pairs))
816 end if
817
818 cas%index(:,:,:) = 0
819
820 ! create pairs
821 jpair = 1
822 do ik = 1, cas%nik
823 do ast = cas%n_occ(ik) + 1, cas%nst
824 do ist = 1, cas%n_occ(ik)
825 if (cas%is_included(ist, ast, ik)) then
826 cas%index(ist, ast, ik) = jpair
827 cas%pair(jpair)%i = ist
828 cas%pair(jpair)%a = ast
829 cas%pair(jpair)%kk = ik
830 jpair = jpair + 1
831 end if
832 end do
833 end do
834 end do
835
836 if (cas%has_photons) then
837 ! create pairs for photon modes (negative number refers to photonic excitation)
838 do ik = 1, cas%pt_nmodes
839 cas%pair(cas%n_pairs + ik)%i = 1
840 cas%pair(cas%n_pairs + ik)%a = -ik
841 cas%pair(cas%n_pairs + ik)%kk = -ik
842 end do
843 end if
844
845 safe_deallocate_a(cas%is_included)
846
847 call restart_init(cas%restart_dump, sys%namespace, restart_casida, restart_type_dump, sys%mc, ierr)
848 call restart_init(cas%restart_load, sys%namespace, restart_casida, restart_type_load, sys%mc, ierr)
849
850 pop_sub(casida_type_init)
851 end subroutine casida_type_init
852
853
854 ! ---------------------------------------------------------
855 subroutine casida_type_end(cas)
856 type(casida_t), intent(inout) :: cas
857
858 push_sub(casida_type_end)
859
860 assert(allocated(cas%pair))
861 safe_deallocate_a(cas%pair)
862 safe_deallocate_a(cas%index)
863 if (cas%states_are_real) then
864 safe_deallocate_a(cas%dmat)
865 safe_deallocate_a(cas%dtm)
866 else
867 safe_deallocate_a(cas%zmat)
868 safe_deallocate_a(cas%ztm)
869 end if
870 safe_deallocate_a(cas%s)
871 safe_deallocate_a(cas%f)
872 safe_deallocate_a(cas%w)
873 safe_deallocate_a(cas%ind)
874
875 if (cas%qcalc) then
876 safe_deallocate_a(cas%qf)
877 safe_deallocate_a(cas%qf_avg)
878 end if
879
880 safe_deallocate_a(cas%n_occ)
881 safe_deallocate_a(cas%n_unocc)
882
883 if (cas%calc_forces) then
884 if (cas%states_are_real) then
885 safe_deallocate_a(cas%dmat_save)
886 else
887 safe_deallocate_a(cas%zmat_save)
888 end if
889 safe_deallocate_a(cas%forces)
890 end if
891
892 call restart_end(cas%restart_dump)
893 call restart_end(cas%restart_load)
894
895 if (cas%distributed_matrix) then
896#ifdef HAVE_SCALAPACK
897 call blacs_proc_grid_end(cas%proc_grid)
898#endif
899 end if
900
901 safe_deallocate_a(cas%qvector)
902
903 pop_sub(casida_type_end)
904 end subroutine casida_type_end
905
906
907 ! ---------------------------------------------------------
910 subroutine casida_work(sys, cas)
911 type(electrons_t), target, intent(inout) :: sys
912 type(casida_t), intent(inout) :: cas
913
914 type(states_elec_t), pointer :: st
915 type(grid_t), pointer :: gr
916
917 real(real64), allocatable :: rho_spin(:, :)
918 real(real64), allocatable :: fxc_spin(:,:,:)
919 character(len=100) :: restart_filename
920
921 push_sub(casida_work)
922
923 ! sanity checks
924 assert(cas%type >= casida_eps_diff .and. cas%type <= casida_casida)
925
926 ! some shortcuts
927 st => sys%st
928 gr => sys%gr
929
930 ! initialize stuff
931 if (cas%states_are_real) then
932 cas%dmat = m_zero
933 cas%dtm = m_zero
934 else
935 cas%zmat = m_zero
936 cas%ztm = m_zero
937 end if
938 cas%f = m_zero
939 cas%w = m_zero
940 cas%s = m_zero
941 if (cas%qcalc) then
942 cas%qf = m_zero
943 cas%qf_avg = m_zero
944 end if
945
946 if (cas%type /= casida_eps_diff .or. cas%calc_forces) then
947 ! We calculate here the kernel, since it will be needed later.
948 safe_allocate(cas%rho(1:gr%np_part, 1:st%d%nspin))
949 safe_allocate(cas%fxc(1:gr%np, 1:st%d%nspin, 1:st%d%nspin))
950 cas%gga = in_family(sys%ks%xc%family, [xc_family_gga])
951 if (cas%gga) then
952 safe_allocate(cas%fxc_grad(1:gr%np, 1:gr%der%dim, 1:gr%der%dim, 1:st%d%nspin, 1:st%d%nspin))
953 if (st%d%ispin == spin_polarized) then
954 safe_allocate(cas%fxc_grad_spin(1:gr%np, 1:gr%der%dim, 1:st%d%nspin, 1:st%d%nspin))
955 else
956 safe_allocate(cas%fxc_grad_spin(0, 0, 0, 0))
957 end if
958 else
959 safe_allocate(cas%fxc_grad(0, 0, 0, 0, 0))
960 safe_allocate(cas%fxc_grad_spin(0, 0, 0, 0))
961 end if
962
963 call states_elec_total_density(st, gr, cas%rho)
964 if (cas%triplet) then
965 safe_allocate(rho_spin(1:gr%np_part, 1:2))
966 safe_allocate(fxc_spin(1:gr%np, 1:2, 1:2))
967
968 rho_spin(:, 1) = m_half * cas%rho(:, 1)
969 rho_spin(:, 2) = m_half * cas%rho(:, 1)
970
971 call xc_get_fxc(sys%ks%xc, gr, sys%namespace, rho_spin, spin_polarized, fxc_spin)
972 cas%fxc(:, 1, 1) = m_half * (fxc_spin(:, 1, 1) - fxc_spin(:, 1, 2))
973
974 safe_deallocate_a(rho_spin)
975 safe_deallocate_a(fxc_spin)
976 else
977 if (cas%gga) then
978 call xc_get_fxc(sys%ks%xc, gr, sys%namespace, cas%rho, st%d%ispin, cas%fxc, cas%fxc_grad, &
979 cas%fxc_grad_spin)
980 else
981 call xc_get_fxc(sys%ks%xc, gr, sys%namespace, cas%rho, st%d%ispin, cas%fxc)
982 end if
983 end if
984
985 if (sys%ks%sic%level == sic_adsic) then
986 call fxc_add_adsic(sys%namespace, sys%ks, st, gr, cas)
987 end if
988
989 end if
990
991 restart_filename = 'kernel'
992 if (cas%triplet) restart_filename = trim(restart_filename)//'_triplet'
993
994 select case (cas%type)
995 case (casida_eps_diff)
996 call solve_eps_diff()
998 if (cas%states_are_real) then
999 call dcasida_get_matrix(cas, sys%namespace, sys%hm, st, sys%ks, gr, cas%dmat, cas%fxc, &
1000 cas%fxc_grad, cas%fxc_grad_spin, restart_filename)
1001 call dcasida_solve(cas, sys)
1002 else
1003 call zcasida_get_matrix(cas, sys%namespace, sys%hm, st, sys%ks, gr, cas%zmat, cas%fxc, &
1004 cas%fxc_grad, cas%fxc_grad_spin, restart_filename)
1005 call zcasida_solve(cas, sys)
1006 end if
1007 end select
1008
1009 ! compute oscillator strengths on all processes for the ScaLAPACK layout
1010 if (mpi_grp_is_root(cas%mpi_grp) .or. cas%distributed_matrix) then
1011 if (cas%states_are_real) then
1012 call doscillator_strengths(cas, gr, st)
1013 else
1014 call zoscillator_strengths(cas, gr, st)
1015 end if
1016 end if
1017
1018 if (cas%calc_forces) then
1019 if (cas%states_are_real) then
1020 call dcasida_forces(cas, sys, gr, st)
1021 else
1022 call zcasida_forces(cas, sys, gr, st)
1023 end if
1024 end if
1025
1026 if (cas%states_are_real) then
1027 call dcasida_write(cas, sys)
1028 else
1029 call zcasida_write(cas, sys)
1030 end if
1031
1032 ! clean up
1033 if (cas%type /= casida_eps_diff .or. cas%calc_forces) then
1034 safe_deallocate_a(cas%fxc)
1035 safe_deallocate_a(cas%fxc_grad)
1036 safe_deallocate_a(cas%fxc_grad_spin)
1037 safe_deallocate_a(cas%rho)
1038 end if
1039
1040 pop_sub(casida_work)
1041
1042 contains
1043
1044 ! ---------------------------------------------------------
1045 subroutine solve_eps_diff
1046
1047 integer :: ia
1048 real(real64), allocatable :: w(:)
1049
1050 push_sub(casida_work.solve_eps_diff)
1051
1052 ! initialize progress bar
1053 if (mpi_grp_is_root(mpi_world)) call loct_progress_bar(-1, cas%n_pairs)
1054
1055 do ia = 1, cas%n_pairs
1056 cas%w(ia) = st%eigenval(cas%pair(ia)%a, cas%pair(ia)%kk) - &
1057 st%eigenval(cas%pair(ia)%i, cas%pair(ia)%kk)
1058 if (cas%w(ia) < -m_epsilon) then
1059 message(1) = "There is a negative unocc-occ KS eigenvalue difference for"
1060 write(message(2),'("states ",I5," and ",I5," of k-point ",I5,".")') cas%pair(ia)%i, cas%pair(ia)%a, cas%pair(ia)%kk
1061 message(3) = "This indicates an inconsistency between gs, unocc, and/or casida calculations."
1062 call messages_fatal(3, only_root_writes = .true., namespace=sys%namespace)
1063 end if
1064 if (mpi_grp_is_root(mpi_world)) call loct_progress_bar(ia, cas%n_pairs)
1065 end do
1066
1067 safe_allocate(w(1:cas%n_pairs))
1068 w = cas%w
1069 call sort(w, cas%ind)
1070 safe_deallocate_a(w)
1071
1072 if (mpi_grp_is_root(mpi_world)) write(*, "(1x)")
1073
1075 end subroutine solve_eps_diff
1076
1077 ! ---------------------------------------------------------
1078 subroutine fxc_add_adsic(namespace, ks, st, gr, cas)
1079 type(namespace_t), intent(in) :: namespace
1080 type(v_ks_t), intent(in) :: ks
1081 type(states_elec_t), intent(in) :: st
1082 type(grid_t), intent(in) :: gr
1083 type(casida_t), intent(inout) :: cas
1084
1085 real(real64), allocatable :: rho(:, :)
1086 real(real64), allocatable :: fxc_sic(:,:,:)
1087
1088 push_sub(casida_work.fxc_add_adsic)
1089
1090 !Check spin and triplets
1091 if (st%d%ispin /= unpolarized) then
1092 message(1) = "Casida calculation with ADSIC not implemented for spin-polarized calculations."
1093 call messages_fatal(1, namespace=sys%namespace)
1094 end if
1095 if (cas%triplet) then
1096 message(1) = "Casida calculation with ADSIC not implemented for triplet excitations."
1097 call messages_fatal(1, namespace=sys%namespace)
1098 end if
1099 if (cas%gga) then
1100 message(1) = "Casida calculation with ADSIC not implemented with GGAs."
1101 call messages_fatal(1, namespace=sys%namespace)
1102 end if
1103
1104 safe_allocate(fxc_sic(1:gr%np, 1:st%d%nspin, 1:st%d%nspin))
1105 safe_allocate(rho(1:gr%np, 1:st%d%nspin))
1106 rho = cas%rho/st%qtot
1107
1108 call xc_get_fxc(ks%xc, gr, namespace, rho, 1, fxc_sic)
1109
1110 cas%fxc = cas%fxc - fxc_sic/st%qtot
1111
1112 safe_deallocate_a(rho)
1113 safe_deallocate_a(fxc_sic)
1114
1115 pop_sub(casida_work.fxc_add_adsic)
1116 end subroutine fxc_add_adsic
1117
1118 end subroutine casida_work
1119
1120 ! ---------------------------------------------------------
1121 real(real64) function casida_matrix_factor(cas, sys)
1122 type(casida_t), intent(in) :: cas
1123 type(electrons_t), intent(in) :: sys
1124
1125 push_sub(casida_matrix_factor)
1126
1127 casida_matrix_factor = m_one
1128
1129 if (cas%type == casida_variational) then
1130 casida_matrix_factor = m_two * casida_matrix_factor
1131 end if
1132
1133 if (sys%st%d%ispin == unpolarized) then
1134 casida_matrix_factor = m_two * casida_matrix_factor
1135 end if
1136
1137 pop_sub(casida_matrix_factor)
1138
1139 end function casida_matrix_factor
1140
1141 ! ---------------------------------------------------------
1142 subroutine qcasida_write(cas, namespace)
1143 type(casida_t), intent(in) :: cas
1144 type(namespace_t), intent(in) :: namespace
1145
1146 integer :: iunit, ia
1147
1148 if (.not. mpi_grp_is_root(mpi_world)) return
1149
1150 push_sub(qcasida_write)
1151
1152 call io_mkdir(casida_dir, namespace)
1153 iunit = io_open(casida_dir//'q'//trim(theory_name(cas)), namespace, action='write')
1154 write(iunit, '(a1,a14,1x,a24,1x,a24,1x,a10,3es15.8,a2)') '#','E' , '|<f|exp(iq.r)|i>|^2', &
1155 '<|<f|exp(iq.r)|i>|^2>','; q = (',cas%qvector(1:cas%space_dim),')'
1156 write(iunit, '(a1,a14,1x,a24,1x,a24,1x,10x,a15)') '#', trim(units_abbrev(units_out%energy)), &
1157 trim('-'), &
1158 trim('-'), &
1159 trim('a.u.')
1160
1161 if (cas%avg_order == 0) then
1162 do ia = 1, cas%n_pairs
1163 write(iunit, '(es15.8,es15.8)') units_from_atomic(units_out%energy, cas%w(cas%ind(ia))), cas%qf(cas%ind(ia))
1164 end do
1165 else
1166 do ia = 1, cas%n_pairs
1167 write(iunit, '(3es15.8)') units_from_atomic(units_out%energy, cas%w(cas%ind(ia))), &
1168 cas%qf (cas%ind(ia)), &
1169 cas%qf_avg(cas%ind(ia))
1170 end do
1171 end if
1172
1173 call io_close(iunit)
1174
1175 pop_sub(qcasida_write)
1176
1177 end subroutine qcasida_write
1178
1179 ! ---------------------------------------------------------
1180 character(len=80) pure function theory_name(cas)
1181 type(casida_t), intent(in) :: cas
1182
1183 select case (cas%type)
1184 case (casida_eps_diff)
1185 theory_name = "eps_diff"
1186 case (casida_petersilka)
1187 theory_name = "petersilka"
1188 case (casida_tamm_dancoff)
1189 theory_name = "tamm_dancoff"
1190 case (casida_variational)
1191 theory_name = "variational"
1192 case (casida_casida)
1193 theory_name = "casida"
1194 case default
1195 theory_name = "unknown"
1196 end select
1197
1198 end function theory_name
1199
1200 logical function isnt_degenerate(cas, st, ia, jb)
1201 type(casida_t), intent(in) :: cas
1202 type(states_elec_t), intent(in) :: st
1203 integer, intent(in) :: ia
1204 integer, intent(in) :: jb
1205
1206 push_sub(isnt_degenerate)
1207
1208 isnt_degenerate = (abs((st%eigenval(cas%pair(ia)%a, cas%pair(ia)%kk) - st%eigenval(cas%pair(ia)%i, cas%pair(ia)%kk)) &
1209 - (st%eigenval(cas%pair(jb)%a, cas%pair(jb)%kk) - st%eigenval(cas%pair(jb)%i, cas%pair(jb)%kk))) > 1e-8_real64)
1210
1211 pop_sub(isnt_degenerate)
1212 end function isnt_degenerate
1213
1214 integer function get_global_row(cas, jb_local) result(jb)
1215 implicit none
1216 type(casida_t), intent(inout) :: cas
1217 integer, intent(in) :: jb_local
1218
1219 if (.not. cas%distributed_matrix) then
1220 jb = jb_local
1221 else
1222#ifdef HAVE_SCALAPACK
1223 jb = indxl2g(jb_local, cas%block_size, cas%proc_grid%myrow, 0, cas%proc_grid%nprow)
1224#endif
1225 end if
1226 end function get_global_row
1227
1228 integer function get_global_col(cas, ia_local) result(ia)
1229 implicit none
1230 type(casida_t), intent(inout) :: cas
1231 integer, intent(in) :: ia_local
1232
1233 if (.not. cas%distributed_matrix) then
1234 ia = ia_local
1235 else
1236#ifdef HAVE_SCALAPACK
1237 ia = indxl2g(ia_local, cas%block_size, cas%proc_grid%mycol, 0, cas%proc_grid%npcol)
1238#endif
1239 end if
1240 end function get_global_col
1241
1242 subroutine local_indices(cas, ia, jb, on_this_processor, ia_local, jb_local)
1243 implicit none
1244 type(casida_t), intent(in) :: cas
1245 integer, intent(in) :: ia, jb
1246 logical, intent(out) :: on_this_processor
1247 integer, intent(out) :: ia_local, jb_local
1248#ifdef HAVE_SCALAPACK
1249 integer :: ia_proc, jb_proc
1250#endif
1251
1252 if (.not. cas%distributed_matrix) then
1253 on_this_processor = .true.
1254 ia_local = ia
1255 jb_local = jb
1256 else
1257#ifdef HAVE_SCALAPACK
1258 ia_proc = indxg2p(ia, cas%block_size, cas%proc_grid%mycol, 0, cas%proc_grid%npcol)
1259 jb_proc = indxg2p(jb, cas%block_size, cas%proc_grid%myrow, 0, cas%proc_grid%nprow)
1260 if (cas%proc_grid%mycol == ia_proc .and. cas%proc_grid%myrow == jb_proc) then
1261 on_this_processor = .true.
1262 ia_local = indxg2l(ia, cas%block_size, cas%proc_grid%mycol, 0, cas%proc_grid%npcol)
1263 jb_local = indxg2l(jb, cas%block_size, cas%proc_grid%myrow, 0, cas%proc_grid%nprow)
1264 else
1265 on_this_processor = .false.
1266 ia_local = -1
1267 jb_local = -1
1268 end if
1269#endif
1270 end if
1271 end subroutine local_indices
1272
1273#include "undef.F90"
1274#include "real.F90"
1275#include "casida_inc.F90"
1276#include "undef.F90"
1277#include "complex.F90"
1278#include "casida_inc.F90"
1279
1280end module casida_oct_m
1281
1282!! Local Variables:
1283!! mode: f90
1284!! coding: utf-8
1285!! End:
subroutine fxc_add_adsic(namespace, ks, st, gr, cas)
Definition: casida.F90:1138
subroutine solve_eps_diff
Definition: casida.F90:1105
This is the common interface to a sorting routine. It performs the shell algorithm,...
Definition: sort.F90:149
double floor(double __x) __attribute__((__nothrow__
This module implements batches of mesh functions.
Definition: batch.F90:133
This module provides the BLACS processor grid.
subroutine, public blacs_proc_grid_init(this, mpi_grp, procdim)
Initializes a blacs context from an MPI communicator with topological information.
subroutine, public blacs_proc_grid_end(this)
This module handles the calculation mode.
integer, parameter, public p_strategy_kpoints
parallelization in k-points
integer, parameter, public p_strategy_other
something else like e-h pairs
integer, parameter, public p_strategy_domains
parallelization in domains
type(calc_mode_par_t), public calc_mode_par
Singleton instance of parallel calculation mode.
This module implements the Casida equations for excited states.
Definition: casida.F90:118
integer function get_global_col(cas, ia_local)
Definition: casida.F90:1288
integer, parameter solver_scalapack
Definition: casida.F90:197
subroutine zoscillator_strengths(cas, mesh, st)
Definition: casida.F90:2995
integer, parameter casida_petersilka
Definition: casida.F90:190
subroutine, public casida_run_init()
Definition: casida.F90:292
integer, parameter casida_casida
Definition: casida.F90:190
subroutine casida_type_init(cas, sys)
allocates stuff, and constructs the arrays pair_i and pair_j
Definition: casida.F90:791
integer function get_global_row(cas, jb_local)
Definition: casida.F90:1274
integer, parameter casida_eps_diff
Definition: casida.F90:190
character(len=80) pure function theory_name(cas)
Definition: casida.F90:1240
subroutine dcasida_forces(cas, sys, gr, st)
Definition: casida.F90:2281
subroutine dcasida_get_matrix(cas, namespace, hm, st, ks, gr, matrix, fxc, fxc_grad, fxc_grad_spin, restart_file, is_forces)
Definition: casida.F90:1804
subroutine local_indices(cas, ia, jb, on_this_processor, ia_local, jb_local)
Definition: casida.F90:1302
subroutine zcasida_solve(cas, sys)
Definition: casida.F90:4147
subroutine casida_work(sys, cas)
this subroutine calculates electronic excitation energies using the matrix formulation of M....
Definition: casida.F90:970
subroutine zcasida_get_matrix(cas, namespace, hm, st, ks, gr, matrix, fxc, fxc_grad, fxc_grad_spin, restart_file, is_forces)
Definition: casida.F90:3409
integer, parameter casida_variational
Definition: casida.F90:190
subroutine qcasida_write(cas, namespace)
Definition: casida.F90:1202
subroutine zcasida_write(cas, sys)
Definition: casida.F90:4284
logical function isnt_degenerate(cas, st, ia, jb)
Definition: casida.F90:1260
subroutine, public casida_run(system, from_scratch)
Definition: casida.F90:321
subroutine doscillator_strengths(cas, mesh, st)
Definition: casida.F90:1390
subroutine dcasida_write(cas, sys)
Definition: casida.F90:2679
real(real64) function casida_matrix_factor(cas, sys)
Definition: casida.F90:1181
subroutine casida_run_legacy(sys, fromScratch)
Definition: casida.F90:339
integer, parameter casida_tamm_dancoff
Definition: casida.F90:190
subroutine zcasida_forces(cas, sys, gr, st)
Definition: casida.F90:3886
subroutine casida_type_end(cas)
Definition: casida.F90:915
subroutine dcasida_solve(cas, sys)
Definition: casida.F90:2542
This module implements a calculator for the density and defines related functions.
Definition: density.F90:120
subroutine, public states_elec_total_density(st, mesh, total_rho)
This routine calculates the total electronic density.
Definition: density.F90:851
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_epsilon
Definition: global.F90:204
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public io_function_read_what_how_when(namespace, space, what, how, output_interval, what_tag_in, how_tag_in, output_interval_tag_in, ignore_error)
Definition: io.F90:114
integer pure function, public kpoints_number(this)
Definition: kpoints.F90:1099
A module to handle KS potential, without the external potential.
integer, parameter, public hartree_fock
integer, parameter, public generalized_kohn_sham_dft
integer, parameter, public dft_u_none
Definition: lda_u.F90:201
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1096
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:530
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1028
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:696
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1068
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:599
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
Definition: mpi.F90:370
type(mpi_comm), parameter, public mpi_comm_undefined
used to indicate a communicator has not been initialized
Definition: mpi.F90:136
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:270
subroutine mpi_grp_init(grp, comm)
Initialize MPI group instance.
Definition: mpi.F90:308
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
logical pure function, public multicomm_strategy_is_parallel(mc, level)
Definition: multicomm.F90:724
This module implements the basic mulsisystem class, a container system for other systems.
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
subroutine, public photon_mode_set_n_electrons(this, qtot)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
integer, parameter, public restart_casida
Definition: restart.F90:200
subroutine, public restart_rm(restart, name)
Remove directory or file "name" that is located inside the current restart directory.
Definition: restart.F90:840
integer, parameter, public restart_gs
Definition: restart.F90:200
subroutine, public restart_init(restart, namespace, data_type, type, mc, ierr, mesh, dir, exact)
Initializes a restart object.
Definition: restart.F90:517
integer, parameter, public restart_type_dump
Definition: restart.F90:246
integer, parameter, public restart_type_load
Definition: restart.F90:246
subroutine, public restart_end(restart)
Definition: restart.F90:723
This module contains interfaces for ScaLAPACK routines Interfaces are from http:
Definition: scalapack.F90:131
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:117
pure logical function, public states_are_complex(st)
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_count_pairs(st, namespace, n_pairs, n_occ, n_unocc, is_included, is_frac_occ)
number of occupied-unoccupied pairs for Casida
This module handles reading and writing restart information for the states_elec_t.
subroutine, public states_elec_look_and_load(restart, namespace, space, st, mesh, kpoints, is_complex, packed)
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
This module defines the unit system, used for input and output.
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:118
subroutine, public v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, calc_eigenval, calc_current)
Definition: v_ks.F90:684
Definition: xc.F90:114
logical pure function, public family_is_mgga_with_exc(xcs)
Is the xc function part of the mGGA family with an energy functional.
Definition: xc.F90:584
subroutine, public xc_get_fxc(xcs, gr, namespace, rho, ispin, fxc, fxc_grad, fxc_grad_spin)
Returns the exchange-correlation kernel.
Definition: xc.F90:1975
pure logical function, public in_family(family, xc_families)
Definition: xc.F90:612
integer, parameter, public sic_adsic
Averaged density SIC.
Definition: xc_sic.F90:148
This class contains all parameters, needed for Casida calculations.
Definition: casida.F90:202
Class describing the electron system.
Definition: electrons.F90:218
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:169
Container class for lists of system_oct_m::system_t.
The states_elec_t class contains all electronic wave functions.
int true(void)