diff --git a/LR_Modules/adddvhubscf.f90 b/LR_Modules/adddvhubscf.f90 index b05b79caf..29cc3008c 100644 --- a/LR_Modules/adddvhubscf.f90 +++ b/LR_Modules/adddvhubscf.f90 @@ -12,14 +12,18 @@ SUBROUTINE adddvhubscf (ipert, ik) ! !! DFPT+U !! This routine calculates the SCF derivative of the Hubbard potential times \(\psi\): - !! $$ |\Delta V_{SCF}(k+q,is) \psi(\text{ibnd},k,is)\rangle = - !! - \sum_{I,m1,m2} \text{Hubbard}_U(I)\cdot \text{dnsscf}(m1,m2,is,I,\text{imode})\cdot - !! |S\phi(I,k+q,m1)\rangle\langle S\phi(I,k,m2|\psi(\text{ibnd},k,is)\rangle $$ + !! \begin{equation}\notag + !! \begin{split} + !! |\Delta V_{SCF}(k+q,is) \psi(\text{ibnd},k,is)\rangle = + !! - &\sum_{I,m1,m2} \text{Hubbard_U}(I)\cdot\text{dnsscf}(m1,m2,is,I,\text{imode})\cdot \\ + !! &|S\phi(I,k+q,m1)\rangle\langle S\phi(I,k,m2|\psi(\text{ibnd},k,is)\rangle + !! \end{split} + !! \end{equation} ! - !! Addition of the \(text{J0}\) terms: + !! Addition of the \(\text{J0}\) terms: ! !! $$ + \sum_{I,m1,m2} \text{Hubbard_J0}(I)\cdot \text{dnsscf}(m1,m2,\text{isi},I,\text{imode})\cdot - !! \|S\phi(I,k+q,m1)\rangle\langle S\phi(I,k,m2)\|\psi(\text{ibnd},k,is)\rangle $$ + !! |S\phi(I,k+q,m1)\rangle\langle S\phi(I,k,m2)|\psi(\text{ibnd},k,is)\rangle $$ ! !! Where: !! \(\text{is}\) = current_spin; diff --git a/PHonon/PH/alpha2f.f90 b/PHonon/PH/alpha2f.f90 index 32edd2bce..28f7af6ad 100644 --- a/PHonon/PH/alpha2f.f90 +++ b/PHonon/PH/alpha2f.f90 @@ -140,7 +140,7 @@ END SUBROUTINE read_polarization !-------------------------------------------------------------------- SUBROUTINE read_lam() !------------------------------------------------------------------ - !! This routine reads \(\text{lambda}_{q nu}\) & \(\text{omega}_{q nu}\) + !! This routine reads \(\text{lambda}_\text{q nu}\) & \(\text{omega}_\text{q nu}\) !! from lambda*.dat ! USE kinds, ONLY : DP @@ -296,7 +296,7 @@ END SUBROUTINE compute_a2F !----------------------------------------------------------------- SUBROUTINE compute_lambda() !--------------------------------------------------------------- - !! This routine computes \(\text{omega_ln}\) & \(\text{lambda}\). + !! This routine computes \(\text{omega_{ln}}\) & \(\text{lambda}\). ! USE kinds, ONLY : DP USE modes, ONLY : nmodes diff --git a/PHonon/PH/delta_sphi.f90 b/PHonon/PH/delta_sphi.f90 index fb77f02ee..78638ed25 100644 --- a/PHonon/PH/delta_sphi.f90 +++ b/PHonon/PH/delta_sphi.f90 @@ -12,7 +12,7 @@ SUBROUTINE delta_sphi (ikk, ikq, na, icart, nah, ihubst, wfcatomk_, wfcatomkpq_, dqsphi, dmqsphi, iflag) !--------------------------------------------------------------------------------- ! - ! DFPT+U: This routine calculates a vector at k: + !! DFPT+U: This routine calculates a vector at k ! ! |\Delta_{-q}(S_{k+q} \phi_(k+q,I,m)) > = S_{k} | \d^{icart}phi_(k,nah,m) > + ! \sum{l1,l2} [ | \dbeta^{icar_}_(k,na_,l1) > * qq_nt(na_, l1 ,l2) * @@ -20,7 +20,7 @@ SUBROUTINE delta_sphi (ikk, ikq, na, icart, nah, ihubst, wfcatomk_, wfcatomkpq_, ! | \beta_(k,na_,l1)> * qq_nt(na_, l1 ,l2) * ! < \dbeta^{icar_}_(k+q,na_,l2) | phi_(k+q,nah,m) > ] ! - !! and also a vector at k+q : + !! and also a vector at k+q. ! ! |\Delta_q(S_{k} \phi_(k,I,m)) > = S_{k+q}| \d^{icart}phi_(k+q,nah,m) > + ! \sum{l1,l2} [ | \dbeta^{icar_}_(k+q,na_,l1) > * qq_nt(na_, l1 ,l2) * @@ -32,8 +32,10 @@ SUBROUTINE delta_sphi (ikk, ikq, na, icart, nah, ihubst, wfcatomk_, wfcatomkpq_, ! |\Delta_{-q}(S_{k+q} \phi_(k+q,I,m)) > ! iflag = 0 : calculate ONLY |\Delta_{-q}(S_{k+q} \phi_(k+q,I,m)) > ! - ! Written by A. Floris - ! Modified by I. Timrov (01.10.2018) + !! See source comment for details on the implemented formulas. + ! + !! Written by A. Floris. + !! Modified by I. Timrov (01.10.2018). ! USE kinds, ONLY : DP USE uspp_param, ONLY : nh, nhm diff --git a/PHonon/PH/dvqpsi_us_only.f90 b/PHonon/PH/dvqpsi_us_only.f90 index 798809513..13328ea3a 100644 --- a/PHonon/PH/dvqpsi_us_only.f90 +++ b/PHonon/PH/dvqpsi_us_only.f90 @@ -9,8 +9,9 @@ !---------------------------------------------------------------------- subroutine dvqpsi_us_only (ik, uact, becp1, alphap) !---------------------------------------------------------------------- - !! This routine calculates dV_bare/dtau * psi for one perturbation - !! with a given q. The displacements are described by a vector uact. + !! This routine calculates \(\text{dV_bare}/\text{dtau}\cdot\text{psi}\) + !! for one perturbation with a given q. + !! The displacements are described by a vector uact. !! The result is stored in \(\text{dvpsi}\). The routine is called for !! each k-point and for each pattern u. It computes simultaneously all !! the bands. diff --git a/PHonon/PH/dvscf_interpolate.f90 b/PHonon/PH/dvscf_interpolate.f90 index a6eb9fdb0..e75ad0bfc 100644 --- a/PHonon/PH/dvscf_interpolate.f90 +++ b/PHonon/PH/dvscf_interpolate.f90 @@ -10,8 +10,8 @@ MODULE dvscf_interpolate !---------------------------------------------------------------------------- !! !! Module for Fourier interpolation of phonon potential dvscf. - !! Uses the output of dvscf_q2r.x program to compute the induced part of - !! dvscf at each q point. + !! Uses the output of \(\texttt{dvscf_q2r.x}\) program to compute the induced + !! part of \(\text{dvscf}\) at each q point. ! ! See header of dvscf_q2r.f90 for details. ! @@ -61,8 +61,8 @@ MODULE dvscf_interpolate !---------------------------------------------------------------------------- SUBROUTINE dvscf_interpol_setup() !-------------------------------------------------------------------------- - !! Do setups for dvscf_r2q. Called in phq_setup at the beginning of each - !! q point calculation. + !! Do setups for \(\texttt{dvscf_r2q}\). Called in \(\texttt{phq_setup}\) + !! at the beginning of each q-point calculation. !-------------------------------------------------------------------------- ! USE kinds, ONLY : DP @@ -234,20 +234,20 @@ MODULE dvscf_interpolate !---------------------------------------------------------------------------- SUBROUTINE dvscf_r2q(xq, u_in, dvscf) !-------------------------------------------------------------------------- - !! !! Read inverse Fourier transformed potential w_pot (written by - !! dvscf_q2r.x) and Fourier transform to compute dvscf at given q point. - !! - !! Originally proposed by [1], long-range part described in [2]. - !! [1] Eiguren and Ambrosch-Draxl, PRB 78, 045124 (2008) + !! \(\texttt{dvscf_q2r.x}\)) and Fourier transform to compute + !! \(\text{dvscf}\) at given q point. + ! + !! Originally proposed by [1], long-range part described in [2]. + !! [1] Eiguren and Ambrosch-Draxl, PRB 78, 045124 (2008) !! [2] Xavier Gonze et al, Comput. Phys. Commun., 107042 (2019) - !! - !! dvscf(r,q) = exp(-iqr) (dvlong(r,q) + sum_R exp(iqR) w_pot(r,R)) - !! + ! + !! $$ \text{dvscf}(r,q) = \text{exp}(-iqr) (\text{dvlong}(r,q) + + !! \sum_R \text{exp}(iqR) \text{w_pot}(r,R)) $$ + ! !! In this subroutine, pool parallelization is used to distribute R points !! to nodes so that the root of each pool reads different w_pot file !! simultaneously, reducing the io time. - !! !-------------------------------------------------------------------------- ! USE kinds, ONLY : DP @@ -473,15 +473,14 @@ MODULE dvscf_interpolate !---------------------------------------------------------------------------- SUBROUTINE dvscf_shift_center(dvscf_q, xq, shift_half, sign) !---------------------------------------------------------------------------- - !! Shift center of phonon potential. - !! If sign = +1, shift center from origin to tau (used in dvscf_r2q). + !! Shift center of phonon potential. + !! If sign = +1, shift center from origin to tau (used in dvscf_r2q). !! If sign = -1, shift center from tau to origin (used in dvscf_q2r). - !! + ! !! For ipol = 1, 2, 3, if shift_half(ipol) is true, the origin is set at !! r = 0.5 for direction ipol. If false, the origin is set at r = 0.0. !! Setting shift_half = .true. is useful when the size of supercell is odd, !! becuase the center of the supercell is at (0.5, 0.5, 0.5). - !! !---------------------------------------------------------------------------- USE kinds, ONLY : DP USE constants, ONLY : tpi @@ -573,19 +572,20 @@ MODULE dvscf_interpolate !---------------------------------------------------------------------------- !! This subroutine calculates the long-range dipole potential for given !! xq and zeu. - !! Input xq: q vector in Cartesian coordinate - !! + !! Input xq: q vector in Cartesian coordinate. + ! !! Currently, only the dipole part (Frohlich) is implemented. The quadrupole !! potential is not implemented. - !! + ! !! [1] Xavier Gonze et al, Comput. Phys. Commun., 107042 (2019) - !! + ! !! Taken from Eq.(13) of Ref. [1] - !! dvlong(G,q)_{a,x} = 1j * 4pi / Omega * e^2 - !! * [ (q+G)_y * Zstar_{a,yx} * exp(-i*(q+G)*tau_a)) ] - !! / [ (q+G)_y * epsilon_yz * (q+G)_z ] + !! $$ \text{dvlong}(G,q)_{a,x} = 1j \cdot 4\pi / \text{omega} \cdot e^2 + !! \cdot [ (q+G)_y \cdot \text{Zstar}_{a,yx} \cdot \text{exp} + !! (-i\cdot(q+G)\cdot\text{tau}_a)) ] + !! / [ (q+G)_y \cdot \epsilon_yz \cdot (q+G)_z ] $$ !! a: atom index, x, y: Cartesian direction index - !! + ! !! Since QE uses Rydberg units, we multiply the e^2 = 2.0 factor. !! The units of q and G are 2pi/a, so we need to divide dvlong by tpiba. !! diff --git a/PHonon/PH/dwfc.f90 b/PHonon/PH/dwfc.f90 index 3863c0544..8cc49c944 100644 --- a/PHonon/PH/dwfc.f90 +++ b/PHonon/PH/dwfc.f90 @@ -27,7 +27,7 @@ SUBROUTINE dwfc (npw_ , igk_ , ik_ , icart_ , func, dfunc) COMPLEX(DP), INTENT(IN) :: func(npwx) COMPLEX(DP), INTENT(OUT) :: dfunc(npwx) ! - ! Local variables + ! ... local variables ! INTEGER :: ig REAL(DP) :: gvec, xk_aux @@ -55,7 +55,7 @@ SUBROUTINE d2wfc (npw_ , igk_ , ik_ , icart_ , jcart_, func, d2func) !! This routine calculates the second derivative of a wave function !! with respect to the position operator \(r_\text{icart}\) and !! \(r_\text{jcart}\): - !! $$ d^2 \psi/d\text{icart}d\text{jcart} = \sum_G [-(k+G)_\text{icart} + !! $$ d^2 \psi/d\text{icart}\ d\text{jcart} = \sum_G [-(k+G)_\text{icart} !! \cdot (-(k+G)_\text{gcart})] \psi(G) $$ ! USE kinds, ONLY : DP @@ -72,7 +72,7 @@ SUBROUTINE d2wfc (npw_ , igk_ , ik_ , icart_ , jcart_, func, d2func) COMPLEX(DP), INTENT(IN) :: func(npwx) COMPLEX(DP), INTENT(OUT) :: d2func(npwx) ! - ! Local variables + ! ... local variables ! INTEGER :: ig REAL(DP) :: gvec, xk_aux, gvec2, xk_aux2 diff --git a/PHonon/PH/dynmat_hub_bare.f90 b/PHonon/PH/dynmat_hub_bare.f90 index 70e2992fb..71c317a7d 100644 --- a/PHonon/PH/dynmat_hub_bare.f90 +++ b/PHonon/PH/dynmat_hub_bare.f90 @@ -9,19 +9,23 @@ SUBROUTINE dynmat_hub_bare !--------------------------------------------------------------------- !! DFPT+U: This routine does two tasks: - !! 1) Computes d2ns_bare (very slow) or reads it; + !! 1) Computes \(\text{d2ns_bare}\) (very slow) or reads it; !! 2) Adds to the dynamical matrix the bare part due the Hubbard term: - !! $$ \sum_\text{nah} \text{Hubbard_U}(\text{nah}) \sum_{is, m1, m2} [ - !! (0.5\delta_{m1m2} - \text{ns}(m1,m2,is,\text{nah})) \cdot \text{d2ns_bare} - !! - \text{conjg}(\text{dnsbare}( m1, m2, is, \text{nah}, \text{icart},na)) - !! \cdot \text{dnsbare}( m1, m2, is, nah, jcart,nap) ] $$ + !! \begin{equation}\notag + !! \begin{split} + !! \sum_\text{nah} \text{Hubbard_U}&(\text{nah}) \sum_{is, m1, m2} [ + !! (0.5\delta_{m1m2} - \text{ns}(m1,m2,is,\text{nah})) \cdot \text{d2ns_bare}\\ + !! &- \text{conjg}(\text{dnsbare}( m1, m2, is, \text{nah}, \text{icart},\text{na})) + !! \cdot \text{dnsbare}( m1, m2, is, \text{nah}, \text{jcart},\text{nap}) ] + !! \end{split} + !! \end{equation} ! !! Addition of the J0 term: !! $$ \sum_\text{nah} \text{Hubbard_J0}(\text{nah}) \sum_{is, m1, m2} [ !! \text{ns}(m1,m2,-is,\text{nah}) \cdot \text{d2ns_bare}(m1, m2, is, !! \text{nah}, \text{nap_jcart}, \text{na_icart}) - !! + \text{conjg}(\text{dnsbare}( m1, m2, is, \text{nah}, \text{icart},na))\cdot - !! \text{dnsbare}( m1, m2, -is, nah, jcart,nap) ] $$ + !! + \text{conjg}(\text{dnsbare}(m1,m2,is,\text{nah}, \text{icart},\text{na}))\cdot + !! \text{dnsbare}(m1,m2,-is, \text{nah}, \text{jcart},\text{nap}) ] $$ ! !! Written by A. Floris. !! Modified by I. Timrov (01.10.2018). diff --git a/PHonon/PH/dynmat_hub_scf.f90 b/PHonon/PH/dynmat_hub_scf.f90 index 9ce2f28c8..91c1fd9eb 100644 --- a/PHonon/PH/dynmat_hub_scf.f90 +++ b/PHonon/PH/dynmat_hub_scf.f90 @@ -13,9 +13,10 @@ SUBROUTINE dynmat_hub_scf (irr, nu_i0, nper) !! It adds the scf derivative of the Hubbard energy (terms 1 and 2). Moreover, !! it adds 3 terms due to the orthogonality constraints in the USPP formalism (terms 3,4,5). !! Note, the orthogonality terms 4 and 5 DO NOT follow simply from the 2nd derivative of - !! the Hubbard energy \(\text{E_U}\); they stem from the coupling of the \(-\epsilon\dS\d\lambda\) in - !! the USPP forces with the part of d\psi\d\mu projected onto the occupied manifold - !! (the variable dpsi stems from the scf linear system and lives in the unoccupied manifold). + !! the Hubbard energy \(\text{E_U}\); they stem from the coupling of the \(-\epsilon\ dS\ d\lambda\) in + !! the USPP forces with the part of \(d\psi\ d\mu\) projected onto the occupied manifold + !! (the variable dpsi stems from the scf linear system and lives in the unoccupied manifold). + !! See comments in the source code for other details on the implemented terms. ! ! Terms implemented: ! dyn_hub_scf (ipert, imode) diff --git a/PHonon/PH/gmressolve_all.f90 b/PHonon/PH/gmressolve_all.f90 index 0c7eafcfe..bdfe2d6b2 100644 --- a/PHonon/PH/gmressolve_all.f90 +++ b/PHonon/PH/gmressolve_all.f90 @@ -11,7 +11,7 @@ subroutine gmressolve_all (h_psi, cg_psi, e, d0psi, dpsi, h_diag, & ndmx, ndim, ethr, ik, kter, conv_root, anorm, nbnd, m) !---------------------------------------------------------------------- !! Iterative solution of the linear system by GMRES(m) method: - !! $$ (h - e + Q) \cdot \(\text{dpsi} = \(\text{d0psi}\)\ \ (1) $$ + !! $$ (h - e + Q) \cdot \text{dpsi} = \text{d0psi}\ \ (1) $$ !! where \(h\) is a complex hermitian matrix, \(e\) is a !! complex scalar and \(\text{dpsi}\) and \(\text{d0psi}\) are !! complex vectors. diff --git a/PHonon/PH/init_representations.f90 b/PHonon/PH/init_representations.f90 index daab7f861..19c53eda2 100644 --- a/PHonon/PH/init_representations.f90 +++ b/PHonon/PH/init_representations.f90 @@ -113,8 +113,8 @@ subroutine initialize_grid_variables() !---------------------------------------------------------------------- !! This subroutine initializes the grid variables by reading the !! modes from file. It uses the routine check_if_partial_dyn to - !! set the modes to compute according to \(\text{start_irr}, - !! \(\text{last_irr} or \(\text{modenum}\) and \(\text{ifat}\) flags. + !! set the modes to compute according to \(\text{start_irr}\), + !! \(\text{last_irr}\) or \(\text{modenum}\) and \(\text{ifat}\) flags. ! USE kinds, ONLY : DP USE ions_base, ONLY : nat diff --git a/PHonon/PH/phq_setup.f90 b/PHonon/PH/phq_setup.f90 index 1fa936239..d3d4206f8 100644 --- a/PHonon/PH/phq_setup.f90 +++ b/PHonon/PH/phq_setup.f90 @@ -11,26 +11,26 @@ subroutine phq_setup !! This subroutine prepares several variables which are needed in the !! \(\texttt{phonon}\) program: !! 1) computes the total local potential (external+scf) on the smooth - !! grid to be used in h_psi and similia; + !! grid to be used in \(\texttt{h_psi}\) and similia; !! 2) computes the local magnetization (if necessary); !! 3) computes dmuxc (with GC if needed); !! 4) set the inverse of every matrix invs; !! 5) for metals sets the occupied bands; - !! 6) computes alpha_pv; + !! 6) computes \(\text{alpha_pv}\); !! 7) computes the variables needed to pass to the pattern representation: !! \(\text{u}\): the patterns; !! \(\text{t}\): the matrices of the small group of q on the pattern basis; - !! \(\text{tmq}\): the matrix of the symmetry which sends \(q\rightarrow-q+G}\); + !! \(\text{tmq}\): the matrix of the symmetry which sends \(q\rightarrow -q+G\); !! \(\text{gi}\): the G associated to each symmetry operation; - !! \(\text{gimq}\): the G of the \(q\rigtharrow-q+G}\) symmetry; - !! \(\text{nsymq}\): the order of the small group of q; - !! \(\text{irotmq}\): the index of the \(q\rightarrow-q+G}\) symmetry; + !! \(\text{gimq}\): the G of the \(q\rightarrow -q+G\) symmetry; + !! \(\text{nsymq}\): the order of the small group of \(q\); + !! \(\text{irotmq}\): the index of the \(q\rightarrow -q+G\) symmetry; !! \(\text{nirr}\): the number of irreducible representation; !! \(\text{npert}\): the dimension of each irreducible representation; !! \(\text{nmodes}\): the number of modes; - !! \(\text{minus_q}\): true if there is a symmetry sending \(q\rightarrow-q+G}\); - !! 8) for testing purposes it sets ubar; - !! 9) set the variables needed to deal with nlcc; + !! \(\text{minus_q}\): true if there is a symmetry sending \(q\rightarrow -q+G\); + !! 8) for testing purposes it sets \(\text{ubar}\); + !! 9) set the variables needed to deal with \(\text{nlcc}\); !! 10) set the variables needed for the partial computation ! of the dynamical matrix. ! diff --git a/PHonon/PH/q2qstar_ph.f90 b/PHonon/PH/q2qstar_ph.f90 index 22a7614a8..2ae8b55b6 100644 --- a/PHonon/PH/q2qstar_ph.f90 +++ b/PHonon/PH/q2qstar_ph.f90 @@ -11,7 +11,7 @@ subroutine q2qstar_ph (dyn, at, bg, nat, nsym, s, invs, irt, rtau, & !----------------------------------------------------------------------- !! Generates the dynamical matrices for the star of q and writes them on !! disk for later use. - !! If there is a symmetry operation such that \(q \rigtharrow -q+G \) then + !! If there is a symmetry operation such that \(q \rightarrow -q+G \) then !! imposes on dynamical matrix those conditions related to time reversal !! symmetry. !