diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 94b6640f90..e9bb2959dd 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -375,7 +375,7 @@ - [out\_wannier\_eig](#out_wannier_eig) - [out\_wannier\_unk](#out_wannier_unk) - [out\_wannier\_wvfn\_formatted](#out_wannier_wvfn_formatted) - - [rt-TDDFT: Real-time time dependent density functional theory](#tddft-time-dependent-density-functional-theory) + - [RT-TDDFT: Real-Time Time-Dependent Density Functional Theory](#rt-tddft-real-time-time-dependent-density-functional-theory) - [estep\_per\_md](#estep_per_md) - [td\_dt](#td_dt) - [td\_edm](#td_edm) @@ -3694,44 +3694,45 @@ These variables are used to control berry phase and wannier90 interface paramete [back to top](#full-list-of-input-keywords) -## TDDFT: time dependent density functional theory +## RT-TDDFT: Real-Time Time-Dependent Density Functional Theory ### estep_per_md - **Type**: Integer -- **Description**: The number of electron propagation steps between two ionic steps. +- **Description**: The number of electronic propagation steps between two ionic steps. - **Default**: 1 ### td_dt - **Type**: Real -- **Description**: The time step used in electron propagation. Setting td_dt will reset the md_dt value to td_dt * estep_per_md. -- **Default**: md_dt/estep_per_md +- **Description**: The time step used in electronic propagation. Setting `td_dt` will reset the value of [`md_dt`](#md_dt) to `td_dt * estep_per_md`. +- **Default**: `md_dt / estep_per_md` +- **Unit**: fs ### td_edm - **Type**: Integer -- **Description**: Method to calculate the energy density matrix - - 0: new method (use the original formula). - - 1: old method (use the formula for ground state). +- **Description**: Method to calculate the energy-density matrix, mainly affects the calculation of force and stress. + - 0: Using the original formula: $\mathrm{EDM}_{\mu\nu}=\frac{1}{2} \sum_{\eta \zeta}\left(S_{\mu \eta}^{-1} H_{\eta \zeta} \rho_{\zeta \nu}+\rho_{\mu \eta} H_{\eta \zeta} S_{\zeta \nu}^{-1}\right)$. + - 1: Using the formula for ground state (deprecated): $\mathrm{EDM}_{\mu\nu}=\sum_n^{\mathrm{occ}} c_{\mu,n} f_n \epsilon_n c_{n,\nu}$. Note that this usually does not hold if wave function is not the eigenstate of the Hamiltonian. - **Default**: 0 ### td_print_eij - **Type**: Real -- **Description**: - - <0: don't print $E_{ij}$. - - \>=0: print the $E_{ij}\ (<\psi_i|H|\psi_j>$) elements which are larger than td_print_eij. +- **Description**: Controls the printing of Hamiltonian matrix elements $E_{ij}=\Braket{\psi_i|\hat{H}|\psi_j}$. + - $<0$: Suppress all $E_{ij}$ output. + - $\geqslant 0$: Print only elements with ​​either $|\text{Re}(E_{ij})|$ or $|\text{Im}(E_{ij})|$​​ exceeding `td_print_eij`. - **Default**: -1 +- **Unit**: Ry ### td_propagator - **Type**: Integer -- **Description**: - Methods of propagator - - 0: Crank-Nicolson, based on matrix inversion. - - 1: 4th Taylor expansions of exponential. - - 2: enforced time-reversal symmetry (ETRS). +- **Description**: Methods of electronic propagation: $\psi_{n\boldsymbol{k}}(\boldsymbol{r},t_2) = U(t_2,t_1) \psi_{n\boldsymbol{k}}(\boldsymbol{r},t_1)$. + - 0: Crank-Nicolson, based on matrix inversion: $U = \dfrac{S_{\boldsymbol{k}}-\mathrm{i} H_{\boldsymbol{k}}((t_1+t_2)/2) \Delta t / 2}{S_{\boldsymbol{k}}+\mathrm{i} H_{\boldsymbol{k}}((t_1+t_2)/2) \Delta t / 2}$. + - 1: 4th-order Taylor expansion of exponential: $U = I + \hat{A} + \frac{1}{2}\hat{A}^2 + \frac{1}{6}\hat{A}^3 + \frac{1}{24}\hat{A}^4$, where $\hat{A} = -\mathrm{i} S_{\boldsymbol{k}}^{-1} H_{\boldsymbol{k}}((t_1+t_2)/2)\Delta t$. + - 2: Enforced time-reversal symmetry (ETRS): $U = \mathrm{e}^{-\mathrm{i} S_{\boldsymbol{k}}^{-1} H_{\boldsymbol{k}}(t_2)\frac{\Delta t}{2}} \mathrm{e}^{-\mathrm{i} S_{\boldsymbol{k}}^{-1} H_{\boldsymbol{k}}(t_1)\frac{\Delta t}{2}}$. - 3: Crank-Nicolson, based on solving linear equation. - **Default**: 0 @@ -3739,53 +3740,75 @@ These variables are used to control berry phase and wannier90 interface paramete - **Type**: Boolean - **Description**: - - True: add a laser material interaction (extern laser field). - - False: no extern laser field. + - True: Add a laser-material interaction (external electric field). + - False: No external electric field. - **Default**: False ### td_vext_dire - **Type**: String - **Description**: - If `td_vext` is True, the td_vext_dire is a string to set the number of electric fields, like `td_vext_dire 1 2` representing external electric field is added to the x and y axis at the same time. Parameters of electric field can also be written as a string, like `td_gauss_phase 0 1.5707963267948966` representing the Gauss field in the x and y directions has a phase delay of Pi/2. See below for more parameters of electric field. - - 1: the direction of external light field is along x axis. - - 2: the direction of external light field is along y axis. - - 3: the direction of external light field is along z axis. + Specifies the direction(s) of the external electric field when `td_vext` is enabled. For example, `td_vext_dire 1 2` indicates that external electric fields are applied to both the x and y directions simultaneously. Electric field parameters can also be written as strings. For example, `td_gauss_phase 0 1.5707963` indicates that the Gaussian type electric fields in the x and y directions have a phase delay of $\pi/2$. See below for more electric field parameters. + - 1: The external field direction is along the x-axis. + - 2: The external field direction is along the y-axis. + - 3: The external field direction is along the z-axis. - **Default**: 1 + > Note: The axes refer to the absolute Cartesian coordinate system, which is independent of lattice vectors or reciprocal space definitions.​ This is different from [`efield_dir`](#efield_dir), which is used in ground-state KSDFT. + ### td_stype - **Type**: Integer -- **Description**: - Type of electric field in space domain - - 0: length gauge. - - 1: velocity gauge. - - 2: hybrid gauge. +- **Description**: Type of electric field in the space domain, i.e. the gauge of the electric field. + - 0: Length gauge. + - 1: Velocity gauge. + - 2: Hybrid gauge. See [*J. Chem. Theory Comput.* 2025, 21, 3335−3341](https://pubs.acs.org/doi/10.1021/acs.jctc.5c00111) for more information. - **Default**: 0 ### td_ttype -- **Type**: Integer +- **Type**: String - **Description**: - Type of electric field in time domain - - 0: Gaussian type function. - - 1: Trapezoid function. - - 2: Trigonometric function. - - 3: Heaviside function. - - 4: HHG function. + Type of electric field in the time domain. + - 0: Gaussian type function: + $$ + E(t) = A \cos\left[2\pi f(t-t_0)+\varphi\right]\exp\left[-\frac{(t-t_0)^2}{2\sigma^2}\right] + $$ + - 1: Trapezoid function: + $$ + E(t) = + \begin{cases} + A \cos(2\pi f t + \varphi) \cdot \dfrac{t}{t_1}, & t < t_1 \\ + A \cos(2\pi f t + \varphi), & t_1 \leqslant t \leqslant t_2 \\ + A \cos(2\pi f t + \varphi) \left(1 - \dfrac{t - t_2}{t_3 - t_2}\right), & t_2 < t < t_3 \\ + 0, & t \geqslant t_3 + \end{cases} + $$ + - 2: Trigonometric function: + $$ + E(t) = A \cos(2\pi f_1 t + \varphi_1) \sin^2(2\pi f_2 t + \varphi_2) + $$ + - 3: Heaviside step function: + $$ + E(t) = + \begin{cases} + A, & t < t_0 \\ + 0, & t \geqslant t_0 + \end{cases} + $$ - **Default**: 0 ### td_tstart - **Type**: Integer -- **Description**: Number of steps where electric field starts +- **Description**: The initial time step when the time-dependent electric field is activated. - **Default**: 1 ### td_tend - **Type**: Integer -- **Description**: Number of steps where electric field ends -- **Default**: 100 +- **Description**: The final time step when the time-dependent electric field is deactivated. The field remains active between `td_tstart` and `td_tend`. +- **Default**: 1000 ### td_lcut1 @@ -3793,9 +3816,9 @@ These variables are used to control berry phase and wannier90 interface paramete - **Description**: The lower bound of the interval in the length gauge RT-TDDFT, where $x$ is the fractional coordinate: - $$ + $ E(x)=\begin{cases}E_0, & \mathtt{cut1}\leqslant x \leqslant \mathtt{cut2} \\-E_0\left(\dfrac{1}{\mathtt{cut1}+1-\mathtt{cut2}}-1\right), & 0 < x < \mathtt{cut1~~or~~cut2} < x < 1 \end{cases} - $$ + $ - **Default**: 0.05 @@ -3805,230 +3828,201 @@ These variables are used to control berry phase and wannier90 interface paramete - **Description**: The upper bound of the interval in the length gauge RT-TDDFT, where $x$ is the fractional coordinate: - $$ + $ E(x)=\begin{cases}E_0, & \mathtt{cut1}\leqslant x \leqslant \mathtt{cut2} \\-E_0\left(\dfrac{1}{\mathtt{cut1}+1-\mathtt{cut2}}-1\right), & 0 < x < \mathtt{cut1~~or~~cut2} < x < 1 \end{cases} - $$ + $ - **Default**: 0.95 ### td_gauss_freq -- **Type**: Real +- **Type**: String - **Description**: - Frequency (freq) of Gauss type electric field (fs^-1)\ - amp\*cos(2pi\*freq(t-t0)+phase)exp(-(t-t0)^2/2sigma^2) + Frequency $f$ of the Gaussian type electric field. - **Default**: 22.13 +- **Unit**: 1/fs ### td_gauss_phase -- **Type**: Real +- **Type**: String - **Description**: - Phase of Gauss type electric field\ - amp\*(2pi\*freq(t-t0)+phase)exp(-(t-t0)^2/2sigma^2) + Phase $\varphi$ of the Gaussian type electric field. - **Default**: 0.0 ### td_gauss_sigma -- **Type**: Real +- **Type**: String - **Description**: - Sigma of Gauss type electric field (fs)\ - amp\*cos(2pi\*freq(t-t0)+phase)exp(-(t-t0)^2/2sigma^2) + Pulse width (standard deviation) $\sigma$ of the Gaussian type electric field. - **Default**: 30.0 +- **Unit**: fs ### td_gauss_t0 -- **Type**: Real +- **Type**: String - **Description**: - Step number of time center (t0) of Gauss type electric field\ - amp\*cos(2pi\*freq(t-t0)+phase)exp(-(t-t0)^2/2sigma^2) + Step number of the time center $t_0$ of the Gaussian type electric field. - **Default**: 100 ### td_gauss_amp -- **Type**: Real +- **Type**: String - **Description**: - Amplitude (amp) of Gauss type electric field (V/Angstrom)\ - amp\*cos(2pi\*freq(t-t0)+phase)exp(-(t-t0)^2/2sigma^2) + Amplitude $A$ of the Gaussian type electric field. - **Default**: 0.25 +- **Unit**: V/Å ### td_trape_freq -- **Type**: Real +- **Type**: String - **Description**: - Frequency (freq) of Trapezoid type electric field (fs^-1)\ - E = amp\*cos(2pi\*freq\*t+phase) t/t1 , tt3 + Frequency $f$ of the trapezoid type electric field. - **Default**: 1.60 +- **Unit**: 1/fs ### td_trape_phase -- **Type**: Real +- **Type**: String - **Description**: - Phase of Trapezoid type electric field\ - E = amp\*cos(2pi\*freq\*t+phase) t/t1 , tt3 + Phase $\varphi$ of the trapezoid type electric field. - **Default**: 0.0 ### td_trape_t1 -- **Type**: Real +- **Type**: String - **Description**: - Step number of time interval 1 (t1) of Trapezoid type electric field\ - E = amp\*cos(2pi\*freq\*t+phase) t/t1 , tt3 + Step number of the time interval $t_1$ of the trapezoid type electric field. - **Default**: 1875 ### td_trape_t2 -- **Type**: Real +- **Type**: String - **Description**: - Step number of time interval 2 (t2) of Trapezoid type electric field\ - E = amp\*cos(2pi\*freq\*t+phase) t/t1 , tt3 + Step number of the time interval $t_2$ of the trapezoid type electric field. - **Default**: 5625 ### td_trape_t3 -- **Type**: Real +- **Type**: String - **Description**: - Step number of time interval 3 (t3) of Trapezoid type electric field\ - E = amp\*cos(2pi\*freq\*t+phase) t/t1 , tt3 + Step number of the time interval $t_3$ of the trapezoid type electric field. - **Default**: 7500 ### td_trape_amp -- **Type**: Real +- **Type**: String - **Description**: - Amplitude (amp) of Trapezoid type electric field (V/Angstrom)\ - E = amp\*cos(2pi\*freq\*t+phase) t/t1 , tt3 + Amplitude $A$ of the trapezoid type electric field. - **Default**: 2.74 +- **Unit**: V/Å ### td_trigo_freq1 -- **Type**: Real +- **Type**: String - **Description**: - Frequency 1 (freq1) of Trigonometric type electric field (fs^-1)\ - amp\*cos(2\*pi\*freq1\*t+phase1)\*sin(2\*pi\*freq2\*t+phase2)^2 + Frequency $f_1$ of the trigonometric type electric field. - **Default**: 1.164656 +- **Unit**: 1/fs ### td_trigo_freq2 -- **Type**: Real +- **Type**: String - **Description**: - Frequency 2 (freq2) of Trigonometric type electric field (fs^-1)\ - amp\*cos(2\*pi\*freq1\*t+phase1)\*sin(2\*pi\*freq2\*t+phase2)^2 + Frequency $f_2$ of the trigonometric type electric field. - **Default**: 0.029116 +- **Unit**: 1/fs ### td_trigo_phase1 -- **Type**:Real +- **Type**:String - **Description**: - Phase 1 (phase1) of Trigonometric type electric field\ - amp\*cos(2\*pi\*freq1\*t+phase1)\*sin(2\*pi\*freq2\*t+phase2)^2 + Phase $\varphi_1$ of the trigonometric type electric field. - **Default**: 0.0 ### td_trigo_phase2 -- **Type**: Real +- **Type**: String - **Description**: - Phase 2 (phase2) of Trigonometric type electric field\ - amp\*cos(2\*pi\*freq1\*t+phase1)\*sin(2\*pi\*freq2\*t+phase2)^2 + Phase $\varphi_2$ of the trigonometric type electric field. - **Default**: 0.0 ### td_trigo_amp -- **Type**: Real +- **Type**: String - **Description**: - Amplitude (amp) of Trigonometric type electric field (V/Angstrom)\ - amp\*cos(2\*pi\*freq1\*t+phase1)\*sin(2\*pi\*freq2\*t+phase2)^2 + Amplitude $A$ of the trigonometric type electric field. - **Default**: 2.74 +- **Unit**: V/Å ### td_heavi_t0 -- **Type**: Real +- **Type**: String - **Description**: - Step number of switch time (t0) of Heaviside type electric field\ - E = amp , tt0 + Step number of the switch time $t_0$ of the Heaviside type electric field. - **Default**: 100 ### td_heavi_amp -- **Type**: Real +- **Type**: String - **Description**: - Amplitude (amp) of Heaviside type electric field (V/Angstrom)\ - E = amp , tt0 -- **Default**: 2.74 + Amplitude $A$ of the Heaviside type electric field. +- **Default**: 1.0 +- **Unit**: V/Å ### out_dipole - **Type**: Boolean - **Description**: - - True: output dipole. - - False: do not output dipole. + - True: Output electric dipole moment. + - False: Do not output electric dipole moment. - **Default**: False ### out_current - **Type**: Boolean -- **Description**: Output current in real time TDDFT simulations with the velocity gauge - - True: output current. - - False: do not output current. +- **Description**: + - True: Output current. + - False: Do not output current. - **Default**: False ### out_current_k - **Type**: Boolean -- **Description**: Output rt-TDDFT current for all k points - - True: output tddft current for all k points - - False: output current in total +- **Description**: + - True: Output current for each k-points separately. + - False: Output current in total. - **Default**: False ### out_efield - **Type**: Boolean -- **Description**: Output TDDFT Efield or not (V/Angstrom) - - True: output efield - - False: do not output efield +- **Description**: Whether to output the electric field data to files. When enabled, writes real-time electric field values (unit: ​​V/Å​​) into files named `efield_[num].txt`, where `[num]` is the ​​sequential index of the electric field ranges from `0` to `N-1` for `N` configured fields. It is noteworthy that the field type sequence follows [`td_ttype`](#td_ttype), while the direction sequence follows [`td_vext_dire`](#td_vext_dire). + - True: Output electric field. + - False: Do not output electric field. - **Default**: False ### out_vecpot - **Type**: Boolean -- **Description**: Output TDDFT Vector potential or not (a.u.) - - True: output Vector potential in file "OUT.suffix/At.dat" - - False: do not output Vector potential +- **Description**: Output vector potential or not (unit: a.u.). + - True: Output vector potential into file `At.dat`. + - False: Do not output vector potential. - **Default**: False ### init_vecpot_file - **Type**: Boolean -- **Description**: Init vector potential through file or not - - True: init vector potential from file "At.dat".(a.u.) It consists of four columns, representing istep and vector potential on each direction. - - False: calculate vector potential by integral of Efield +- **Description**: Initialize vector potential through file or not. + - True: Initialize vector potential from file `At.dat` (unit: a.u.). It consists of four columns, representing the step number and vector potential on each direction. + - False: Calculate vector potential by integrating the electric field. - **Default**: False ### ocp - **Type**: Boolean - **Availability**: - - For PW and LCAO codes: If set to 1, the band occupations will be determined by `ocp_set`. - - For RT-TDDFT in LCAO codes: If set to 1, same as above, but the occupations will be constrained starting from the second ionic step. + - For PW and LCAO codes: If enabled, the band occupations will be determined by `ocp_set`. + - For RT-TDDFT in LCAO codes: If enabled, same as above, but the occupations will be constrained starting from the second ionic step. - For OFDFT: This feature is not available. - **Description**: - True: Fixes the band occupations based on the values specified in `ocp_set`. diff --git a/source/source_base/matrix3.cpp b/source/source_base/matrix3.cpp index 7fa4398e58..a811d113bb 100644 --- a/source/source_base/matrix3.cpp +++ b/source/source_base/matrix3.cpp @@ -7,179 +7,182 @@ Matrix3::Matrix3(const double &r11, const double &r12, const double &r13, const double &r21, const double &r22, const double &r23, const double &r31, const double &r32, const double &r33) { - e11 = r11;e12 = r12;e13 = r13; - e21 = r21;e22 = r22;e23 = r23; - e31 = r31;e32 = r32;e33 = r33; + e11 = r11; e12 = r12; e13 = r13; + e21 = r21; e22 = r22; e23 = r23; + e31 = r31; e32 = r32; e33 = r33; } void Matrix3::Identity(void) { - e11 = 1;e12 = 0;e13 = 0; - e21 = 0;e22 = 1;e23 = 0; - e31 = 0;e32 = 0;e33 = 1; + e11 = 1; e12 = 0; e13 = 0; + e21 = 0; e22 = 1; e23 = 0; + e31 = 0; e32 = 0; e33 = 1; } void Matrix3::Zero(void) { - e11 = 0;e12 = 0;e13 = 0; - e21 = 0;e22 = 0;e23 = 0; - e31 = 0;e32 = 0;e33 = 0; + e11 = 0; e12 = 0; e13 = 0; + e21 = 0; e22 = 0; e23 = 0; + e31 = 0; e32 = 0; e33 = 0; } -double Matrix3::Det(void) const +double Matrix3::Det(void) const { - return e11*e22*e33 - - e11*e32*e23 + - e21*e32*e13 - - e21*e12*e33 + - e31*e12*e23 - - e31*e22*e13; + return e11*e22*e33 - + e11*e32*e23 + + e21*e32*e13 - + e21*e12*e33 + + e31*e12*e23 - + e31*e22*e13; } Matrix3 Matrix3::Transpose(void) const { - return Matrix3(e11, e21, e31, e12, e22, e32, e13, e23, e33); + return Matrix3(e11, e21, e31, e12, e22, e32, e13, e23, e33); } Matrix3 Matrix3::Inverse(void) const { - double d = this->Det(); + double d = this->Det(); - if(d == 0)d = 1; + if(d == 0) + { + d = 1; + } - return Matrix3((e22*e33 - e23*e32) / d, - -(e12*e33 - e13*e32) / d, - (e12*e23 - e13*e22) / d, - -(e21*e33 - e23*e31) / d, - (e11*e33 - e13*e31) / d, - -(e11*e23 - e13*e21) / d, - (e21*e32 - e22*e31) / d, - -(e11*e32 - e12*e31) / d, - (e11*e22 - e12*e21) / d); + return Matrix3((e22*e33 - e23*e32) / d, + -(e12*e33 - e13*e32) / d, + (e12*e23 - e13*e22) / d, + -(e21*e33 - e23*e31) / d, + (e11*e33 - e13*e31) / d, + -(e11*e23 - e13*e21) / d, + (e21*e32 - e22*e31) / d, + -(e11*e32 - e12*e31) / d, + (e11*e22 - e12*e21) / d); } Matrix3& Matrix3::operator = (const Matrix3 &m) { - e11 = m.e11;e12 = m.e12;e13 = m.e13; - e21 = m.e21;e22 = m.e22;e23 = m.e23; - e31 = m.e31;e32 = m.e32;e33 = m.e33; - return *this; + e11 = m.e11; e12 = m.e12; e13 = m.e13; + e21 = m.e21; e22 = m.e22; e23 = m.e23; + e31 = m.e31; e32 = m.e32; e33 = m.e33; + return *this; } Matrix3& Matrix3::operator +=(const Matrix3 &m) { - e11 += m.e11;e12 += m.e12;e13 += m.e13; - e21 += m.e21;e22 += m.e22;e23 += m.e23; - e31 += m.e31;e32 += m.e32;e33 += m.e33; - return *this; + e11 += m.e11; e12 += m.e12; e13 += m.e13; + e21 += m.e21; e22 += m.e22; e23 += m.e23; + e31 += m.e31; e32 += m.e32; e33 += m.e33; + return *this; } Matrix3& Matrix3::operator -=(const Matrix3 &m) { - e11 -= m.e11;e12 -= m.e12;e13 -= m.e13; - e21 -= m.e21;e22 -= m.e22;e23 -= m.e23; - e31 -= m.e31;e32 -= m.e32;e33 -= m.e33; - return *this; + e11 -= m.e11; e12 -= m.e12; e13 -= m.e13; + e21 -= m.e21; e22 -= m.e22; e23 -= m.e23; + e31 -= m.e31; e32 -= m.e32; e33 -= m.e33; + return *this; } Matrix3& Matrix3::operator *=(const double &s) { - e11 *= s;e12 *= s;e13 *= s; - e21 *= s;e22 *= s;e23 *= s; - e31 *= s;e32 *= s;e33 *= s; - return *this; + e11 *= s; e12 *= s; e13 *= s; + e21 *= s; e22 *= s; e23 *= s; + e31 *= s; e32 *= s; e33 *= s; + return *this; } Matrix3& Matrix3::operator /=(const double &s) { - e11 /= s;e12 /= s;e13 /= s; - e21 /= s;e22 /= s;e23 /= s; - e31 /= s;e32 /= s;e33 /= s; - return *this; + e11 /= s; e12 /= s; e13 /= s; + e21 /= s; e22 /= s; e23 /= s; + e31 /= s; e32 /= s; e33 /= s; + return *this; } //m1+m2 Matrix3 operator +(const Matrix3 &m1, const Matrix3 &m2) { - return Matrix3(m1.e11 + m2.e11, - m1.e12 + m2.e12, - m1.e13 + m2.e13, - m1.e21 + m2.e21, - m1.e22 + m2.e22, - m1.e23 + m2.e23, - m1.e31 + m2.e31, - m1.e32 + m2.e32, - m1.e33 + m2.e33); + return Matrix3(m1.e11 + m2.e11, + m1.e12 + m2.e12, + m1.e13 + m2.e13, + m1.e21 + m2.e21, + m1.e22 + m2.e22, + m1.e23 + m2.e23, + m1.e31 + m2.e31, + m1.e32 + m2.e32, + m1.e33 + m2.e33); } //m1-m2 Matrix3 operator -(const Matrix3 &m1, const Matrix3 &m2) { - return Matrix3(m1.e11 - m2.e11, // Zujian Dai fix bug 2019-01-21 - m1.e12 - m2.e12, - m1.e13 - m2.e13, - m1.e21 - m2.e21, - m1.e22 - m2.e22, - m1.e23 - m2.e23, - m1.e31 - m2.e31, - m1.e32 - m2.e32, - m1.e33 - m2.e33); + return Matrix3(m1.e11 - m2.e11, // Zujian Dai fix bug 2019-01-21 + m1.e12 - m2.e12, + m1.e13 - m2.e13, + m1.e21 - m2.e21, + m1.e22 - m2.e22, + m1.e23 - m2.e23, + m1.e31 - m2.e31, + m1.e32 - m2.e32, + m1.e33 - m2.e33); } //m/s Matrix3 operator /(const Matrix3 &m, const double &s) { - return Matrix3(m.e11 / s, m.e12 / s, m.e13 / s, - m.e21 / s, m.e22 / s, m.e23 / s, - m.e31 / s, m.e32 / s, m.e33 / s); + return Matrix3(m.e11 / s, m.e12 / s, m.e13 / s, + m.e21 / s, m.e22 / s, m.e23 / s, + m.e31 / s, m.e32 / s, m.e33 / s); } //m1*m2 Matrix3 operator *(const Matrix3 &m1, const Matrix3 &m2) { - return Matrix3(m1.e11*m2.e11 + m1.e12*m2.e21 + m1.e13*m2.e31, - m1.e11*m2.e12 + m1.e12*m2.e22 + m1.e13*m2.e32, - m1.e11*m2.e13 + m1.e12*m2.e23 + m1.e13*m2.e33, - m1.e21*m2.e11 + m1.e22*m2.e21 + m1.e23*m2.e31, - m1.e21*m2.e12 + m1.e22*m2.e22 + m1.e23*m2.e32, - m1.e21*m2.e13 + m1.e22*m2.e23 + m1.e23*m2.e33, - m1.e31*m2.e11 + m1.e32*m2.e21 + m1.e33*m2.e31, - m1.e31*m2.e12 + m1.e32*m2.e22 + m1.e33*m2.e32, - m1.e31*m2.e13 + m1.e32*m2.e23 + m1.e33*m2.e33); + return Matrix3(m1.e11*m2.e11 + m1.e12*m2.e21 + m1.e13*m2.e31, + m1.e11*m2.e12 + m1.e12*m2.e22 + m1.e13*m2.e32, + m1.e11*m2.e13 + m1.e12*m2.e23 + m1.e13*m2.e33, + m1.e21*m2.e11 + m1.e22*m2.e21 + m1.e23*m2.e31, + m1.e21*m2.e12 + m1.e22*m2.e22 + m1.e23*m2.e32, + m1.e21*m2.e13 + m1.e22*m2.e23 + m1.e23*m2.e33, + m1.e31*m2.e11 + m1.e32*m2.e21 + m1.e33*m2.e31, + m1.e31*m2.e12 + m1.e32*m2.e22 + m1.e33*m2.e32, + m1.e31*m2.e13 + m1.e32*m2.e23 + m1.e33*m2.e33); } //m*s Matrix3 operator *(const Matrix3 &m,const double &s) { - return Matrix3(m.e11*s, m.e12*s, m.e13*s, - m.e21*s, m.e22*s, m.e23*s, - m.e31*s, m.e32*s, m.e33*s); + return Matrix3(m.e11*s, m.e12*s, m.e13*s, + m.e21*s, m.e22*s, m.e23*s, + m.e31*s, m.e32*s, m.e33*s); } //s*m Matrix3 operator *(const double &s, const Matrix3 &m) { - return Matrix3(m.e11*s, m.e12*s, m.e13*s, - m.e21*s, m.e22*s, m.e23*s, - m.e31*s, m.e32*s, m.e33*s); + return Matrix3(m.e11*s, m.e12*s, m.e13*s, + m.e21*s, m.e22*s, m.e23*s, + m.e31*s, m.e32*s, m.e33*s); } // whether m1==m2 bool operator==(const Matrix3 &m1, const Matrix3 &m2) { - if(m1.e11 == m2.e11 && - m1.e12 == m2.e12 && - m1.e13 == m2.e13 && - m1.e21 == m2.e21 && - m1.e22 == m2.e22 && - m1.e23 == m2.e23 && - m1.e31 == m2.e31 && - m1.e32 == m2.e32 && - m1.e33 == m2.e33) - { - return true; - } - return false; + if(m1.e11 == m2.e11 && + m1.e12 == m2.e12 && + m1.e13 == m2.e13 && + m1.e21 == m2.e21 && + m1.e22 == m2.e22 && + m1.e23 == m2.e23 && + m1.e31 == m2.e31 && + m1.e32 == m2.e32 && + m1.e33 == m2.e33) + { + return true; + } + return false; } //whether m1 != m2 @@ -189,21 +192,21 @@ bool operator!=(const Matrix3 &m1, const Matrix3 &m2) } -void Matrix3::print(void)const +void Matrix3::print(void) const { - std::cout << e11 << std::setw(15) << e12 << std::setw(15) << e13 << std::endl ; - std::cout << e21 << std::setw(15) << e22 << std::setw(15) << e23 << std::endl ; - std::cout << e31 << std::setw(15) << e32 << std::setw(15) << e33 << std::endl ; - return; + std::cout << std::setw(15) << e11 << std::setw(15) << e12 << std::setw(15) << e13 << std::endl; + std::cout << std::setw(15) << e21 << std::setw(15) << e22 << std::setw(15) << e23 << std::endl; + std::cout << std::setw(15) << e31 << std::setw(15) << e32 << std::setw(15) << e33 << std::endl; + return; } -ModuleBase::matrix Matrix3::to_matrix(void)const // Peize Lin add 2021.03.09 +ModuleBase::matrix Matrix3::to_matrix(void)const // Peize Lin add 2021.03.09 { - ModuleBase::matrix m(3,3); - m(0,0)=e11; m(0,1)=e12; m(0,2)=e13; - m(1,0)=e21; m(1,1)=e22; m(1,2)=e23; - m(2,0)=e31; m(2,1)=e32; m(2,2)=e33; - return m; + ModuleBase::matrix m(3,3); + m(0,0)=e11; m(0,1)=e12; m(0,2)=e13; + m(1,0)=e21; m(1,1)=e22; m(1,2)=e23; + m(2,0)=e31; m(2,1)=e32; m(2,2)=e33; + return m; } } \ No newline at end of file diff --git a/source/source_estate/module_pot/H_TDDFT_pw.cpp b/source/source_estate/module_pot/H_TDDFT_pw.cpp index 165f1fe5b4..3f6c8a7fcd 100644 --- a/source/source_estate/module_pot/H_TDDFT_pw.cpp +++ b/source/source_estate/module_pot/H_TDDFT_pw.cpp @@ -3,10 +3,10 @@ #include "source_base/constants.h" #include "source_base/math_integral.h" #include "source_base/timer.h" -#include "source_lcao/module_rt/evolve_elec.h" -#include "source_pw/module_pwdft/global.h" #include "source_io/input_conv.h" #include "source_io/module_parameter/parameter.h" +#include "source_lcao/module_rt/evolve_elec.h" +#include "source_pw/module_pwdft/global.h" namespace elecstate { @@ -15,20 +15,17 @@ int H_TDDFT_pw::istep = -1; bool H_TDDFT_pw::is_initialized = false; double H_TDDFT_pw::amp; -double H_TDDFT_pw::bmod; -double H_TDDFT_pw::bvec[3]; -// Used for calculating electric field force on ions +// Used for calculating electric field force on ions, summing over directions vector H_TDDFT_pw::global_vext_time = {0.0, 0.0, 0.0}; int H_TDDFT_pw::stype; // 0 : length gauge 1: velocity gauge std::vector H_TDDFT_pw::ttype; -// 0 Gauss type function. -// 1 trapezoid type function. -// 2 Trigonometric functions, sin^2. -// 3 heaviside function. -// 4 HHG function. +// 0: Gaussian type function. +// 1: Trapezoid type function. +// 2: Trigonometric functions, sin^2. +// 3: Heaviside step function. int H_TDDFT_pw::tstart; int H_TDDFT_pw::tend; @@ -94,9 +91,9 @@ void H_TDDFT_pw::current_step_info(const std::string& file_dir, int& istep) } file >> istep; - file >> At[0] >> At[1] >>At[2]; + file >> At[0] >> At[1] >> At[2]; file >> At_laststep[0] >> At_laststep[1] >> At_laststep[2]; - At_laststep=-At_laststep; + At_laststep = -At_laststep; file.close(); } @@ -119,7 +116,6 @@ void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo) { return; } - //std::cout << "calculate electric potential" << std::endl; ModuleBase::timer::tick("H_TDDFT_pw", "cal_fixed_v"); @@ -131,20 +127,6 @@ void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo) global_vext_time = {0.0, 0.0, 0.0}; - if (PARAM.inp.td_vext_dire.size() != 1) - { - ModuleBase::WARNING("H_TDDFT_pw::cal_fixed_v", - "Multiple electric fields detected. This feature may have potential issues and is not " - "recommended for use!"); - } - if (PARAM.inp.td_vext_dire.size() > 2) - { - // To avoid breaking the integration test 601_NO_TDDFT_H2_len_hhg, a maximum of 2 electric fields are allowed - ModuleBase::WARNING_QUIT("H_TDDFT_pw::cal_fixed_v", - "For the sake of program stability, the feature of applying multiple electric fields " - "simultaneously has been temporarily disabled. Thank you for your understanding!"); - } - for (auto direc: PARAM.inp.td_vext_dire) { std::vector vext_space(this->rho_basis_->nrxx, 0.0); @@ -155,7 +137,7 @@ void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo) if (PARAM.inp.out_efield && GlobalV::MY_RANK == 0) { std::stringstream as; - as << PARAM.globalv.global_out_dir << "efield_" << count << ".dat"; + as << PARAM.globalv.global_out_dir << "efield_" << count << ".txt"; std::ofstream ofs(as.str().c_str(), std::ofstream::app); ofs << H_TDDFT_pw::istep * dt * ModuleBase::AU_to_FS << "\t" << vext_time * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << std::endl; @@ -198,8 +180,6 @@ void H_TDDFT_pw::cal_v_space_length(std::vector& vext_space, int direc) ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space_length"); ModuleBase::timer::tick("H_TDDFT_pw", "cal_v_space_length"); - prepare(ucell_->G, direc); - for (int ir = 0; ir < this->rho_basis_->nrxx; ++ir) { int i = ir / (this->rho_basis_->ny * this->rho_basis_->nplane); @@ -212,15 +192,21 @@ void H_TDDFT_pw::cal_v_space_length(std::vector& vext_space, int direc) switch (direc) { case 1: - vext_space[ir] = cal_v_space_length_potential(x) / bmod; + vext_space[ir] = cal_v_space_length_potential(x) * this->ucell_->latvec.e11 + + cal_v_space_length_potential(y) * this->ucell_->latvec.e21 + + cal_v_space_length_potential(z) * this->ucell_->latvec.e31; break; case 2: - vext_space[ir] = cal_v_space_length_potential(y) / bmod; + vext_space[ir] = cal_v_space_length_potential(x) * this->ucell_->latvec.e12 + + cal_v_space_length_potential(y) * this->ucell_->latvec.e22 + + cal_v_space_length_potential(z) * this->ucell_->latvec.e32; break; case 3: - vext_space[ir] = cal_v_space_length_potential(z) / bmod; + vext_space[ir] = cal_v_space_length_potential(x) * this->ucell_->latvec.e13 + + cal_v_space_length_potential(y) * this->ucell_->latvec.e23 + + cal_v_space_length_potential(z) * this->ucell_->latvec.e33; break; default: @@ -272,10 +258,6 @@ int H_TDDFT_pw::check_ncut(int t_type) ncut = 2; break; - // case 4: - // vext_time = cal_v_time_HHG(); - // break; - default: std::cout << "time_domain_type of electric field is wrong" << std::endl; break; @@ -285,11 +267,10 @@ int H_TDDFT_pw::check_ncut(int t_type) void H_TDDFT_pw::update_At() { - //std::cout << "calculate electric potential" << std::endl; // time evolve H_TDDFT_pw::istep++; // midpoint rule should be used both in Hamiltonian and here. - At = At + At_laststep/2.0; + At = At + At_laststep / 2.0; At_laststep.set(0.0, 0.0, 0.0); Et.set(0.0, 0.0, 0.0); @@ -341,8 +322,8 @@ void H_TDDFT_pw::update_At() At_laststep[direc - 1] -= out; break; case 2: - At_laststep[direc-1] -= out; - Et[direc-1] += vext_time[0]; + At_laststep[direc - 1] -= out; + Et[direc - 1] += vext_time[0]; break; default: std::cout << "space_domain_type of electric field is wrong" << std::endl; @@ -353,7 +334,7 @@ void H_TDDFT_pw::update_At() if (PARAM.inp.out_efield && GlobalV::MY_RANK == 0) { std::stringstream as; - as << PARAM.globalv.global_out_dir << "efield_" << count << ".dat"; + as << PARAM.globalv.global_out_dir << "efield_" << count << ".txt"; std::ofstream ofs(as.str().c_str(), std::ofstream::app); ofs << H_TDDFT_pw::istep * dt * ModuleBase::AU_to_FS << "\t" << vext_time[0] * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << std::endl; @@ -362,7 +343,7 @@ void H_TDDFT_pw::update_At() // total count++ count++; } - At = At + At_laststep/2.0; + At = At + At_laststep / 2.0; ModuleBase::timer::tick("H_TDDFT_pw", "update_At"); return; @@ -390,10 +371,6 @@ double H_TDDFT_pw::cal_v_time(int t_type, const bool last) vext_time = cal_v_time_heaviside(last); break; - // case 4: - // vext_time = cal_v_time_HHG(); - // break; - default: std::cout << "time_domain_type of electric field is wrong" << std::endl; break; @@ -487,36 +464,12 @@ double H_TDDFT_pw::cal_v_time_heaviside(const bool last) { vext_time = 0.0; } - if(last)heavi_count++; - - return vext_time; -} - -void H_TDDFT_pw::prepare(const ModuleBase::Matrix3& G, int& dir) -{ - if (dir == 1) - { - bvec[0] = G.e11; - bvec[1] = G.e12; - bvec[2] = G.e13; - } - else if (dir == 2) - { - bvec[0] = G.e21; - bvec[1] = G.e22; - bvec[2] = G.e23; - } - else if (dir == 3) - { - bvec[0] = G.e31; - bvec[1] = G.e32; - bvec[2] = G.e33; - } - else + if (last) { - ModuleBase::WARNING_QUIT("H_TDDFT_pw::prepare", "direction is wrong!"); + heavi_count++; } - bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2)); + + return vext_time; } void H_TDDFT_pw::compute_force(const UnitCell& cell, ModuleBase::matrix& fe) @@ -526,16 +479,14 @@ void H_TDDFT_pw::compute_force(const UnitCell& cell, ModuleBase::matrix& fe) { for (int ia = 0; ia < cell.atoms[it].na; ++ia) { - for (int jj = 0; jj < 3; ++jj) + for (int direc = 0; direc < 3; ++direc) { // No need to multiply ModuleBase::e2, since the unit of force is Ry/Bohr - fe(iat, jj) - = (std::abs(bmod) > 1e-10 ? global_vext_time[jj] * cell.atoms[it].ncpp.zv * bvec[jj] / bmod : 0); + fe(iat, direc) = global_vext_time[direc] * cell.atoms[it].ncpp.zv; } ++iat; } } } -} // namespace elecstate - +} // namespace elecstate \ No newline at end of file diff --git a/source/source_estate/module_pot/H_TDDFT_pw.h b/source/source_estate/module_pot/H_TDDFT_pw.h index 6c69d63566..8f2064eb35 100644 --- a/source/source_estate/module_pot/H_TDDFT_pw.h +++ b/source/source_estate/module_pot/H_TDDFT_pw.h @@ -1,9 +1,9 @@ #ifndef H_TDDFT_PW_H #define H_TDDFT_PW_H +#include "pot_base.h" #include "source_io/input_conv.h" #include "source_io/module_parameter/parameter.h" // PARAM.globalv.global_readin_dir, PARAM.inp.mdp.md_restart -#include "pot_base.h" namespace elecstate { @@ -40,22 +40,21 @@ class H_TDDFT_pw : public PotBase void cal_fixed_v(double* vl_pseudo) override; /** - * @brief compute force of electric field + * @brief Compute ionic force of electric field * - * @param[in] cell information of cell - * @param[out] fe force of electric field F=qE + * @param[in] cell Information of cell + * @param[out] fe Force of electric field F = qE */ static void compute_force(const UnitCell& cell, ModuleBase::matrix& fe); // parameters - static int stype; // 0 : length gauge 1: velocity gauge + static int stype; // 0: length gauge; 1: velocity gauge; 2: hybrid gauge static std::vector ttype; - // 0 Gauss type function. - // 1 trapezoid type function. - // 2 Trigonometric functions, sin^2. - // 3 heaviside function. - // 4 HHG function. + // 0: Gaussian type function. + // 1: Trapezoid type function. + // 2: Trigonometric functions, sin^2. + // 3: Heaviside step function. static int tstart; static int tend; @@ -64,18 +63,18 @@ class H_TDDFT_pw : public PotBase static double dt_int; static int istep_int; - // space domain parameters + // Space domain parameters // length gauge static double lcut1; static double lcut2; - // velocity gauge, vector magnetic potential + // velocity gauge, vector potential static ModuleBase::Vector3 At; static ModuleBase::Vector3 At_laststep; static ModuleBase::Vector3 Et; - // time domain parameters + // Time domain parameters // Gauss static int gauss_count; @@ -84,11 +83,11 @@ class H_TDDFT_pw : public PotBase static std::vector gauss_sigma; // time(a.u.) static std::vector gauss_t0; static std::vector gauss_amp; // Ry/bohr - // add for velocity gauge, recut dt into n pieces to make sure the integral is accurate Enough - // must be even, thus would get odd number of points for simpson integral + // add for velocity gauge, recut dt into n pieces to make sure the integral is accurate enough + // must be even, thus would get odd number of points for Simpson integral static std::vector gauss_ncut; - // trapezoid + // Trapezoid static int trape_count; static std::vector trape_omega; // time(a.u.)^-1 static std::vector trape_phase; @@ -96,7 +95,7 @@ class H_TDDFT_pw : public PotBase static std::vector trape_t2; static std::vector trape_t3; static std::vector trape_amp; // Ry/bohr - // add for velocity gauge, recut dt into n pieces to make sure the integral is accurate Enough + // add for velocity gauge, recut dt into n pieces to make sure the integral is accurate enough static std::vector trape_ncut; // Trigonometric @@ -106,7 +105,7 @@ class H_TDDFT_pw : public PotBase static std::vector trigo_phase1; static std::vector trigo_phase2; static std::vector trigo_amp; // Ry/bohr - // add for velocity gauge, recut dt into n pieces to make sure the integral is accurate Enough + // add for velocity gauge, recut dt into n pieces to make sure the integral is accurate enough static std::vector trigo_ncut; // Heaviside @@ -118,44 +117,33 @@ class H_TDDFT_pw : public PotBase static void update_At(); private: - // internal time-step, - //-------hypothesis------- - // Vext will evolve by time, every time cal_fixed_v() is called, istep++ - //------------------------ static int istep; static bool is_initialized; // static flag variable, used to ensure initialization only once static double amp; static vector global_vext_time; - static double bmod; - static double bvec[3]; - const UnitCell* ucell_ = nullptr; // Obtain the current MD step information, used for restart calculation void current_step_info(const std::string& file_dir, int& istep); - // potential of electric field in space domain : length gauge and velocity gauge + // Potential of electric field in space domain: for length gauge only void cal_v_space(std::vector& vext_space, int direc); void cal_v_space_length(std::vector& vext_space, int direc); double cal_v_space_length_potential(double i); - // potential of electric field in time domain : Gauss , trapezoid, trigonometric, heaviside, HHG + // Potential of electric field in time domain: Gaussian, trapezoid, trigonometric, Heaviside static double cal_v_time(int t_type, const bool last); static double cal_v_time_Gauss(const bool last); static double cal_v_time_trapezoid(const bool last); static double cal_v_time_trigonometric(const bool last); static double cal_v_time_heaviside(const bool last); - // double cal_v_time_HHG(); - // get ncut number for At integral + // Get ncut number for At integral static int check_ncut(int t_type); - - void prepare(const ModuleBase::Matrix3& G, int& dir); }; } // namespace elecstate #endif - diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 2222cd4c19..7c38b47372 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -297,11 +297,10 @@ struct Input_para int propagator = 0; ///< method of propagator int td_stype = 0; ///< type of space domain 0 : length gauge 1: velocity gauge std::string td_ttype = "0"; ///< type of time domain - ///< 0 Gauss type function. - ///< 1 trapezoid type function. - ///< 2 Trigonometric functions, sin^2. - ///< 3 heaviside function. - ///< 4 HHG function. + ///< 0: Gaussian type function. + ///< 1: Trapezoid type function. + ///< 2: Trigonometric functions, sin^2. + ///< 3: Heaviside step function. int td_tstart = 1; int td_tend = 1000; diff --git a/source/source_io/test/for_testing_input_conv.h b/source/source_io/test/for_testing_input_conv.h index 29fcac72d0..8cca07361f 100644 --- a/source/source_io/test/for_testing_input_conv.h +++ b/source/source_io/test/for_testing_input_conv.h @@ -45,34 +45,24 @@ double elecstate::Efield::efield_pos_max; double elecstate::Efield::efield_pos_dec; double elecstate::Efield::efield_amp; -// parameters of electric field for tddft +// Parameters of electric field for RT-TDDFT -int elecstate::H_TDDFT_pw::stype; // 0 : length gauge 1: velocity gauge +int elecstate::H_TDDFT_pw::stype; std::vector elecstate::H_TDDFT_pw::ttype; -// 0 Gauss type function. -// 1 trapezoid type function. -// 2 Trigonometric functions, sin^2. -// 3 heaviside function. -// 4 HHG function. int elecstate::H_TDDFT_pw::tstart; int elecstate::H_TDDFT_pw::tend; double elecstate::H_TDDFT_pw::dt; double elecstate::H_TDDFT_pw::dt_int; -// space domain parameters - -// length gauge double elecstate::H_TDDFT_pw::lcut1; double elecstate::H_TDDFT_pw::lcut2; -// time domain parameters - bool TD_Velocity::tddft_velocity; bool TD_Velocity::out_mat_R; -// Gauss +// Gaussian int elecstate::H_TDDFT_pw::gauss_count; std::vector elecstate::H_TDDFT_pw::gauss_omega; // time(a.u.)^-1 std::vector elecstate::H_TDDFT_pw::gauss_phase; @@ -81,7 +71,7 @@ std::vector elecstate::H_TDDFT_pw::gauss_t0; std::vector elecstate::H_TDDFT_pw::gauss_amp; // Ry/bohr std::vector elecstate::H_TDDFT_pw::gauss_ncut; -// trapezoid +// Trapezoid int elecstate::H_TDDFT_pw::trape_count; std::vector elecstate::H_TDDFT_pw::trape_omega; // time(a.u.)^-1 std::vector elecstate::H_TDDFT_pw::trape_phase; diff --git a/source/source_lcao/module_rt/band_energy.cpp b/source/source_lcao/module_rt/band_energy.cpp index 8767d37405..63b2fad71b 100644 --- a/source/source_lcao/module_rt/band_energy.cpp +++ b/source/source_lcao/module_rt/band_energy.cpp @@ -96,9 +96,15 @@ void compute_ekb(const Parallel_Orbitals* pv, { bb = 0.0; } - if (aa > 0.0 || bb > 0.0) + if (std::abs(aa) > 0.0 || std::abs(bb) > 0.0) { - ofs_running << i << " " << j << " " << aa << "+" << bb << "i " << std::endl; + std::streamsize original_precision = ofs_running.precision(); + ofs_running << std::fixed << std::setprecision(8); + ofs_running << "i = " << std::setw(2) << i << ", j = " << std::setw(2) << j + << ", Eij = " << std::setw(12) << aa << " + " << std::setw(12) << bb << " i" + << std::endl; + ofs_running.unsetf(std::ios_base::fixed); + ofs_running.precision(original_precision); } } } @@ -233,9 +239,15 @@ void compute_ekb_tensor(const Parallel_Orbitals* pv, { bb = 0.0; } - if (aa > 0.0 || bb > 0.0) + if (std::abs(aa) > 0.0 || std::abs(bb) > 0.0) { - ofs_running << i << " " << j << " " << aa << "+" << bb << "i " << std::endl; + std::streamsize original_precision = ofs_running.precision(); + ofs_running << std::fixed << std::setprecision(8); + ofs_running << "i = " << std::setw(2) << i << ", j = " << std::setw(2) << j + << ", Eij = " << std::setw(12) << aa << " + " << std::setw(12) << bb << " i" + << std::endl; + ofs_running.unsetf(std::ios_base::fixed); + ofs_running.precision(original_precision); } } } @@ -371,9 +383,15 @@ void compute_ekb_tensor_lapack(const Parallel_Orbitals* pv, { bb = 0.0; } - if (aa > 0.0 || bb > 0.0) + if (std::abs(aa) > 0.0 || std::abs(bb) > 0.0) { - ofs_running << i << " " << j << " " << aa << "+" << bb << "i " << std::endl; + std::streamsize original_precision = ofs_running.precision(); + ofs_running << std::fixed << std::setprecision(8); + ofs_running << "i = " << std::setw(2) << i << ", j = " << std::setw(2) << j + << ", Eij = " << std::setw(12) << aa << " + " << std::setw(12) << bb << " i" + << std::endl; + ofs_running.unsetf(std::ios_base::fixed); + ofs_running.precision(original_precision); } } }