diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index b0a33902d6..f709710d33 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp @@ -1,6 +1,7 @@ #include "H_TDDFT_pw.h" -#include "module_base/math_integral.h" + #include "module_base/constants.h" +#include "module_base/math_integral.h" #include "module_base/timer.h" #include "module_hamilt_lcao/module_tddft/evolve_elec.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" @@ -27,7 +28,7 @@ std::vector H_TDDFT_pw::ttype; int H_TDDFT_pw::tstart; int H_TDDFT_pw::tend; double H_TDDFT_pw::dt; -//cut dt for integral +// cut dt for integral double H_TDDFT_pw::dt_int; int H_TDDFT_pw::istep_int; // space domain parameters @@ -36,7 +37,7 @@ int H_TDDFT_pw::istep_int; double H_TDDFT_pw::lcut1; double H_TDDFT_pw::lcut2; -//velocity gauge +// velocity gauge ModuleBase::Vector3 H_TDDFT_pw::At; // time domain parameters @@ -48,7 +49,7 @@ std::vector H_TDDFT_pw::gauss_phase; std::vector H_TDDFT_pw::gauss_sigma; // time(a.u.) std::vector H_TDDFT_pw::gauss_t0; std::vector H_TDDFT_pw::gauss_amp; // Ry/bohr -std::vector H_TDDFT_pw::gauss_ncut; //cut for integral +std::vector H_TDDFT_pw::gauss_ncut; // cut for integral // trapezoid int H_TDDFT_pw::trape_count; @@ -58,7 +59,7 @@ std::vector H_TDDFT_pw::trape_t1; std::vector H_TDDFT_pw::trape_t2; std::vector H_TDDFT_pw::trape_t3; std::vector H_TDDFT_pw::trape_amp; // Ry/bohr -std::vector H_TDDFT_pw::trape_ncut; //cut for integral +std::vector H_TDDFT_pw::trape_ncut; // cut for integral // Trigonometric int H_TDDFT_pw::trigo_count; @@ -67,20 +68,22 @@ std::vector H_TDDFT_pw::trigo_omega2; // time(a.u.)^-1 std::vector H_TDDFT_pw::trigo_phase1; std::vector H_TDDFT_pw::trigo_phase2; std::vector H_TDDFT_pw::trigo_amp; // Ry/bohr -std::vector H_TDDFT_pw::trigo_ncut; //cut for integral +std::vector H_TDDFT_pw::trigo_ncut; // cut for integral // Heaviside int H_TDDFT_pw::heavi_count; std::vector H_TDDFT_pw::heavi_t0; std::vector H_TDDFT_pw::heavi_amp; // Ry/bohr -void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) +void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo) { ModuleBase::TITLE("H_TDDFT_pw", "cal_fixed_v"); - //skip if velocity_gague - if(stype==1) {return; -} + // skip if velocity_gague + if (stype == 1) + { + return; + } // time evolve H_TDDFT_pw::istep++; @@ -104,7 +107,7 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) for (auto direc: module_tddft::Evolve_elec::td_vext_dire_case) { std::vector vext_space(this->rho_basis_->nrxx, 0.0); - double vext_time = cal_v_time(ttype[count],true); + double vext_time = cal_v_time(ttype[count], true); if (module_tddft::Evolve_elec::out_efield && GlobalV::MY_RANK == 0) { @@ -118,7 +121,9 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) cal_v_space(vext_space, direc); for (size_t ir = 0; ir < this->rho_basis_->nrxx; ++ir) + { vl_pseudo[ir] += vext_space[ir] * vext_time; + } count++; } @@ -126,7 +131,7 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) return; } -void H_TDDFT_pw::cal_v_space(std::vector &vext_space, int direc) +void H_TDDFT_pw::cal_v_space(std::vector& vext_space, int direc) { ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space"); ModuleBase::timer::tick("H_TDDFT_pw", "cal_v_space"); @@ -145,7 +150,7 @@ void H_TDDFT_pw::cal_v_space(std::vector &vext_space, int direc) return; } -void H_TDDFT_pw::cal_v_space_length(std::vector &vext_space, int direc) +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"); @@ -190,21 +195,22 @@ double H_TDDFT_pw::cal_v_space_length_potential(double i) double vext_space = 0.0; if (i < lcut1) { - vext_space = ((i - lcut1) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut1) * this->ucell_->lat0; + vext_space = -((i - lcut1) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut1) * this->ucell_->lat0; } else if (i >= lcut1 && i < lcut2) { - vext_space = -i * this->ucell_->lat0; + vext_space = i * this->ucell_->lat0; } else if (i >= lcut2) { - vext_space = ((i - lcut2) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut2) * this->ucell_->lat0; + vext_space = -((i - lcut2) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut2) * this->ucell_->lat0; } return vext_space; } + int H_TDDFT_pw::check_ncut(int t_type) { - int ncut=0; + int ncut = 0; switch (t_type) { case 0: @@ -233,12 +239,13 @@ int H_TDDFT_pw::check_ncut(int t_type) } return ncut; } + void H_TDDFT_pw::update_At() { std::cout << "calculate electric potential" << std::endl; // time evolve H_TDDFT_pw::istep++; - + // judgement to skip vext if (!module_tddft::Evolve_elec::td_vext || istep > tend || istep < tstart) { @@ -249,7 +256,7 @@ void H_TDDFT_pw::update_At() trape_count = 0; trigo_count = 0; heavi_count = 0; - //parameters for integral + // parameters for integral int ncut = 1; bool last = false; double out = 0.0; @@ -257,33 +264,38 @@ void H_TDDFT_pw::update_At() for (auto direc: module_tddft::Evolve_elec::td_vext_dire_case) { last = false; - //cut the integral space and initial relvant parameters + // cut the integral space and initialize relevant parameters ncut = check_ncut(ttype[count]); - istep_int = istep*ncut; - dt_int = dt/double(ncut); + istep_int = istep * ncut; + dt_int = dt / double(ncut); - //store vext_time for each time point, include the first and last point - double* vext_time = new double[ncut+1](); - for(int i=0; i<=ncut; i++) + // store vext_time for each time point, include the first and last point + std::vector vext_time(ncut + 1, 0.0); // Use std::vector to manage memory + for (int i = 0; i <= ncut; i++) { - //if this is the last point, type_count++ - if(i==ncut)last=true; - vext_time[i] = cal_v_time(ttype[count],last); + // if this is the last point, type_count++ + if (i == ncut) + { + last = true; + } + vext_time[i] = cal_v_time(ttype[count], last); istep_int++; } - ModuleBase::Integral::Simpson_Integral(ncut+1, vext_time, dt_int, out); - //update At value for its direction + // Call the Simpson's rule integration using std::vector data + ModuleBase::Integral::Simpson_Integral(ncut + 1, vext_time.data(), dt_int, out); + + // update At value for its direction switch (stype) { case 1: - At[direc-1] -= out; + At[direc - 1] -= out; break; default: std::cout << "space_domain_type of electric field is wrong" << std::endl; break; } - //output Efield + // output Efield if (module_tddft::Evolve_elec::out_efield && GlobalV::MY_RANK == 0) { std::stringstream as; @@ -293,8 +305,7 @@ void H_TDDFT_pw::update_At() << vext_time[0] * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << std::endl; ofs.close(); } - //total count++ - delete[] vext_time; + // total count++ count++; } return; @@ -345,8 +356,10 @@ double H_TDDFT_pw::cal_v_time_Gauss(const bool last) double gauss_t = (istep_int - t0 * ncut) * dt_int; vext_time = cos(omega * gauss_t + phase) * exp(-gauss_t * gauss_t * 0.5 / (sigma * sigma)) * amp; - if(last) {gauss_count++; -} + if (last) + { + gauss_count++; + } return vext_time; } @@ -376,8 +389,10 @@ double H_TDDFT_pw::cal_v_time_trapezoid(const bool last) } vext_time = vext_time * amp * cos(omega * istep_int * dt_int + phase); - if(last) {trape_count++; -} + if (last) + { + trape_count++; + } return vext_time; } @@ -394,8 +409,10 @@ double H_TDDFT_pw::cal_v_time_trigonometric(const bool last) const double timenow = istep_int * dt_int; vext_time = amp * cos(omega1 * timenow + phase1) * sin(omega2 * timenow + phase2) * sin(omega2 * timenow + phase2); - if(last) {trigo_count++; -} + if (last) + { + trigo_count++; + } return vext_time; } @@ -405,11 +422,14 @@ double H_TDDFT_pw::cal_v_time_heaviside() double t0 = *(heavi_t0.begin() + heavi_count); double amp = *(heavi_amp.begin() + heavi_count); double vext_time = 0.0; - if (istep < t0) { + if (istep < t0) + { vext_time = amp; - } else if (istep >= t0) { + } + else if (istep >= t0) + { vext_time = 0.0; -} + } heavi_count++; return vext_time; diff --git a/source/module_io/write_dipole.cpp b/source/module_io/write_dipole.cpp index eb1ca6188c..57794e45ed 100644 --- a/source/module_io/write_dipole.cpp +++ b/source/module_io/write_dipole.cpp @@ -37,21 +37,26 @@ void ModuleIO::write_dipole(const double* rho_save, #ifndef __MPI double dipole_elec_x = 0.0, dipole_elec_y = 0.0, dipole_elec_z = 0.0; + double lat_factor_x = GlobalC::ucell.lat0 * 0.529177 / rhopw->nx; + double lat_factor_y = GlobalC::ucell.lat0 * 0.529177 / rhopw->ny; + double lat_factor_z = GlobalC::ucell.lat0 * 0.529177 / rhopw->nz; + for (int k = 0; k < rhopw->nz; k++) { for (int j = 0; j < rhopw->ny; j++) { for (int i = 0; i < rhopw->nx; i++) { - dipole_elec_x += rho_save[i * rhopw->ny * rhopw->nz + j * rhopw->nz + k] * i - * GlobalC::ucell.lat0 * 0.529177 / rhopw->nx; - dipole_elec_y += rho_save[i * rhopw->ny * rhopw->nz + j * rhopw->nz + k] * j - * GlobalC::ucell.lat0 * 0.529177 / rhopw->ny; - dipole_elec_z += rho_save[i * rhopw->ny * rhopw->nz + j * rhopw->nz + k] * k - * GlobalC::ucell.lat0 * 0.529177 / rhopw->nz; + int index = i * rhopw->ny * rhopw->nz + j * rhopw->nz + k; + double rho_val = rho_save[index]; + + dipole_elec_x += rho_val * i * lat_factor_x; + dipole_elec_y += rho_val * j * lat_factor_y; + dipole_elec_z += rho_val * k * lat_factor_z; } } } + dipole_elec_x *= GlobalC::ucell.omega / static_cast(rhopw->nxyz); dipole_elec_y *= GlobalC::ucell.omega / static_cast(rhopw->nxyz); dipole_elec_z *= GlobalC::ucell.omega / static_cast(rhopw->nxyz); @@ -103,7 +108,6 @@ void ModuleIO::write_dipole(const double* rho_save, double sum = 0; for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ++ia) { - sum += GlobalC::ucell.atoms[it].taud[ia][i]; } dipole_ion[i] += sum * GlobalC::ucell.atoms[it].ncpp.zv; @@ -138,7 +142,7 @@ void ModuleIO::write_dipole(const double* rho_save, return; } -double ModuleIO::prepare(const UnitCell &cell, int &dir) +double ModuleIO::prepare(const UnitCell& cell, int& dir) { double bvec[3] = {0.0}; double bmod = 0.0; diff --git a/tests/integrate/601_NO_TDDFT_H2_len_current/refcurrent_total.dat b/tests/integrate/601_NO_TDDFT_H2_len_current/refcurrent_total.dat index 487e061c1e..5f225fd90f 100644 --- a/tests/integrate/601_NO_TDDFT_H2_len_current/refcurrent_total.dat +++ b/tests/integrate/601_NO_TDDFT_H2_len_current/refcurrent_total.dat @@ -1,4 +1,4 @@ 0 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -2 4.2068003950947370e-04 -3.6201227751435700e-05 -5.7838914079840259e-05 -3 -1.7465510685220912e-04 -9.7934211131431987e-05 -1.5648656468267278e-04 +2 -4.2482488063038011e-04 3.6208093141272442e-05 5.7880862345902912e-05 +3 1.7506474069499989e-04 9.7936761289276890e-05 1.5650847511644919e-04 \ No newline at end of file diff --git a/tests/integrate/601_NO_TDDFT_H2_len_current/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_current/result.ref index 9d4f28fe1a..dbe64306cd 100644 --- a/tests/integrate/601_NO_TDDFT_H2_len_current/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_len_current/result.ref @@ -1,6 +1,6 @@ -etotref -30.91255703869593 -etotperatomref -15.4562785193 -totalforceref 0.459640 -totalstressref 0.945730 +etotref -30.91255706614114 +etotperatomref -15.4562785331 +totalforceref 0.459938 +totalstressref 0.946457 CompareCurrent_pass 0 totaltimeref 1.72 diff --git a/tests/integrate/601_NO_TDDFT_H2_len_gauss/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_gauss/result.ref index 8c19628585..fddba87a4f 100644 --- a/tests/integrate/601_NO_TDDFT_H2_len_gauss/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_len_gauss/result.ref @@ -1,5 +1,5 @@ -etotref -31.36940219279552 -etotperatomref -15.6847010964 -totalforceref 0.574354 -totalstressref 5.532849 +etotref -30.45477750879012 +etotperatomref -15.2273887544 +totalforceref 0.574814 +totalstressref 4.287440 totaltimeref 1.57 diff --git a/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/result.ref index f189567a70..5d2e57b6ef 100644 --- a/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/result.ref @@ -1,5 +1,5 @@ -etotref -31.36928614218948 -etotperatomref -15.6846430711 -totalforceref 0.572958 -totalstressref 5.533528 +etotref -30.45466274119030 +etotperatomref -15.2273313706 +totalforceref 0.573436 +totalstressref 4.285714 totaltimeref 1.55 diff --git a/tests/integrate/601_NO_TDDFT_H2_len_heavi/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_heavi/result.ref index 8213942bfa..de8d31b7db 100644 --- a/tests/integrate/601_NO_TDDFT_H2_len_heavi/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_len_heavi/result.ref @@ -1,5 +1,5 @@ -etotref -38.85688376698273 -etotperatomref -19.4284418835 -totalforceref 0.564906 -totalstressref 77.567005 +etotref -22.97718184915547 +etotperatomref -11.4885909246 +totalforceref 0.564908 +totalstressref 76.073961 totaltimeref 1.56 diff --git a/tests/integrate/601_NO_TDDFT_H2_len_hhg/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_hhg/result.ref index a80a55e0f9..70aa5fdf09 100644 --- a/tests/integrate/601_NO_TDDFT_H2_len_hhg/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_len_hhg/result.ref @@ -1,5 +1,5 @@ -etotref -82.62629897647498 -etotperatomref -41.3131494882 -totalforceref 0.748884 -totalstressref 501.373588 +etotref 20.38731422853930 +etotperatomref 10.1936571143 +totalforceref 0.738256 +totalstressref 489.377111 totaltimeref 1.56 diff --git a/tests/integrate/601_NO_TDDFT_H2_len_trape/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_trape/result.ref index df5e8f25da..b8cbe20aa0 100644 --- a/tests/integrate/601_NO_TDDFT_H2_len_trape/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_len_trape/result.ref @@ -1,5 +1,5 @@ -etotref -30.92452072277197 -etotperatomref -15.4622603614 -totalforceref 0.574664 -totalstressref 1.256903 +etotref -30.89965242483882 +etotperatomref -15.4498262124 +totalforceref 0.574662 +totalstressref 1.034216 totaltimeref 1.44 diff --git a/tests/integrate/601_NO_TDDFT_H2_len_trigo/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_trigo/result.ref index bfb7a3007e..b0cdd8ce31 100644 --- a/tests/integrate/601_NO_TDDFT_H2_len_trigo/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_len_trigo/result.ref @@ -1,5 +1,5 @@ -etotref -30.91505170216714 -etotperatomref -15.4575258511 -totalforceref 0.574664 -totalstressref 1.165899 +etotref -30.90912142263002 +etotperatomref -15.4545607113 +totalforceref 0.574662 +totalstressref 1.108904 totaltimeref 1.43