PROGRESS  master
prg_response_mod.F90
Go to the documentation of this file.
1 
7 module prg_response_mod
8 
9  use bml
10  use prg_kernelparser_mod
11  use prg_densitymatrix_mod
12  use prg_openfiles_mod
13 
14  implicit none
15 
16  private
17 
18  integer, parameter :: dp = kind(1.0d0)
19  real(dp), parameter :: pi = 3.14159265358979323846264338327950_dp
20 
21  type, public :: respdata_type
22  character(20) :: respmode
23  character(20) :: typeofpert
24  character(20) :: bmltype
25  integer :: mdim
26  real(dp) :: numthresh
27  logical :: computedipole
28  logical :: getresponse
29  real(dp) :: fieldintensity
30  real(dp) :: field(3)
31  end type respdata_type
32 
33  public :: prg_pert_from_file, prg_compute_dipole, prg_pert_constant_field
34  public :: prg_compute_response_rs, prg_parse_response, prg_compute_response_sp2
35  public :: prg_compute_response_fd, prg_compute_polarizability, prg_write_dipole_tcl
36  public :: prg_project_response, prg_pert_sin_pot, prg_pert_cos_pot, prg_canon_response, prg_canon_response_orig
37  public :: prg_canon_response_vector, prg_canon_response_p1_dpdmu
38 
39 contains
40 
45  subroutine prg_parse_response(RespData, filename)
46  implicit none
47  type(respdata_type) :: respdata
48  integer, parameter :: nkey_char = 3, nkey_int = 4, nkey_re = 5, nkey_log = 2
49  character(len=*) :: filename
50 
51  !Library of keywords with the respective defaults.
52  character(len=50), parameter :: keyvector_char(nkey_char) = [character(len=100) :: &
53  'TypeOfPert=', 'BMLType=','RespMode=']
54  character(len=100) :: valvector_char(nkey_char) = [character(len=100) :: &
55  'None', 'Dense','RayleighSchrodinger']
56 
57  character(len=50), parameter :: keyvector_int(nkey_int) = [character(len=50) :: 'Mdim=' &
58  , 'p2=', 'p3=', 'p4=']
59  integer :: valvector_int(nkey_int) = (/0 &
60  , 0, 0, 1 /)
61 
62  character(len=50), parameter :: keyvector_re(nkey_re) = [character(len=50) :: &
63  'Fieldx=','Fieldy=','Fieldz=' &
64  ,'FieldIntensity=','NumThresh=' ]
65  real(dp) :: valvector_re(nkey_re) = (/0.0,0.0,0.0&
66  ,0.001,0.0 /)
67 
68  character(len=50), parameter :: keyvector_log(nkey_log) = [character(len=100) :: &
69  'ComputeDipole=', 'GetResponse=']
70  logical :: valvector_log(nkey_log) = (/ .false., .false./)
71  !End of the library
72 
73  !Start and stop characters
74  character(len=50), parameter :: startstop(2) = [character(len=50) :: &
75  'RESPONSE{', '}']
76 
77  call prg_parsing_kernel(keyvector_char,valvector_char&
78  ,keyvector_int,valvector_int,keyvector_re,valvector_re,&
79  keyvector_log,valvector_log,trim(filename),startstop)
80 
81  !Characters
82  respdata%TypeofPert = valvector_char(1)
83 
84  if(valvector_char(2) == "Dense")then
85  respdata%BmlType = bml_matrix_dense
86  elseif(valvector_char(2) == "Ellpack")then
87  respdata%BmlType = bml_matrix_ellpack
88  elseif(valvector_char(2) == "Ellblock")then
89  respdata%BmlType = bml_matrix_ellblock
90  endif
91 
92  respdata%RespMode = valvector_char(3)
93 
94  !Reals
95  respdata%Field(1) = valvector_re(1)
96  respdata%Field(2) = valvector_re(2)
97  respdata%Field(3) = valvector_re(3)
98  respdata%FieldIntensity = valvector_re(4)
99  respdata%NumThresh = valvector_re(5)
100 
101  !Logicals
102  respdata%ComputeDipole = valvector_log(1)
103  respdata%GetResponse = valvector_log(2)
104 
105  !Integers
106  respdata%Mdim = valvector_int(1)
107 
108  end subroutine prg_parse_response
109 
122  subroutine prg_compute_dipole(charges,coordinate,dipoleMoment,factor,verbose)
123 
124  integer :: cont, i, nats, verbose
125  real(dp), intent(in) :: charges(:), coordinate(:,:), factor
126  real(dp), intent(inout) :: dipolemoment(3)
127 
128  if(verbose.gt.1)write(*,*)"Computing dipole moment ..."
129 
130  cont=0
131 
132  nats = size(charges,dim=1)
133 
134  dipolemoment = 0.0_dp
135 
136  do i=1,nats
137  dipolemoment(1)=dipolemoment(1)+coordinate(1,i)*charges(i)
138  dipolemoment(2)=dipolemoment(2)+coordinate(2,i)*charges(i)
139  dipolemoment(3)=dipolemoment(3)+coordinate(3,i)*charges(i)
140  enddo
141 
142  dipolemoment=factor*dipolemoment
143 
144  write(*,*)"dipoleMomentx=",dipolemoment(1)
145  write(*,*)"dipoleMomenty=",dipolemoment(2)
146  write(*,*)"dipoleMomentz=",dipolemoment(3)
147 
148  end subroutine prg_compute_dipole
149 
160  subroutine prg_write_dipole_tcl(dipoleMoment,file,factor,verbose)
161  implicit none
162  integer :: verbose, io_unit
163  real(dp), intent(in) :: factor
164  real(dp), intent(in) :: dipolemoment(3)
165  character(*), intent(in) :: file
166 
167  call prg_open_file(io_unit, "dipole.tcl")
168 
169  write(io_unit,*)"mol new ",file
170  write(io_unit,*)"set molid 0"
171 
172  write(io_unit,*)"proc vmd_draw_arrow {mol start end} {"
173  write(io_unit,*)" set middle [vecadd $start [vecscale 0.9 [vecsub $end $start]]]"
174  write(io_unit,*)" graphics $mol cylinder $start $middle radius 0.15"
175  write(io_unit,*)" graphics $mol cone $middle $end radius 0.25"
176  write(io_unit,*)"}"
177 
178  write(io_unit,*)"draw arrow {0 0 0} {",&
179  factor*dipolemoment(1), factor*dipolemoment(2), factor*dipolemoment(3),"}"
180 
181  close(io_unit)
182 
183  end subroutine prg_write_dipole_tcl
184 
200  subroutine prg_compute_polarizability(rsp_bml,prt_bml,polarizability,factor,verbose)
201  implicit none
202  integer :: cont, i, nats, verbose
203  real(dp), intent(in) :: factor
204  real(dp), intent(inout) :: polarizability
205  real(dp) :: dipolez, dipolex, dipoley
206  type(bml_matrix_t), intent(in) :: prt_bml, rsp_bml
207  type(bml_matrix_t) :: aux_bml
208 
209  if(verbose.gt.1)write(*,*)"Computing polarizability ..."
210 
211  call bml_copy_new(prt_bml,aux_bml)
212  call bml_multiply(rsp_bml,prt_bml,aux_bml,1.0_dp,0.0_dp,0.0_dp)
213 
214  polarizability = -factor*bml_trace(aux_bml)
215 
216  write(*,*)"Polarizability=",polarizability
217 
218  call bml_deallocate(aux_bml)
219 
220  end subroutine prg_compute_polarizability
221 
225  subroutine prg_pert_from_file(prt_bml,norb)
226  implicit none
227  integer :: norb, verbose
228  type(bml_matrix_t), intent(inout) :: prt_bml
229  if(verbose.gt.1) write(*,*)'Reading perturbation V from file ...'
230 
231  end subroutine prg_pert_from_file
232 
250  subroutine prg_compute_response_rs(ham_bml,prt_bml,rsp_bml,lambda&
251  ,bndfil,threshold,verbose)
252 
253  character(20) :: bml_type
254  integer :: i, j, l, norb
255  integer :: verbose
256  real(dp) :: lambda, nocc
257  real(dp), allocatable :: aux(:,:), evals(:), row(:)
258  real(dp), intent(in) :: bndfil, threshold
259  type(bml_matrix_t) :: aux1_bml, aux_bml, occupation_bml, umat_bml
260  type(bml_matrix_t) :: aux2_bml
261  type(bml_matrix_t), intent(in) :: ham_bml, prt_bml
262  type(bml_matrix_t), intent(inout) :: rsp_bml
263 
264  norb = bml_get_n(ham_bml)
265 
266  if(verbose.ge.1) write(*,*)'In prg_compute_response_RS... Norbs = ', norb
267 
268  allocate(evals(norb))
269 
270  bml_type = bml_get_type(ham_bml)
271 
272  !Ensure dense type to diagonalize
273  allocate(aux(norb,norb))
274  call bml_export_to_dense(ham_bml,aux)
275  call bml_zero_matrix(bml_matrix_dense,bml_element_real,dp,norb,norb,aux_bml)
276  call bml_import_from_dense(bml_matrix_dense,aux,aux_bml,threshold,norb)
277  deallocate(aux)
278 
279  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,umat_bml)
280 
281  !Eigenvectors and eigevalues of H
282  call bml_diagonalize(aux_bml,evals,umat_bml)
283  call bml_deallocate(aux_bml)
284 
285  !Ensure the umat type is in bml_type.
286  allocate(aux(norb,norb))
287  call bml_export_to_dense(umat_bml,aux)
288  call bml_deallocate(umat_bml)
289  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,umat_bml)
290  call bml_import_from_dense(bml_type,aux,umat_bml,threshold,norb)
291  deallocate(aux)
292 
293  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux_bml)
294  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux1_bml)
295 
296  !Obtaining matrix $\f V = C^{\dagger} H^{(1)} C^{\dagger} $\f
297  call bml_transpose(umat_bml,aux_bml) !C^{\dagger}
298  !First product C^{\dagger} H^{(1)}
299  call bml_multiply(aux_bml,prt_bml,aux1_bml,1.0_dp,0.0_dp,threshold)
300  !Seccond product to obtain V (aux_bml)
301  call bml_multiply(aux1_bml,umat_bml,aux_bml,1.0_dp,0.0_dp,threshold)
302 
303  !\f$ \tilde{V}_{ij} = \frac{V_{ij}}{\epsilon_j - \epsilon_i} \f$
304  allocate(row(norb))
305  do i=1,norb
306  call bml_get_row(aux_bml,i,row)
307  do j=1,norb
308  if(j.ne.i)then
309  row(j) = row(j)/(evals(j)-evals(i))
310  else
311  row(j) = 0.0_dp
312  endif
313  enddo
314  call bml_set_row(aux_bml,i,row,threshold)
315  enddo
316  deallocate(row)
317 
318  !\f$ C^{(1)} = C \tilde{V} \f$
319  call bml_multiply(umat_bml,aux_bml,aux1_bml,1.0_dp,0.0_dp,threshold)
320 
321  !Get matrix f
322  nocc = norb*bndfil
323  write(*,*)"Trace CV",bml_trace(aux_bml), nocc
324 
325  do i=1,norb !Reusing eigenvalues to apply the theta function.
326  if(i.le.nocc) then
327  evals(i) = 2.0_dp
328  else
329  evals(i) = 0.0_dp
330  endif
331  enddo
332  if(abs(nocc - int(nocc)).gt.0.01_dp)then
333  evals(int(nocc)+1) = 1.0_dp
334  endif
335 
336  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,occupation_bml)
337  call bml_set_diagonal(occupation_bml, evals) !eps(i,i) = eps(i)
338  deallocate(evals)
339 
340  !Cf(C^{(1)})^{t}$
341  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,aux2_bml)
342  call bml_multiply(umat_bml,occupation_bml,aux_bml,1.0_dp,0.0_dp,threshold)
343  call bml_transpose(aux1_bml,aux2_bml)
344  call bml_multiply(aux_bml,aux2_bml,rsp_bml,1.0_dp,0.0_dp,threshold)
345 
346  !CfC^{(1)}$ + C^{(1)}fC
347  call bml_transpose(umat_bml,aux2_bml)
348  call bml_multiply(occupation_bml,aux2_bml,aux_bml,1.0_dp,0.0_dp,threshold)
349 
350  call bml_deallocate(occupation_bml)
351  call bml_deallocate(umat_bml)
352 
353  call bml_multiply(aux1_bml,aux_bml,rsp_bml,1.0_dp,1.0_dp,threshold)
354  call bml_scale(1.0_dp,rsp_bml)
355 
356  write(*,*)"Trace",bml_trace(rsp_bml)
357 
358  call bml_deallocate(aux_bml)
359  call bml_deallocate(aux1_bml)
360 
361  end subroutine prg_compute_response_rs
362 
380  subroutine prg_compute_response_fd(ham_bml,prt_bml,rsp_bml,prg_delta&
381  ,bndfil,threshold,verbose)
382 
383  character(20) :: bml_type
384  integer :: i, j, l, norb
385  integer :: verbose
386  real(dp) :: prg_delta, nocc
387  real(dp), allocatable :: resp(:,:), evals(:), row(:)
388  real(dp), intent(in) :: bndfil, threshold
389  type(bml_matrix_t) :: rhoplus_bml, rhominus_bml, occupation_bml, umat_bml
390  type(bml_matrix_t) :: aux_bml
391  type(bml_matrix_t), intent(in) :: ham_bml, prt_bml
392  type(bml_matrix_t), intent(inout) :: rsp_bml
393 
394  norb = bml_get_n(ham_bml)
395 
396  if(verbose.ge.1) write(*,*)'In prg_compute_response_RS... Norbs = ', norb
397 
398  allocate(resp(norb,norb))
399 
400  call bml_copy_new(ham_bml,aux_bml)
401 
402  bml_type = bml_get_type(ham_bml)
403 
404  !! - \f$ H^+ = H^{(1)} + \delta I \f$
405  call bml_add_deprecated(1.0_dp, aux_bml, prg_delta, prt_bml, threshold)
406 
407  call bml_copy_new(ham_bml,rhoplus_bml)
408 
409  !! - \f$ \rho^+ = f(H + H^+) \f$
410  call prg_build_density_t0(aux_bml,rhoplus_bml,threshold,bndfil)
411 
412  !! - \f$ H^- = H^{(1)} - \delta I \f$
413  call bml_add_deprecated(1.0_dp, aux_bml, -2.0_dp*prg_delta, prt_bml, threshold)
414 
415  call bml_copy_new(ham_bml,rhominus_bml)
416 
417  !! - \f$ \rho^- = f(H + H^-) \f$
418  call prg_build_density_t0(aux_bml,rhominus_bml,threshold,bndfil)
419 
420  !! - \f$ \rho^{(1)} = (\rho^+ - \rho^-)/(2\delta) \f$.
421  call bml_add_deprecated(1.0_dp, rhoplus_bml, -1.0_dp, rhominus_bml, threshold)
422  call bml_scale(1.0_dp/(2.0_dp*prg_delta),rhoplus_bml,rsp_bml)
423 
424  end subroutine prg_compute_response_fd
425 
447  subroutine prg_pert_constant_field(field,intensity,coordinate,lambda&
448  ,prt_bml,threshold,spindex,norbi,verbose,over_bml)
449 
450  integer :: cont, i, ii, nats
451  integer :: norb
452  integer, intent(in) :: spindex(:), norbi(:), verbose
453  real(dp) :: fieldnorm, fieldx, fieldy, fieldz
454  real(dp) :: intensity, lambda, threshold
455  real(dp), allocatable :: diag(:)
456  real(dp), intent(in) :: coordinate(:,:), field(3)
457  type(bml_matrix_t) :: aux_bml
458  type(bml_matrix_t), intent(inout) :: prt_bml
459  type(bml_matrix_t), optional, intent(in) :: over_bml
460 
461  write(*,*)'Applying a constant field perturbation ...'
462 
463  fieldx = field(1)
464  fieldy = field(2)
465  fieldz = field(3)
466 
467  nats=size(spindex,dim=1)
468  norb=sum(norbi(spindex(:)))
469 
470  write(*,*)nats,norb
471 
472  fieldnorm=sqrt(fieldx**2+fieldy**2+fieldz**2) !Get the norm
473 
474  fieldx=intensity*fieldx/fieldnorm !Normalize and add intensity
475  fieldy=intensity*fieldy/fieldnorm
476  fieldz=intensity*fieldz/fieldnorm
477 
478  cont=0;
479 
480  allocate(diag(norb))
481 
482  do i=1,nats
483  do ii=1,norbi(spindex(i))
484  cont=cont+1
485  diag(cont)=lambda*(fieldx*coordinate(1,i) + &
486  fieldy*coordinate(2,i) + fieldz*coordinate(3,i))
487  enddo
488  enddo
489 
490  if(bml_get_n(prt_bml) < 0) stop "bml_pert not allocated"
491 
492  call bml_set_diagonal(prt_bml,diag,threshold)
493 
494  deallocate(diag)
495 
496  if(present(over_bml))then !(S*V + V*S)/2
497  call bml_copy_new(prt_bml,aux_bml)
498  call bml_multiply(over_bml,aux_bml,prt_bml,0.5_dp, 0.0_dp,threshold)
499  call bml_multiply(aux_bml,over_bml,prt_bml,0.5_dp,1.0_dp,threshold)
500  call bml_deallocate(aux_bml)
501  endif
502 
503  end subroutine prg_pert_constant_field
504 
524  subroutine prg_pert_sin_pot(direction,lx,coordinate,lambda&
525  ,prt_bml,threshold,spindex,norbi,verbose,over_bml)
526 
527  integer :: cont, i, ii, nats
528  integer :: norb
529  integer, intent(in) :: spindex(:), norbi(:), verbose
530  real(dp) :: fieldnorm, fieldx, fieldy, fieldz
531  real(dp) :: lx
532  real(dp) :: intensity, lambda, threshold
533  real(dp), allocatable :: diag(:)
534  real(dp), intent(in) :: coordinate(:,:)
535  type(bml_matrix_t) :: aux_bml
536  type(bml_matrix_t), intent(inout) :: prt_bml
537  type(bml_matrix_t), optional, intent(in) :: over_bml
538  character :: direction
539 
540  write(*,*)'Applying a constant field perturbation ...'
541 
542  nats=size(spindex,dim=1)
543  norb=sum(norbi(spindex(:)))
544 
545  write(*,*)nats,norb
546 
547  cont=0;
548 
549  allocate(diag(norb))
550 
551  do i=1,nats
552  do ii=1,norbi(spindex(i))
553  cont=cont+1
554  diag(cont)=lambda*sin(2.0_dp*pi*((coordinate(1,i)-lx)/lx)-pi)
555  enddo
556  enddo
557 
558  if(bml_get_n(prt_bml) < 0) stop "bml_pert not allocated"
559 
560  call bml_set_diagonal(prt_bml,diag,threshold)
561 
562  deallocate(diag)
563 
564  if(present(over_bml))then !(S*V + V*S)/2
565  call bml_copy_new(prt_bml,aux_bml)
566  call bml_multiply(over_bml,aux_bml,prt_bml,0.5_dp, 0.0_dp,threshold)
567  call bml_multiply(aux_bml,over_bml,prt_bml,0.5_dp,1.0_dp,threshold)
568  call bml_deallocate(aux_bml)
569  endif
570 
571  end subroutine prg_pert_sin_pot
572 
592  subroutine prg_pert_cos_pot(direction,lx,coordinate,lambda&
593  ,prt_bml,threshold,spindex,norbi,verbose,over_bml)
594 
595  integer :: cont, i, ii, nats
596  integer :: norb
597  integer, intent(in) :: spindex(:), norbi(:), verbose
598  real(dp) :: fieldnorm, fieldx, fieldy, fieldz
599  real(dp) :: lx
600  real(dp) :: intensity, lambda, threshold
601  real(dp), allocatable :: diag(:)
602  real(dp), intent(in) :: coordinate(:,:)
603  type(bml_matrix_t) :: aux_bml
604  type(bml_matrix_t), intent(inout) :: prt_bml
605  type(bml_matrix_t), optional, intent(in) :: over_bml
606  character :: direction
607 
608  write(*,*)'Applying a constant field perturbation ...'
609 
610  nats=size(spindex,dim=1)
611  norb=sum(norbi(spindex(:)))
612 
613  write(*,*)nats,norb
614 
615  cont=0;
616 
617  allocate(diag(norb))
618 
619  do i=1,nats
620  do ii=1,norbi(spindex(i))
621  cont=cont+1
622  diag(cont)=lambda*cos(2.0_dp*pi*((coordinate(1,i)-lx)/lx)-pi)
623  enddo
624  enddo
625 
626  if(bml_get_n(prt_bml) < 0) stop "bml_pert not allocated"
627 
628  call bml_set_diagonal(prt_bml,diag,threshold)
629 
630  deallocate(diag)
631 
632  if(present(over_bml))then !(S*V + V*S)/2
633  call bml_copy_new(prt_bml,aux_bml)
634  call bml_multiply(over_bml,aux_bml,prt_bml,0.5_dp, 0.0_dp,threshold)
635  call bml_multiply(aux_bml,over_bml,prt_bml,0.5_dp,1.0_dp,threshold)
636  call bml_deallocate(aux_bml)
637  endif
638 
639  end subroutine prg_pert_cos_pot
640 
642  ! and a perturbation V by purification. The method implemented here is the SP2 method
643  ! [A. Niklasson, Phys. Rev. B, 66, 155115 (2002)].
644  !! \param ham_bml Hamiltonian in bml format (\f$ H^{(0)} \f$).
645  !! \param prt_bml Perturbation in bml format (\f$ H^{(1)} \f$).
646  !! \param rsp_bml First order response to the perturbation (\f$ \rho^{(1)} \f$).
647  !! \param bndfil Filing factor.
648  !! \param threshold Threshold value for matrix elements.
649  !! \param verbose Different levels of verbosity.
650  !! \warning This works only for the prg_orthogonalized form of ham_bml.
651  !! \warning The response must be in the prg_orthogonalized form.
652  !!
653  subroutine prg_compute_response_sp2(ham_bml,prt_bml,rsp_bml,rho_bml,lambda&
654  ,bndfil,minsp2iter,maxsp2iter,sp2conv,idemtol,threshold,verbose)
655 
656  character(20) :: bml_type
657  integer :: i, j, l, norb
658  integer :: verbose, mdim
659  integer, intent(in) :: minsp2iter, maxsp2iter
660  real(dp) :: lambda, occ, maxminusmin
661  real(dp) :: alpha, beta, trx, trd
662  real(dp), allocatable :: aux(:,:), evals(:), row(:), gbnd(:)
663  real(dp), intent(in) :: bndfil, threshold,idemtol
664  character(len=*), intent(in) :: sp2conv
665  type(bml_matrix_t) :: aux1_bml, aux_bml, occupation_bml
666  type(bml_matrix_t) :: x2_bml
667  type(bml_matrix_t), intent(in) :: ham_bml, prt_bml
668  type(bml_matrix_t), intent(inout) :: rsp_bml, rho_bml
669  real(dp) :: idemperr, idemperr1, idemperr2
670  integer :: iter, breakloop
671  real(dp) :: trx2, trxold, tr2xx2
672  real(dp) :: limdiff
673  real(dp), allocatable :: trace(:)
674 
675  if(verbose.ge.1) write(*,*)'In prg_compute_response_SP2... '
676 
677  norb = bml_get_n(ham_bml)
678  bml_type = bml_get_type(ham_bml)
679  mdim = bml_get_m(ham_bml)
680 
681  idemperr = 0.0_dp
682  idemperr1 = 0.0_dp
683  idemperr2 = 0.0_dp
684 
685  if(bml_get_n(rsp_bml).le.0)then
686  call bml_zero_matrix(bml_matrix_dense,bml_element_real,dp,norb,norb,rsp_bml)
687  endif
688 
689  occ = bndfil*float(norb)
690 
691  call bml_copy(ham_bml, rho_bml)
692  call bml_copy(prt_bml,rsp_bml)
693 
694  allocate(gbnd(2))
695  call bml_gershgorin(ham_bml, gbnd)
696  maxminusmin = gbnd(2) - gbnd(1)
697  alpha = -1.00_dp/maxminusmin
698  beta = gbnd(2)/maxminusmin
699  call bml_scale_add_identity(rho_bml, alpha, beta, 0.00_dp)
700  call bml_scale(alpha, rsp_bml) !Same as SP2 but with the perturbation.
701 
702  allocate(trace(2))
703  trx = bml_trace(rho_bml)
704  trd = bml_trace(rsp_bml)
705 
706  iter = 0
707  breakloop = 0
708 
709  ! X2 <- X
710  call bml_copy_new(rho_bml, x2_bml)
711  call bml_copy_new(rsp_bml, aux_bml)
712 
713  do while (breakloop .eq. 0 .and. iter .lt. maxsp2iter)
714  iter = iter + 1
715 
716  call bml_print_matrix("rsp_bml",rsp_bml,0,1,0,1)
717 
718  ! X2 <- X * X
719  call bml_multiply_x2(rho_bml, x2_bml, threshold, trace)
720  trd = bml_trace(rsp_bml)
721 
722  !Compute anticonmutator {X_n^0,Delta_n}
723  call bml_multiply(rho_bml, rsp_bml, aux_bml, 1.0d0, 0.0d0)
724  call bml_multiply(rsp_bml,rho_bml, aux_bml, 1.0d0, 1.0d0)
725 
726  trx2 = trace(2)
727  write(*,*) 'iter = ', iter, 'trx = ', trx, ' trx2 = ', trx2
728  write(*,*) 'iter = ', iter, 'tr(resp) = ', trd
729 
730  tr2xx2 = 2.0_dp*trx - trx2
731  trxold = trx
732  limdiff = abs(trx2 - occ) - abs(tr2xx2 - occ)
733 
734  if (limdiff .ge. idemtol) then
735 
736  ! X <- 2 * X - X2
737  call bml_add_deprecated(2.0_dp, rho_bml, -1.0_dp, x2_bml, threshold)
738  call bml_add_deprecated(2.0_dp,rsp_bml,-1.0_dp, aux_bml, threshold)
739 
740  trx = 2.0_dp * trx - trx2
741 
742  elseif(limdiff .lt. -idemtol) then
743 
744  ! X <- X2
745  call bml_copy(x2_bml, rho_bml)
746  call bml_copy(aux_bml, rsp_bml)
747 
748  trx = trx2
749 
750  else
751 
752  trx = trxold
753  breakloop = 1
754 
755  end if
756 
757  idemperr2 = idemperr1
758  idemperr1 = idemperr
759  idemperr = abs(trx - trxold)
760 
761  if (sp2conv .eq. "Rel" .and. iter .ge. minsp2iter .and. &
762  (idemperr .ge. idemperr2 .or. idemperr .lt. idemtol)) then
763  breakloop = 1
764  end if
765 
766  if (iter .eq. maxsp2iter) then
767  write(*,*) "SP2 purification is not converging: STOP!"
768  stop
769  end if
770 
771  enddo
772 
773  ! X <- 2 * X
774  call bml_scale(2.0_dp, rho_bml) !D = 2.0_dp*D
775  call bml_scale(2.0_dp, rsp_bml)
776 
777  call bml_deallocate(x2_bml)
778  call bml_deallocate(aux_bml)
779 
780  deallocate(trace)
781 
782  end subroutine prg_compute_response_sp2
783 
797  subroutine prg_project_response(rsp_bml,over_bml,spindex,norbi,coordinates,rspfunc,verbose)
798 
799  integer :: cont, i, j, nats
800  integer :: norbs
801  integer, intent(in) :: norbi(:), spindex(:), verbose
802  real(dp), allocatable :: diagonal(:)
803  real(dp), allocatable, intent(inout) :: rspfunc(:)
804  real(dp), intent(in) :: coordinates(:,:)
805  type(bml_matrix_t) :: aux_bml
806  type(bml_matrix_t), intent(in) :: over_bml
807  type(bml_matrix_t), intent(inout) :: rsp_bml
808 
809  nats=size(spindex, dim=1)
810  norbs=bml_get_n(rsp_bml)
811  call bml_copy_new(rsp_bml,aux_bml)
812  call bml_multiply(rsp_bml,over_bml,aux_bml,1.0_dp,0.0_dp,0.0_dp)
813 
814  if(.not. allocated(rspfunc))allocate(rspfunc(nats))
815  allocate(diagonal(norbs))
816  diagonal=0.0_dp
817  call bml_get_diagonal(aux_bml,diagonal)
818  ! call bml_get_diagonal(rsp_bml,diagonal)
819 
820  cont=0
821  rspfunc = 0.0_dp
822  do i=1,nats
823  do j=1,norbi(spindex(i))
824  cont=cont+1
825  rspfunc(i) = rspfunc(i) + diagonal(cont)
826  enddo
827  enddo
828 
829  deallocate(diagonal)
830 
831  call bml_deallocate(aux_bml)
832 
833  end subroutine prg_project_response
834 
835 
849  subroutine prg_canon_response_vector(P1_bml,H1_bml,Nocc,beta,evals,mu0,m,thresh,HDIM)
850 
851  real(dp), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
852  type(bml_matrix_t), intent(inout) :: p1_bml, h1_bml
853  type(bml_matrix_t) :: x_bml
854  integer, intent(in) :: hdim, m ! Nocc = Number of occupied orbitals, m= Number of recursion steps
855  real(dp), intent(in) :: evals(hdim), thresh !Q and e are eigenvectors and eigenvalues of H0
856  real(dp), intent(in) :: beta, mu0 ! Electronic temperature and chemical potential
857  real(dp), intent(in) :: nocc
858  real(dp) :: x(hdim,hdim), dx1(hdim,hdim), y(hdim,hdim) !Temporary matrices
859  real(dp) :: h_0(hdim), p_0(hdim), dpdmu(hdim), p_02(hdim),id0(hdim)
860  real(dp) :: cnst, kb, mu1
861  real(dp), allocatable :: focc(:), row1(:)
862  real(dp) :: tmp, tmp2, dmu
863  integer :: i, j, k
864  character(20) :: bml_type
865 
866  kb = 8.61739e-5 ! (eV/K)
867 
868  allocate(focc(hdim))
869  do i = 1, hdim
870  focc(i) = 1.0/(1+exp(beta*(evals(i) - mu0)))
871  enddo
872  do i = 1, hdim
873  do j = 1, hdim
874  x(j, i) = evals(i) - evals(j)
875  y(j, i) = focc(i) - focc(j)
876  enddo
877  enddo
878  deallocate(focc)
879 
880  do i = 1, hdim
881  tmp = exp((evals(i) - mu0) * beta)
882  tmp2 = (1.d0 + tmp) * (1.d0 + tmp)
883  dpdmu(i) = -tmp/tmp2 * beta
884  enddo
885 
886  do i = 1, hdim
887  x(i,i) = x(i,i) + 1.d0
888  enddo
889 
890  do i = 1, hdim
891  do j = 1, hdim
892  x(j,i) = y(j,i) / x(j,i)
893  enddo
894  x(i,i) = x(i,i) + dpdmu(i)
895  enddo
896 
897  ! H1_bml in this subroutine is already in transformed representation
898  allocate(row1(hdim))
899  bml_type = bml_get_type(h1_bml)
900 
901  ! element_wise multiplication (in bml)
902  call bml_zero_matrix(bml_type, bml_element_real, dp, hdim, hdim, x_bml)
903 
904  call bml_import_from_dense(bml_type, x, x_bml, thresh,hdim)
905  call bml_element_multiply_ab(x_bml, h1_bml, p1_bml, thresh)
906  call bml_get_diagonal(p1_bml,row1)
907 
908  tmp = 0.0
909  tmp2 = 0.0
910  do i = 1, hdim
911  tmp = tmp + x(i,i)
912  tmp2 = tmp2 + row1(i)
913  enddo
914 
915  if (abs(tmp) > 1.d-12) then
916  dmu = tmp2 / tmp
917  else
918  dmu = 0.d0
919  endif
920 
921  do i = 1, hdim
922  row1(i) = row1(i) - dmu * x(i,i)
923  enddo
924 
925  !set Y diag
926  call bml_set_diagonal(p1_bml, row1)
927  call bml_deallocate(x_bml)
928 
929  deallocate(row1)
930 
931  end subroutine prg_canon_response_vector
932 
933 
947  subroutine prg_canon_response(P1_bml,H1_bml,Nocc,beta,evals,mu0,m,HDIM)
948 
949  type(bml_matrix_t), intent(inout) :: p1_bml, h1_bml
950  type(bml_matrix_t) :: dx1_bml
951  real(dp), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
952  integer, intent(in) :: hdim, m ! Nocc = Number of occupied orbitals, m= Number of recursion steps
953  real(dp), intent(in) :: evals(hdim) !Q and e are eigenvectors and eigenvalues of H0
954  real(dp), intent(in) :: beta, mu0 ! Electronic temperature and chemical potential
955  real(dp) :: x(hdim,hdim), dx1(hdim,hdim), y(hdim,hdim) !Temporary matrices
956  real(dp) :: h_0(hdim), p_0(hdim), dpdmu(hdim), p_02(hdim),id0(hdim)
957  real(dp) :: cnst, kb, mu1
958  real(dp), allocatable :: row1(:), row2(:)
959  real(dp), intent(in) :: nocc
960  integer :: i, j, k
961 
962  kb = 8.61739e-5 ! (eV/K)
963  h_0 = evals ! Diagonal Hamiltonian H0 respresented in the eigenbasis Q
964  cnst = beta/(1.d0*2**(m+2)) ! Scaling constant
965  p_0 = 0.5d0 + cnst*(h_0-mu0) ! Initialization for P0 represented in eigenbasis Q
966 
967  call bml_copy_new(h1_bml,p1_bml)
968  call bml_copy_new(h1_bml,dx1_bml)
969  call bml_scale(-cnst,p1_bml) !(set mu1 = 0 for simplicity) ! Initialization of DM response in Q representation (not diagonal in Q)
970 
971  allocate(row1(hdim))
972  allocate(row2(hdim))
973 
974  do i = 1,m ! Loop over m recursion steps
975  p_02 = p_0*p_0
976  !$omp parallel do default(none) private(k) &
977  !$omp private(j,row1,row2) &
978  !$omp shared(HDIM,P1_bml,p_0,DX1_bml)
979  do k = 1,hdim
980  call bml_get_row(p1_bml,k,row1)
981  do j = 1,hdim
982  row2(j) = p_0(k)*row1(j) + row1(j)*p_0(j)
983  enddo
984  call bml_set_row(dx1_bml,k,row2)
985  enddo
986  !$omp end parallel do
987  id0 = 1.d0/(2.d0*(p_02-p_0)+1.d0)
988  p_0 = id0*p_02
989 
990  !$omp parallel do default(none) private(k) &
991  !$omp private(j,row1,row2) &
992  !$omp shared(HDIM,P1_bml,p_0,DX1_bml,iD0)
993  do k = 1,hdim
994  call bml_get_row(p1_bml,k,row1)
995  call bml_get_row(dx1_bml,k,row2)
996  do j = 1,hdim
997  row1(j) = id0(k)*(row2(j) + 2.d0*(row1(j)-row2(j))*p_0(j))
998  enddo
999  call bml_set_row(p1_bml,k,row1)
1000  enddo
1001  !$omp end parallel do
1002 
1003  enddo
1004 
1005  deallocate(row2)
1006  call bml_deallocate(dx1_bml)
1007 
1008  dpdmu = beta*p_0*(1.d0-p_0)
1009  mu1 = 0.d0
1010  call bml_get_diagonal(p1_bml,row1)
1011  do i = 1,hdim
1012  mu1 = mu1 + row1(i)
1013  enddo
1014 
1015  mu1 = -mu1/sum(dpdmu)
1016  do i = 1,hdim
1017  row1(i) = row1(i) + mu1*dpdmu(i) ! Trace correction by adding (dP/dmu)*(dmu/dH1) to dP/dH1
1018  enddo
1019  call bml_set_diagonal(p1_bml,row1)
1020  deallocate(row1)
1021 
1022  end subroutine prg_canon_response
1023 
1024 
1038  subroutine prg_canon_response_orig(P1_bml,H1_bml,Nocc,beta,evals,mu0,m,thresh,HDIM)
1039 
1040  type(bml_matrix_t), intent(inout) :: p1_bml, h1_bml
1041  type(bml_matrix_t) :: dx1_bml
1042  real(dp), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
1043  integer, intent(in) :: hdim, m ! Nocc = Number of occupied orbitals, m=Number of recursion steps
1044  real(dp), intent(in) :: evals(hdim), thresh !Q and e are eigenvectors and eigenvalues of H0
1045  real(dp), intent(in) :: beta, mu0 ! Electronic temperature and chemical potential
1046  real(dp) :: x(hdim,hdim), dx1(hdim,hdim), y(hdim,hdim) !Temporary matrices
1047  real(dp) :: h_0(hdim), p_0(hdim), dpdmu(hdim),p_02(hdim),id0(hdim)
1048  real(dp) :: cnst, kb, mu1
1049  real(dp), allocatable :: p1(:,:), h1(:,:)
1050  real(dp), intent(in) :: nocc
1051  integer :: i, j, k
1052  character(20) :: bml_type
1053 
1054  kb = 8.61739e-5 ! (eV/K)
1055  h_0 = evals ! Diagonal Hamiltonian H0 respresented in the eigenbasis Q
1056  cnst = beta/(1.d0*2**(m+2)) ! Scaling constant
1057  p_0 = 0.5d0 + cnst*(h_0-mu0) ! Initialization for P0 represented in eigenbasis Q
1058  bml_type = bml_get_type(h1_bml)
1059  if(.not.allocated(h1))allocate(h1(hdim,hdim))
1060  call bml_export_to_dense(h1_bml,h1)
1061  if(.not.allocated(p1))allocate(p1(hdim,hdim))
1062  p1 = -cnst*h1
1063  dx1 = h1
1064 
1065  do i = 1,m ! Loop over m recursion steps
1066  p_02 = p_0*p_0
1067  !$omp parallel do default(none) private(k) &
1068  !$omp private(j) &
1069  !$omp shared(HDIM,P1,p_0,DX1)
1070  do k = 1,hdim
1071  do j = 1,hdim
1072  dx1(k,j) = p_0(k)*p1(k,j) + p1(k,j)*p_0(j)
1073  enddo
1074  enddo
1075  !$omp end parallel do
1076  id0 = 1.d0/(2.d0*(p_02-p_0)+1.d0)
1077  p_0 = id0*p_02
1078 
1079  !$omp parallel do default(none) private(k) &
1080  !$omp private(j) &
1081  !$omp shared(HDIM,P1,p_0,DX1,iD0)
1082  do k = 1,hdim
1083  do j = 1,hdim
1084  p1(k,j) = id0(k)*(dx1(k,j) + 2.d0*(p1(k,j)-dx1(k,j))*p_0(j))
1085  enddo
1086  enddo
1087  !$omp end parallel do
1088  enddo
1089 
1090  dpdmu = beta*p_0*(1.d0-p_0)
1091  mu1 = 0.d0
1092  do i = 1,hdim
1093  mu1 = mu1 + p1(i,i)
1094  enddo
1095 
1096  mu1 = -mu1/sum(dpdmu)
1097  do i = 1,hdim
1098  p1(i,i) = p1(i,i) + mu1*dpdmu(i) ! Trace correction by adding (dP/dmu)*(dmu/dH1) to dP/dH1
1099  enddo
1100 
1101  call bml_import_from_dense(bml_type, p1, p1_bml, thresh,hdim)
1102  deallocate(p1)
1103  deallocate(h1)
1104 
1105  end subroutine prg_canon_response_orig
1106 
1112  subroutine prg_canon_response_p1_dpdmu(P1_bml,dPdMu,H1_bml,Norbs,beta,Q_bml,evals,mu0,m,HDIM)
1113 
1114  type(bml_matrix_t), intent(inout) :: p1_bml, h1_bml
1115  type(bml_matrix_t) :: dx1_bml
1116  type(bml_matrix_t) :: q_bml
1117  real(dp), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0
1118  integer, intent(in) :: hdim, m, norbs ! = Number of occupied orbitals,m= Number of recursion steps
1119  real(dp), intent(in) :: evals(hdim) !Q and e are eigenvectors and eigenvalues of H0
1120  real(dp), intent(inout) :: dpdmu(hdim) !Q and e are eigenvectors and eigenvalues of H0
1121  ! real(PREC) :: P1(HDIM,HDIM) ! Density matrix response
1122  ! derivative with respect to perturbation H1 to H0
1123  real(dp), intent(in) :: beta, mu0 ! Electronic temperature and chemicalpotential
1124  real(dp), allocatable :: p1(:,:), dx1(:,:)
1125  real(dp) :: h_0(hdim), p_0(hdim),p_02(hdim),id0(hdim)
1126  real(dp) :: cnst, kb, mu1, sumdpdmu
1127  integer :: i, j, k
1128  character(20) :: bml_type
1129 
1130  kb = 8.61739e-5 ! (eV/K)
1131  h_0 = evals ! Diagonal Hamiltonian H0 respresented in the eigenbasis Q
1132  cnst = beta/(1.d0*2**(m+2)) ! Scaling constant
1133  p_0 = 0.5d0 + cnst*(h_0-mu0) ! Initialization for P0 represented in eigenbasis Q
1134  call bml_copy_new(h1_bml,p1_bml)
1135  call bml_copy_new(h1_bml,dx1_bml)
1136  call bml_scale(-cnst,p1_bml) !(set mu1 = 0 for simplicity) !Initialization of DM response in Q representation (not diagonal in Q)
1137 
1138  allocate(p1(hdim,hdim))
1139  allocate(dx1(hdim,hdim))
1140 
1141  p1 = 0.0_dp
1142  dx1 = 0.0_dp
1143 
1144  call bml_export_to_dense(p1_bml,p1)
1145  call bml_export_to_dense(dx1_bml,dx1)
1146 
1147 
1148  do i = 1,m ! Loop over m recursion steps
1149  p_02 = p_0*p_0
1150  !$omp parallel do default(none) private(k) &
1151  !$omp private(j) &
1152  !$omp shared(HDIM,P1,p_0,DX1)
1153  do k = 1,hdim
1154  dx1(k,:) = (p_0(k) + p_0(:))*p1(k,:)
1155  enddo
1156  !$omp end parallel do
1157  id0 = 1.d0/(2.d0*(p_02-p_0)+1.d0)
1158  p_0 = id0*p_02
1159 
1160  !$omp parallel do default(none) private(k) &
1161  !$omp private(j) &
1162  !$omp shared(HDIM,P1,p_0,DX1,iD0)
1163  do k = 1,hdim
1164  p1(k,:) = id0(k)*(dx1(k,:) + 2.d0*(p1(k,:)-dx1(k,:))*p_0(:))
1165  enddo
1166  !$omp end parallel do
1167  enddo
1168 
1169  bml_type = bml_get_type(p1_bml)
1170  call bml_import_from_dense(bml_type,p1,p1_bml,zero,hdim) !Dense to dense_bml
1171 
1172  deallocate(p1)
1173  deallocate(dx1)
1174 
1175  call bml_deallocate(dx1_bml)
1176 
1177  dpdmu = beta*p_0*(1.d0-p_0)
1178 
1179  end subroutine prg_canon_response_p1_dpdmu
1180 
1181 end module prg_response_mod
prg_densitymatrix_mod
Module to obtain the density matrix by diagonalizing an orthogonalized Hamiltonian.
Definition: prg_densitymatrix_mod.F90:5
prg_densitymatrix_mod::dp
integer, parameter dp
Definition: prg_densitymatrix_mod.F90:14
prg_densitymatrix_mod::prg_build_density_t0
subroutine, public prg_build_density_t0(ham_bml, rho_bml, threshold, bndfil, eigenvalues_out)
Builds the density matrix from for zero electronic temperature. Where, is the matrix eigenvector a...
Definition: prg_densitymatrix_mod.F90:37
prg_kernelparser_mod
Some general parsing functions.
Definition: prg_kernelparser_mod.F90:4
prg_kernelparser_mod::prg_parsing_kernel
subroutine, public prg_parsing_kernel(keyvector_char, valvector_char, keyvector_int, valvector_int, keyvector_re, valvector_re, keyvector_log, valvector_log, filename, startstop)
The general parsing function. It is used to vectorize a set of "keywords" "value" pairs as included i...
Definition: prg_kernelparser_mod.F90:30
prg_openfiles_mod
Module to handle input output files for the PROGRESS lib.
Definition: prg_openfiles_mod.F90:4
prg_openfiles_mod::prg_open_file
subroutine, public prg_open_file(io, name)
Opens a file to write.
Definition: prg_openfiles_mod.F90:38
prg_response_mod
Module to compute the density matrix response and related quantities.
Definition: prg_response_mod.F90:7
prg_response_mod::prg_canon_response_orig
subroutine, public prg_canon_response_orig(P1_bml, H1_bml, Nocc, beta, evals, mu0, m, thresh, HDIM)
First-order Canonical Density Matrix Perturbation Theory.
Definition: prg_response_mod.F90:1039
prg_response_mod::prg_canon_response
subroutine, public prg_canon_response(P1_bml, H1_bml, Nocc, beta, evals, mu0, m, HDIM)
First-order Canonical Density Matrix Perturbation Theory.
Definition: prg_response_mod.F90:948
prg_response_mod::prg_project_response
subroutine, public prg_project_response(rsp_bml, over_bml, spindex, norbi, coordinates, rspfunc, verbose)
Project the response onto atomic positions. First order response to the perturbation ( ) projected on...
Definition: prg_response_mod.F90:798
prg_response_mod::pi
real(dp), parameter pi
Definition: prg_response_mod.F90:19
prg_response_mod::prg_compute_response_rs
subroutine, public prg_compute_response_rs(ham_bml, prt_bml, rsp_bml, lambda, bndfil, threshold, verbose)
Computes the first order response density matrix using Rayleigh Schrodinger Perturbation theory The t...
Definition: prg_response_mod.F90:252
prg_response_mod::prg_compute_dipole
subroutine, public prg_compute_dipole(charges, coordinate, dipoleMoment, factor, verbose)
To compute the dipole moment of the system. The units of the dipole moment are determined by the unit...
Definition: prg_response_mod.F90:123
prg_response_mod::prg_canon_response_vector
subroutine, public prg_canon_response_vector(P1_bml, H1_bml, Nocc, beta, evals, mu0, m, thresh, HDIM)
First-order Canonical Density Matrix Perturbation Theory.
Definition: prg_response_mod.F90:850
prg_response_mod::prg_parse_response
subroutine, public prg_parse_response(RespData, filename)
The parser for the calculation of the DM response.
Definition: prg_response_mod.F90:46
prg_response_mod::prg_canon_response_p1_dpdmu
subroutine, public prg_canon_response_p1_dpdmu(P1_bml, dPdMu, H1_bml, Norbs, beta, Q_bml, evals, mu0, m, HDIM)
First-order Canonical Density Matrix Perturbation Theory.
Definition: prg_response_mod.F90:1113
prg_response_mod::prg_write_dipole_tcl
subroutine, public prg_write_dipole_tcl(dipoleMoment, file, factor, verbose)
To visualize a dipole moment using VMD. This will prg_generate a .tcl script that could be run using ...
Definition: prg_response_mod.F90:161
prg_response_mod::prg_compute_polarizability
subroutine, public prg_compute_polarizability(rsp_bml, prt_bml, polarizability, factor, verbose)
To compute the polarizability of the system. The units of the directional polarizability are determin...
Definition: prg_response_mod.F90:201
prg_response_mod::prg_compute_response_sp2
subroutine, public prg_compute_response_sp2(ham_bml, prt_bml, rsp_bml, rho_bml, lambda, bndfil, minsp2iter, maxsp2iter, sp2conv, idemtol, threshold, verbose)
Finds the first order response matrix from a Hamiltonian matrix.
Definition: prg_response_mod.F90:655
prg_response_mod::prg_pert_sin_pot
subroutine, public prg_pert_sin_pot(direction, lx, coordinate, lambda, prt_bml, threshold, spindex, norbi, verbose, over_bml)
Apply a sinusoidal length dependent potential ( ) where is the x coordinate. The Hamiltonian gets mo...
Definition: prg_response_mod.F90:526
prg_response_mod::prg_compute_response_fd
subroutine, public prg_compute_response_fd(ham_bml, prt_bml, rsp_bml, prg_delta, bndfil, threshold, verbose)
Computes the first order response density matrix using finite differences. The transformation hereby ...
Definition: prg_response_mod.F90:382
prg_response_mod::prg_pert_cos_pot
subroutine, public prg_pert_cos_pot(direction, lx, coordinate, lambda, prt_bml, threshold, spindex, norbi, verbose, over_bml)
Apply a cosine length dependent potential ( ) where is the x coordinate. The Hamiltonian gets modifi...
Definition: prg_response_mod.F90:594
prg_response_mod::prg_pert_from_file
subroutine, public prg_pert_from_file(prt_bml, norb)
Read perturbation from file.
Definition: prg_response_mod.F90:226
prg_response_mod::prg_pert_constant_field
subroutine, public prg_pert_constant_field(field, intensity, coordinate, lambda, prt_bml, threshold, spindex, norbi, verbose, over_bml)
Apply a constant field perturbation through the dipole moment operator ( ). In the matrix representat...
Definition: prg_response_mod.F90:449
prg_response_mod::respdata_type
Definition: prg_response_mod.F90:21