@brief Compute derivative coupling vectors (DCV) using the finite difference method, typically denoted as d_IJ between states I and J.
d_IJ = <I| d/dR |J> = h_IJ / (E_I - E_J),
where h_IJ is the nonadiabatic coupling, defined as
h_IJ = <I| dH/dR |J>.
Note that R is an entire vector, and so is d_IJ.
This routine computes d_IJ^a = <I| d/da |J>, which is a single component of d_IJ by TLF.
The component a can be either a geometric or a time derivative.
The geometric derivative is used to construct the DCV, while
the time derivative is used as NACME for nonadiabatic MD.
The d_IJ^a is computed by numerical differentiation using
a first or second-order formula. In the case of first-order,
O_IJ = <I(a)|J(a+da)> = F,
O_IJ = <I(a+da)|J(a)> = B,
d_IJ^a = (F - B) / a,
where O_IJ_F and O_IJ_B are the state overlaps.
This routine returns d_IJ^a * a = F - B.
The denominator MUST be provided in subsequent calculations
because it can be 2 * a in the case of second-order
numerical differentiation formula. Also, a can be either a time
or geometric parameter.
https://pubs.acs.org/doi/full/10.1021/acs.jctc.8b01049,
@author Seunghoon Lee, Konstantin Komarov
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=dp), | dimension(nstates,*) | :: | nact | |||
real(kind=dp), | dimension(nstates,*) | :: | s_st | |||
integer | :: | nstates |
subroutine get_dcv(nact, s_st, nstates) use precision, only: dp use io_constants, only: iw implicit none integer :: nstates real(kind=dp), dimension(nstates,*) :: nact, s_st integer :: i, j logical :: debug = .true. if (debug) write(iw, & fmt='(/x/,a,/29x," F B F - B")') & ' F = <I(a)|J(a+da)>, B = <I(a+da)|J(a)>, where a is variable.' do i = 1, nstates do j = 1, nstates nact(i,j) = (s_st(i,j)-s_st(j,i)) if (debug) write(iw, & fmt='(x,a,i0,a,i0,a,i0,a,i0,a,f12.8,a,f12.8,f12.8)') & "<S",i,"|S",j,"> and <S",j,"|S",i,"> = ",s_st(I,J), & " and", s_st(J,I), nact(i,j) end do end do end subroutine get_dcv