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