From 67d1fd7a833709fe9388b9ccaee80cc3c8a6dcd4 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Thu, 2 Mar 2023 15:12:37 +0800 Subject: [PATCH 1/6] fix bug of electric field in tddft --- .../potentials/H_TDDFT_pw.cpp | 4 + source/module_io/dipole_io.h | 13 +- source/module_io/write_dipole.cpp | 410 +++--------------- .../601_NO_TDDFT_H2_len_trape/result.ref | 8 +- 4 files changed, 77 insertions(+), 358 deletions(-) diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index 6f7502d7e2..d817dd58b6 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp @@ -229,6 +229,10 @@ double H_TDDFT_pw::cal_v_time(int t_type) vext_time = cal_v_time_Gauss(); break; + case 1: + vext_time = cal_v_time_trapezoid(); + break; + case 2: vext_time = cal_v_time_trigonometric(); break; diff --git a/source/module_io/dipole_io.h b/source/module_io/dipole_io.h index fbd3d8f0a2..4817f0e3c1 100644 --- a/source/module_io/dipole_io.h +++ b/source/module_io/dipole_io.h @@ -6,11 +6,14 @@ namespace ModuleIO { void write_dipole(const double *rho_save, - const int &is, - const int &istep, - const std::string &fn, - const int &precision = 11, - const bool for_plot = false); + const int &is, + const int &istep, + const std::string &fn, + const int &precision = 11, + const bool for_plot = false); + +double prepare(const UnitCell &cell, int &dir); + } // namespace ModuleIO #endif diff --git a/source/module_io/write_dipole.cpp b/source/module_io/write_dipole.cpp index e6e76a2f10..383268512d 100644 --- a/source/module_io/write_dipole.cpp +++ b/source/module_io/write_dipole.cpp @@ -28,6 +28,12 @@ void ModuleIO::write_dipole(const double *rho_save, } } + double bmod[3]; + for (int i = 0; i < 3; i++) + { + bmod[i] = prepare(GlobalC::ucell, i); + } + #ifndef __MPI double dipole_elec_x = 0.0, dipole_elec_y = 0.0, dipole_elec_z = 0.0; for (int k = 0; k < GlobalC::rhopw->nz; k++) @@ -54,258 +60,6 @@ void ModuleIO::write_dipole(const double *rho_save, ofs << istep << " " << dipole_elec_x << " " << dipole_elec_y << dipole_elec_z; #else - // double dipole_elec_x = 0.0, dipole_elec_y = 0.0, dipole_elec_z = 0.0; - - // // only do in the first pool. - // if (GlobalV::MY_POOL == 0) - // { - // // num_z: how many planes on processor 'ip' - // int *num_z = new int[GlobalV::NPROC_IN_POOL]; - // ModuleBase::GlobalFunc::ZEROS(num_z, GlobalV::NPROC_IN_POOL); - // for (int iz = 0; iz < GlobalC::bigpw->nbz; iz++) - // { - // int ip = iz % GlobalV::NPROC_IN_POOL; - // num_z[ip] += GlobalC::bigpw->bz; - // } - - // // start_z: start position of z in - // // processor ip. - // int *start_z = new int[GlobalV::NPROC_IN_POOL]; - // ModuleBase::GlobalFunc::ZEROS(start_z, GlobalV::NPROC_IN_POOL); - // for (int ip = 1; ip < GlobalV::NPROC_IN_POOL; ip++) - // { - // start_z[ip] = start_z[ip - 1] + num_z[ip - 1]; - // } - - // // which_ip: found iz belongs to which ip. - // int *which_ip = new int[GlobalC::rhopw->nz]; - // ModuleBase::GlobalFunc::ZEROS(which_ip, GlobalC::rhopw->nz); - // for (int iz = 0; iz < GlobalC::rhopw->nz; iz++) - // { - // for (int ip = 0; ip < GlobalV::NPROC_IN_POOL; ip++) - // { - // if (iz >= start_z[GlobalV::NPROC_IN_POOL - 1]) - // { - // which_ip[iz] = GlobalV::NPROC_IN_POOL - 1; - // break; - // } - // else if (iz >= start_z[ip] && iz < start_z[ip + 1]) - // { - // which_ip[iz] = ip; - // break; - // } - // } - // // GlobalV::ofs_running << "\n iz=" << iz << " ip=" << which_ip[iz]; - // } - - // // int count=0; - // int nxy = GlobalC::rhopw->nx * GlobalC::rhopw->ny; - // double *zpiece = new double[nxy]; - - // // save the rho one z by one z. - // for (int iz = 0; iz < GlobalC::rhopw->nz; iz++) - // { - // // std::cout << "\n iz=" << iz << std::endl; - // // tag must be different for different iz. - // ModuleBase::GlobalFunc::ZEROS(zpiece, nxy); - // int tag = iz; - // MPI_Status ierror; - - // // case 1: the first part of rho in processor 0. - // if (which_ip[iz] == 0 && GlobalV::RANK_IN_POOL == 0) - // { - // for (int ir = 0; ir < nxy; ir++) - // { - // // mohan change to rho_save on 2012-02-10 - // // because this can make our next restart calculation lead - // // to the same scf_thr as the one saved. - // zpiece[ir] = rho_save[ir * GlobalC::rhopw->nplane + iz - GlobalC::rhopw->startz_current]; - // // GlobalV::ofs_running << "\n get zpiece[" << ir << "]=" << zpiece[ir] << " - // // ir*GlobalC::pw.nczp+iz=" << ir*GlobalC::pw.nczp+iz; - // } - // } - // // case 2: > first part rho: send the rho to - // // processor 0. - // else if (which_ip[iz] == GlobalV::RANK_IN_POOL) - // { - // for (int ir = 0; ir < nxy; ir++) - // { - // // zpiece[ir] = rho[is][ir*num_z[GlobalV::RANK_IN_POOL]+iz]; - // zpiece[ir] = rho_save[ir * GlobalC::rhopw->nplane + iz - GlobalC::rhopw->startz_current]; - // // GlobalV::ofs_running << "\n get zpiece[" << ir << "]=" << zpiece[ir] << " - // // ir*GlobalC::pw.nczp+iz=" << ir*GlobalC::pw.nczp+iz; - // } - // MPI_Send(zpiece, nxy, MPI_DOUBLE, 0, tag, POOL_WORLD); - // } - - // // case 2: > first part rho: processor 0 receive the rho - // // from other processors - // else if (GlobalV::RANK_IN_POOL == 0) - // { - // MPI_Recv(zpiece, nxy, MPI_DOUBLE, which_ip[iz], tag, POOL_WORLD, &ierror); - // // GlobalV::ofs_running << "\n Receieve First number = " << zpiece[0]; - // } - - // // write data - // if (GlobalV::MY_RANK == 0) - // { - // // ofs << "\niz=" << iz; - // // mohan update 2011-03-30 - // for (int iy = 0; iy < GlobalC::rhopw->ny; iy++) - // { - // for (int ix = 0; ix < GlobalC::rhopw->nx; ix++) - // { - // /* - // if(ixny + iy] * ix * GlobalC::ucell.lat0 * 0.529177 - // / GlobalC::rhopw->nx; - // dipole_elec_y += zpiece[ix * GlobalC::rhopw->ny + iy] * iy * GlobalC::ucell.lat0 * 0.529177 - // / GlobalC::rhopw->ny; - // dipole_elec_z += zpiece[ix * GlobalC::rhopw->ny + iy] * iz * GlobalC::ucell.lat0 * 0.529177 - // / GlobalC::rhopw->nz; - // } - // } - // } - // } // end iz - - // delete[] zpiece; - - // double nxyz = GlobalC::rhopw->nx * GlobalC::rhopw->ny * GlobalC::rhopw->nz; - // dipole_elec_x *= GlobalC::ucell.omega / static_cast(nxyz); - // dipole_elec_y *= GlobalC::ucell.omega / static_cast(nxyz); - // dipole_elec_z *= GlobalC::ucell.omega / static_cast(nxyz); - - // std::cout << std::setprecision(8) << "dipole_elec_x: " << dipole_elec_x << std::endl; - // std::cout << std::setprecision(8) << "dipole_elec_y: " << dipole_elec_y << std::endl; - // std::cout << std::setprecision(8) << "dipole_elec_z: " << dipole_elec_z << std::endl; - - // ofs << std::setprecision(8)<< istep << " " << dipole_elec_x << " " << dipole_elec_y << " " << dipole_elec_z - // << std::endl; - - // /* - // double dipole_ion_x = 0.0, dipole_ion_y = 0.0, dipole_ion_z = 0.0, dipole_sum = 0.0; - // if (GlobalC::ucell.ntype == 1) - // { - // for (int ia = 0; ia < GlobalC::ucell.atoms[0].na; ia++) - // { - // dipole_ion_x += GlobalC::ucell.atoms[0].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_01; - // dipole_ion_y += GlobalC::ucell.atoms[0].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_01; - // dipole_ion_z += GlobalC::ucell.atoms[0].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_01; - // } - // } - // else if (GlobalC::ucell.ntype == 2) - // { - // for (int ia = 0; ia < GlobalC::ucell.atoms[0].na; ia++) - // { - // dipole_ion_x += GlobalC::ucell.atoms[0].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_01; - // dipole_ion_y += GlobalC::ucell.atoms[0].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_01; - // dipole_ion_z += GlobalC::ucell.atoms[0].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_01; - // } - // for (int ia = 0; ia < GlobalC::ucell.atoms[1].na; ia++) - // { - // dipole_ion_x += GlobalC::ucell.atoms[1].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_02; - // dipole_ion_y += GlobalC::ucell.atoms[1].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_02; - // dipole_ion_z += GlobalC::ucell.atoms[1].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_02; - // } - // } - // else if (GlobalC::ucell.ntype == 3) - // { - // for (int ia = 0; ia < GlobalC::ucell.atoms[0].na; ia++) - // { - // dipole_ion_x += GlobalC::ucell.atoms[0].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_01; - // dipole_ion_y += GlobalC::ucell.atoms[0].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_01; - // dipole_ion_z += GlobalC::ucell.atoms[0].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_01; - // } - // for (int ia = 0; ia < GlobalC::ucell.atoms[1].na; ia++) - // { - // dipole_ion_x += GlobalC::ucell.atoms[1].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_02; - // dipole_ion_y += GlobalC::ucell.atoms[1].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_02; - // dipole_ion_z += GlobalC::ucell.atoms[1].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_02; - // } - // for (int ia = 0; ia < GlobalC::ucell.atoms[2].na; ia++) - // { - // dipole_ion_x += GlobalC::ucell.atoms[2].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_03; - // dipole_ion_y += GlobalC::ucell.atoms[2].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_03; - // dipole_ion_z += GlobalC::ucell.atoms[2].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - // * ELEC_evolve::td_val_elec_03; - // } - // } - // else - // { - // std::cout << "atom ntype is too large!" << std::endl; - // } - // for(int it=1; it<(GlobalC::ucell.ntype); it++) - // { - // for(int ia=0; ianxyz; + dipole_elec[i] *= GlobalC::ucell.lat0 / bmod[i] * ModuleBase::FOUR_PI / GlobalC::rhopw->nxyz; } std::cout << std::setprecision(8) << "dipole_elec_x: " << dipole_elec[0] << std::endl; @@ -338,109 +92,37 @@ void ModuleIO::write_dipole(const double *rho_save, ofs << std::setprecision(8) << istep << " " << dipole_elec[0] << " " << dipole_elec[1] << " " << dipole_elec[2] << std::endl; - double dipole_ion_x = 0.0, dipole_ion_y = 0.0, dipole_ion_z = 0.0, dipole_sum = 0.0; - if (GlobalC::ucell.ntype == 1) - { - for (int ia = 0; ia < GlobalC::ucell.atoms[0].na; ia++) - { - dipole_ion_x - += GlobalC::ucell.atoms[0].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_01; - dipole_ion_y - += GlobalC::ucell.atoms[0].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_01; - dipole_ion_z - += GlobalC::ucell.atoms[0].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_01; - } - } - else if (GlobalC::ucell.ntype == 2) - { - for (int ia = 0; ia < GlobalC::ucell.atoms[0].na; ia++) - { - dipole_ion_x - += GlobalC::ucell.atoms[0].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 *INPUT.td_val_elec_01; - dipole_ion_y - += GlobalC::ucell.atoms[0].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_01; - dipole_ion_z - += GlobalC::ucell.atoms[0].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_01; - } - for (int ia = 0; ia < GlobalC::ucell.atoms[1].na; ia++) - { - dipole_ion_x - += GlobalC::ucell.atoms[1].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_02; - dipole_ion_y - += GlobalC::ucell.atoms[1].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_02; - dipole_ion_z - += GlobalC::ucell.atoms[1].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_02; - } - } - else if (GlobalC::ucell.ntype == 3) - { - for (int ia = 0; ia < GlobalC::ucell.atoms[0].na; ia++) - { - dipole_ion_x - += GlobalC::ucell.atoms[0].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_01; - dipole_ion_y - += GlobalC::ucell.atoms[0].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_01; - dipole_ion_z - += GlobalC::ucell.atoms[0].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_01; - } - for (int ia = 0; ia < GlobalC::ucell.atoms[1].na; ia++) - { - dipole_ion_x - += GlobalC::ucell.atoms[1].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_02; - dipole_ion_y - += GlobalC::ucell.atoms[1].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_02; - dipole_ion_z - += GlobalC::ucell.atoms[1].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_02; - } - for (int ia = 0; ia < GlobalC::ucell.atoms[2].na; ia++) - { - dipole_ion_x - += GlobalC::ucell.atoms[2].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_03; - dipole_ion_y - += GlobalC::ucell.atoms[2].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_03; - dipole_ion_z - += GlobalC::ucell.atoms[2].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 * INPUT.td_val_elec_03; - } - } - else - { - std::cout << "atom ntype is too large!" << std::endl; - } - for (int it = 1; it < (GlobalC::ucell.ntype); it++) - { - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) - { - dipole_ion_x += GlobalC::ucell.atoms[it].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 * 6; - dipole_ion_y += GlobalC::ucell.atoms[it].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 * 6; - dipole_ion_z += GlobalC::ucell.atoms[it].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 * 6; - } - dipole_ion_x += GlobalC::ucell.atoms[it - 1].taud[0].x * GlobalC::ucell.lat0 * 0.529177 * 1; - dipole_ion_y += GlobalC::ucell.atoms[it - 1].taud[0].y * GlobalC::ucell.lat0 * 0.529177 * 1; - dipole_ion_z += GlobalC::ucell.atoms[it - 1].taud[0].z * GlobalC::ucell.lat0 * 0.529177 * 1; - } + double dipole_ion[3] = {0.0}; + double dipole_sum = 0.0; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int i = 0; i < 3; ++i) { - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + for (int it = 0; it < GlobalC::ucell.ntype; ++it) { - dipole_ion_x += GlobalC::ucell.atoms[it].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 * 5; - dipole_ion_y += GlobalC::ucell.atoms[it].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 * 5; - dipole_ion_z += GlobalC::ucell.atoms[it].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 * 5; + 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; } + dipole_ion[i] *= GlobalC::ucell.lat0 / bmod[i] * ModuleBase::FOUR_PI / GlobalC::ucell.omega; } - std::cout << std::setprecision(8) << "dipole_ion_x: " << dipole_ion_x << std::endl; - std::cout << std::setprecision(8) << "dipole_ion_y: " << dipole_ion_y << std::endl; - std::cout << std::setprecision(8) << "dipole_ion_z: " << dipole_ion_z << std::endl; + std::cout << std::setprecision(8) << "dipole_ion_x: " << dipole_ion[0] << std::endl; + std::cout << std::setprecision(8) << "dipole_ion_y: " << dipole_ion[1] << std::endl; + std::cout << std::setprecision(8) << "dipole_ion_z: " << dipole_ion[2] << std::endl; - double dipole_x = 0.0, dipole_y = 0.0, dipole_z = 0.0; - dipole_x = dipole_ion_x - dipole_elec[0]; - dipole_y = dipole_ion_y - dipole_elec[1]; - dipole_z = dipole_ion_z - dipole_elec[2]; - std::cout << std::setprecision(8) << "dipole_x: " << dipole_x << std::endl; - std::cout << std::setprecision(8) << "dipole_y: " << dipole_y << std::endl; - std::cout << std::setprecision(8) << "dipole_z: " << dipole_z << std::endl; - dipole_sum = sqrt(dipole_x * dipole_x + dipole_y * dipole_y + dipole_z * dipole_z); + double dipole[3] = {0.0}; + for (int i = 0; i < 3; ++i) + { + dipole[i] = dipole_ion[i] - dipole_elec[i]; + } + std::cout << std::setprecision(8) << "dipole_x: " << dipole[0] << std::endl; + std::cout << std::setprecision(8) << "dipole_y: " << dipole[1] << std::endl; + std::cout << std::setprecision(8) << "dipole_z: " << dipole[2] << std::endl; + dipole_sum = sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]); std::cout << std::setprecision(8) << "dipole_sum: " << dipole_sum << std::endl; #endif @@ -454,3 +136,33 @@ void ModuleIO::write_dipole(const double *rho_save, return; } + +double ModuleIO::prepare(const UnitCell &cell, int &dir) +{ + double bvec[3] = {0.0}; + double bmod = 0.0; + if (dir == 0) + { + bvec[0] = cell.G.e11; + bvec[1] = cell.G.e12; + bvec[2] = cell.G.e13; + } + else if (dir == 1) + { + bvec[0] = cell.G.e21; + bvec[1] = cell.G.e22; + bvec[2] = cell.G.e23; + } + else if (dir == 2) + { + bvec[0] = cell.G.e31; + bvec[1] = cell.G.e32; + bvec[2] = cell.G.e33; + } + else + { + ModuleBase::WARNING_QUIT("ModuleIO::prepare", "direction is wrong!"); + } + bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2)); + return bmod; +} \ No newline at end of file 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 21b3f6cf3e..117f16dc78 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.90854261947219 -etotperatomref -15.4542713097 +etotref -30.92097334578261 +etotperatomref -15.4604866729 totalforceref 0.578804 -totalstressref 1.040362 -totaltimeref +4.6094 +totalstressref 1.159829 +totaltimeref +4.0557 From e6947a2409d4c1eccb864501a25b805b62034573 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Thu, 2 Mar 2023 18:44:01 +0800 Subject: [PATCH 2/6] fix bug of dipole and electric field --- .../potentials/H_TDDFT_pw.cpp | 73 ++++++++++++++----- .../module_elecstate/potentials/H_TDDFT_pw.h | 2 + source/module_io/write_dipole.cpp | 4 +- 3 files changed, 59 insertions(+), 20 deletions(-) diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index d817dd58b6..5e25dc4a21 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp @@ -1,10 +1,12 @@ #include "H_TDDFT_pw.h" -#include "module_io/input.h" #include "module_base/constants.h" #include "module_base/timer.h" #include "module_hamilt_lcao/module_tddft/ELEC_evolve.h" +#include "module_io/input.h" #include "module_io/input_conv.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" + namespace elecstate { @@ -16,7 +18,7 @@ int H_TDDFT_pw::istep = -1; // fuxiang add in 2017-05 //========================================================== -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"); @@ -39,18 +41,18 @@ void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo) for (auto direc: ELEC_evolve::td_vext_dire_case) { std::vector vext_space(this->rho_basis_->nrxx, 0.0); - double vext_time = cal_v_time(ttype[count]); + double vext_time = cal_v_time(ttype[count]); if (ELEC_evolve::out_efield && GlobalV::MY_RANK == 0) { std::stringstream as; - as << GlobalV::global_out_dir << "efield_"<rho_basis_->nrxx; ++ir) vl_pseudo[ir] += vext_space[ir] * vext_time; count++; @@ -94,7 +96,7 @@ void H_TDDFT_pw::read_parameters(Input *in) gauss_count = 0; gauss_omega = set_parameters(in->td_gauss_freq, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 gauss_phase = set_parameters(in->td_gauss_phase, 1.0); - gauss_sigma = set_parameters(in->td_gauss_sigma, 1/ModuleBase::AU_to_FS); + gauss_sigma = set_parameters(in->td_gauss_sigma, 1 / ModuleBase::AU_to_FS); gauss_t0 = set_parameters(in->td_gauss_t0, 1.0); gauss_amp = set_parameters(in->td_gauss_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr @@ -163,6 +165,12 @@ 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"); + double bmod[3]; + for (int i = 0; i < 3; i++) + { + bmod[i] = prepare(GlobalC::ucell, i); + } + for (int ir = 0; ir < this->rho_basis_->nrxx; ++ir) { int i = ir / (this->rho_basis_->ny * this->rho_basis_->nplane); @@ -175,15 +183,15 @@ 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); + vext_space[ir] = cal_v_space_length_potential(x) / bmod[0]; break; case 2: - vext_space[ir] = cal_v_space_length_potential(y); + vext_space[ir] = cal_v_space_length_potential(y) / bmod[1]; break; case 3: - vext_space[ir] = cal_v_space_length_potential(z); + vext_space[ir] = cal_v_space_length_potential(z) / bmod[2]; break; default: @@ -198,10 +206,10 @@ void H_TDDFT_pw::cal_v_space_length(std::vector &vext_space, int direc) double H_TDDFT_pw::cal_v_space_length_potential(double i) { - double vext_space=0.0; + 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) { @@ -209,7 +217,7 @@ double H_TDDFT_pw::cal_v_space_length_potential(double i) } 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; } @@ -241,9 +249,9 @@ double H_TDDFT_pw::cal_v_time(int t_type) vext_time = cal_v_time_heaviside(); break; - // case 4: - // vext_time = cal_v_time_HHG(); - // break; + // case 4: + // vext_time = cal_v_time_HHG(); + // break; default: std::cout << "time_domain_type of electric field is wrong" << endl; @@ -308,8 +316,7 @@ double H_TDDFT_pw::cal_v_time_trigonometric() const double timenow = istep * dt; - vext_time = amp * cos(omega1 * timenow + phase1) * sin(omega2 * timenow + phase2) - * sin(omega2 * timenow + phase2); + vext_time = amp * cos(omega1 * timenow + phase1) * sin(omega2 * timenow + phase2) * sin(omega2 * timenow + phase2); trigo_count++; return vext_time; @@ -349,4 +356,34 @@ double H_TDDFT_pw::cal_v_time_heaviside() // return vext_time; // } +double H_TDDFT_pw::prepare(const UnitCell &cell, int &dir) +{ + double bvec[3] = {0.0}; + double bmod = 0.0; + if (dir == 0) + { + bvec[0] = cell.G.e11; + bvec[1] = cell.G.e12; + bvec[2] = cell.G.e13; + } + else if (dir == 1) + { + bvec[0] = cell.G.e21; + bvec[1] = cell.G.e22; + bvec[2] = cell.G.e23; + } + else if (dir == 2) + { + bvec[0] = cell.G.e31; + bvec[1] = cell.G.e32; + bvec[2] = cell.G.e33; + } + else + { + ModuleBase::WARNING_QUIT("H_TDDFT_pw::prepare", "direction is wrong!"); + } + bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2)); + return bmod; +} + } // namespace elecstate diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.h b/source/module_elecstate/potentials/H_TDDFT_pw.h index df9ea2b3d3..e600219b24 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.h +++ b/source/module_elecstate/potentials/H_TDDFT_pw.h @@ -109,6 +109,8 @@ class H_TDDFT_pw : public PotBase double cal_v_time_trigonometric(); double cal_v_time_heaviside(); // double cal_v_time_HHG(); + + double prepare(const UnitCell &cell, int &dir); }; } // namespace elecstate diff --git a/source/module_io/write_dipole.cpp b/source/module_io/write_dipole.cpp index 383268512d..1b9dc884af 100644 --- a/source/module_io/write_dipole.cpp +++ b/source/module_io/write_dipole.cpp @@ -82,7 +82,7 @@ void ModuleIO::write_dipole(const double *rho_save, Parallel_Reduce::reduce_double_pool(dipole_elec[2]); for (int i = 0; i < 3; ++i) { - dipole_elec[i] *= GlobalC::ucell.lat0 / bmod[i] * ModuleBase::FOUR_PI / GlobalC::rhopw->nxyz; + dipole_elec[i] *= GlobalC::ucell.lat0 / bmod[i] * GlobalC::ucell.omega / GlobalC::rhopw->nxyz; } std::cout << std::setprecision(8) << "dipole_elec_x: " << dipole_elec[0] << std::endl; @@ -107,7 +107,7 @@ void ModuleIO::write_dipole(const double *rho_save, } dipole_ion[i] += sum * GlobalC::ucell.atoms[it].ncpp.zv; } - dipole_ion[i] *= GlobalC::ucell.lat0 / bmod[i] * ModuleBase::FOUR_PI / GlobalC::ucell.omega; + dipole_ion[i] *= GlobalC::ucell.lat0 / bmod[i]; //* ModuleBase::FOUR_PI / GlobalC::ucell.omega; } std::cout << std::setprecision(8) << "dipole_ion_x: " << dipole_ion[0] << std::endl; From 0b6b618d85fbfd21eb97c74bb19c1698c50dde66 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Tue, 7 Mar 2023 22:30:25 +0800 Subject: [PATCH 3/6] fix bug of multi-kpoint for tddft --- source/module_esolver/esolver_ks_lcao_tddft.cpp | 6 ++++-- source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp | 1 + tests/integrate/601_NO_TDDFT_CO/result.ref | 8 ++++---- tests/integrate/601_NO_TDDFT_CO_occ/result.ref | 8 ++++---- tests/integrate/601_NO_TDDFT_H2/result.ref | 4 ++-- tests/integrate/601_NO_TDDFT_H2_kpoint/result.ref | 10 +++++----- tests/integrate/601_NO_TDDFT_H2_len_gauss/result.ref | 10 +++++----- .../601_NO_TDDFT_H2_len_gauss_dire/result.ref | 10 +++++----- tests/integrate/601_NO_TDDFT_H2_len_heavi/result.ref | 10 +++++----- tests/integrate/601_NO_TDDFT_H2_len_hhg/result.ref | 10 +++++----- tests/integrate/601_NO_TDDFT_H2_len_trape/result.ref | 10 +++++----- tests/integrate/601_NO_TDDFT_H2_len_trigo/result.ref | 10 +++++----- tests/integrate/601_NO_TDDFT_H2_oldedm/result.ref | 10 +++++----- tests/integrate/601_NO_TDDFT_O3/result.ref | 6 +++--- .../integrate/601_NO_TDDFT_graphene_kpoint/result.ref | 10 +++++----- 15 files changed, 63 insertions(+), 60 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 5533622d3e..ee0af9aeb3 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -306,8 +306,10 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) #endif for (int ik = 0; ik < GlobalC::kv.nks; ++ik) { - psi->fix_k(ik); - for (int index = 0; index < psi[0].size(); ++index) + this->psi->fix_k(ik); + this->psi_laststep->fix_k(ik); + int size0 = psi->get_nbands() * psi->get_nbasis(); + for (int index = 0; index < size0; ++index) psi_laststep[0].get_pointer()[index] = psi[0].get_pointer()[index]; } if (istep > 1 && ELEC_evolve::td_edm == 0) diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp index 39a4a18e69..ade328fce4 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp +++ b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp @@ -41,6 +41,7 @@ void ELEC_evolve::evolve_psi(const int& istep, ModuleBase::timer::tick("Efficience", "evolve_k"); Evolve_LCAO_Matrix ELM(lowf.ParaV); psi->fix_k(ik); + psi_laststep->fix_k(ik); ELM.evolve_complex_matrix(ik, phm, psi, psi_laststep, &(ekb(ik, 0))); ModuleBase::timer::tick("Efficience", "evolve_k"); } // end k diff --git a/tests/integrate/601_NO_TDDFT_CO/result.ref b/tests/integrate/601_NO_TDDFT_CO/result.ref index e7281e1e71..27a1c2fe89 100644 --- a/tests/integrate/601_NO_TDDFT_CO/result.ref +++ b/tests/integrate/601_NO_TDDFT_CO/result.ref @@ -1,5 +1,5 @@ -etotref -602.8643309778789 -etotperatomref -301.4321654889 +etotref -602.8643309993128 +etotperatomref -301.4321654997 totalforceref 16.350762 -totalstressref 30.184707 -totaltimeref +8.5807 +totalstressref 30.184735 +totaltimeref +7.1335 diff --git a/tests/integrate/601_NO_TDDFT_CO_occ/result.ref b/tests/integrate/601_NO_TDDFT_CO_occ/result.ref index 2b3b7e9030..271e8953a8 100644 --- a/tests/integrate/601_NO_TDDFT_CO_occ/result.ref +++ b/tests/integrate/601_NO_TDDFT_CO_occ/result.ref @@ -1,5 +1,5 @@ -etotref -602.8643309786496 -etotperatomref -301.4321654893 +etotref -602.8643311305311 +etotperatomref -301.4321655653 totalforceref 16.350762 -totalstressref 30.184729 -totaltimeref +8.5178 +totalstressref 30.185015 +totaltimeref +7.0590 diff --git a/tests/integrate/601_NO_TDDFT_H2/result.ref b/tests/integrate/601_NO_TDDFT_H2/result.ref index 0fcd92098e..9bffc30214 100644 --- a/tests/integrate/601_NO_TDDFT_H2/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2/result.ref @@ -1,5 +1,5 @@ -etotref -18.05421597489365 +etotref -18.05421597483996 etotperatomref -9.0271079874 totalforceref 44.953238 totalstressref 79.612084 -totaltimeref +8.0145 +totaltimeref +3.7961 diff --git a/tests/integrate/601_NO_TDDFT_H2_kpoint/result.ref b/tests/integrate/601_NO_TDDFT_H2_kpoint/result.ref index a4417260f7..40dca31f18 100644 --- a/tests/integrate/601_NO_TDDFT_H2_kpoint/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_kpoint/result.ref @@ -1,5 +1,5 @@ -etotref -18.05387244304529 -etotperatomref -9.0269362215 -totalforceref 45.003706 -totalstressref 79.777732 -totaltimeref +6.7899 +etotref -18.05362763740304 +etotperatomref -9.0268138187 +totalforceref 44.990644 +totalstressref 79.754356 +totaltimeref +3.8596 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 9f9f5c8814..f7f4110f52 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.36567840538823 -etotperatomref -15.6828392027 -totalforceref 0.578716 -totalstressref 5.434008 -totaltimeref +5.1729 +etotref -31.36567798859614 +etotperatomref -15.6828389943 +totalforceref 0.578604 +totalstressref 5.433996 +totaltimeref +3.9632 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 5c53d29af7..59c1620fa8 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.36554134768115 -etotperatomref -15.6827706738 -totalforceref 0.577344 -totalstressref 5.434175 -totaltimeref +5.1576 +etotref -31.36554049797190 +etotperatomref -15.6827702490 +totalforceref 0.577232 +totalstressref 5.434161 +totaltimeref +4.0295 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 65dc2a84cc..3e52430584 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.85116741796947 -etotperatomref -19.4255837090 -totalforceref 0.568680 -totalstressref 77.448522 -totaltimeref +5.0160 +etotref -38.85116666973314 +etotperatomref -19.4255833349 +totalforceref 0.568498 +totalstressref 77.448537 +totaltimeref +3.8784 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 4a5f661121..b6d67e9d5b 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.60900810491984 -etotperatomref -41.3045040525 -totalforceref 0.757204 -totalstressref 500.961261 -totaltimeref +5.5047 +etotref -82.60904651293552 +etotperatomref -41.3045232565 +totalforceref 0.754930 +totalstressref 500.961654 +totaltimeref +4.0425 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 117f16dc78..4f21088bfc 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.92097334578261 -etotperatomref -15.4604866729 -totalforceref 0.578804 -totalstressref 1.159829 -totaltimeref +4.0557 +etotref -30.92097257037962 +etotperatomref -15.4604862852 +totalforceref 0.578690 +totalstressref 1.159819 +totaltimeref +3.5993 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 39130ba6ee..542503a9f7 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.91150694066410 -etotperatomref -15.4557534703 -totalforceref 0.578804 -totalstressref 1.068851 -totaltimeref +4.6392 +etotref -30.91150616526584 +etotperatomref -15.4557530826 +totalforceref 0.578692 +totalstressref 1.068838 +totaltimeref +3.6067 diff --git a/tests/integrate/601_NO_TDDFT_H2_oldedm/result.ref b/tests/integrate/601_NO_TDDFT_H2_oldedm/result.ref index c03640d32c..ff03624d80 100644 --- a/tests/integrate/601_NO_TDDFT_H2_oldedm/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_oldedm/result.ref @@ -1,5 +1,5 @@ -etotref -18.05387244303919 -etotperatomref -9.0269362215 -totalforceref 45.003706 -totalstressref 79.777732 -totaltimeref +5.0273 +etotref -18.05362763740304 +etotperatomref -9.0268138187 +totalforceref 44.990644 +totalstressref 79.754356 +totaltimeref +3.8057 diff --git a/tests/integrate/601_NO_TDDFT_O3/result.ref b/tests/integrate/601_NO_TDDFT_O3/result.ref index e1a6292780..d9794538f9 100644 --- a/tests/integrate/601_NO_TDDFT_O3/result.ref +++ b/tests/integrate/601_NO_TDDFT_O3/result.ref @@ -1,5 +1,5 @@ -etotref -1335.872362262912 -etotperatomref -445.2907874210 +etotref -1335.872362260576 +etotperatomref -445.2907874202 totalforceref 13.319759 totalstressref 54.934057 -totaltimeref +27.907 +totaltimeref +23.556 diff --git a/tests/integrate/601_NO_TDDFT_graphene_kpoint/result.ref b/tests/integrate/601_NO_TDDFT_graphene_kpoint/result.ref index d7581b6f31..2314206da8 100644 --- a/tests/integrate/601_NO_TDDFT_graphene_kpoint/result.ref +++ b/tests/integrate/601_NO_TDDFT_graphene_kpoint/result.ref @@ -1,5 +1,5 @@ -etotref -257.4519769034942 -etotperatomref -85.8173256345 -totalforceref 10.905441 -totalstressref 2485.841436 -totaltimeref +84.603 +etotref -321.2175762983792 +etotperatomref -107.0725254328 +totalforceref 13.778433 +totalstressref 1806.432527 +totaltimeref +35.306 From ed13f3681077e2c11eead65a1958b78831b5d462 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Tue, 7 Mar 2023 22:53:09 +0800 Subject: [PATCH 4/6] unit conversion of out_efield --- source/module_elecstate/potentials/H_TDDFT_pw.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index 5e25dc4a21..a915af5b7b 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp @@ -3,10 +3,9 @@ #include "module_base/constants.h" #include "module_base/timer.h" #include "module_hamilt_lcao/module_tddft/ELEC_evolve.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/input.h" #include "module_io/input_conv.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" - namespace elecstate { @@ -48,7 +47,8 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) std::stringstream as; as << GlobalV::global_out_dir << "efield_" << count << ".dat"; std::ofstream ofs(as.str().c_str(), std::ofstream::app); - ofs << H_TDDFT_pw::istep * dt * ModuleBase::AU_to_FS << "\t" << vext_time << endl; + ofs << H_TDDFT_pw::istep * dt * ModuleBase::AU_to_FS << "\t" + << vext_time * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << endl; ofs.close(); } From 8047e37f36e4edaa870bcf99eea56a99e3822c4b Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Tue, 7 Mar 2023 22:56:29 +0800 Subject: [PATCH 5/6] delete useless paramters for tddft (td_val_elec*) --- .../module_tddft/ELEC_evolve.cpp | 3 --- .../module_tddft/ELEC_evolve.h | 3 --- source/module_io/input.cpp | 18 ------------------ source/module_io/input.h | 3 --- source/module_io/input_conv.cpp | 3 --- source/module_io/test/INPUT | 3 --- source/module_io/test/input_test.cpp | 6 ------ source/module_io/write_input.cpp | 3 --- 8 files changed, 42 deletions(-) diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp index ade328fce4..66a8fcb1a7 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp +++ b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp @@ -12,9 +12,6 @@ ELEC_evolve::ELEC_evolve(){}; ELEC_evolve::~ELEC_evolve(){}; double ELEC_evolve::td_force_dt; -int ELEC_evolve::td_val_elec_01; -int ELEC_evolve::td_val_elec_02; -int ELEC_evolve::td_val_elec_03; bool ELEC_evolve::td_vext; std::vector ELEC_evolve::td_vext_dire_case; bool ELEC_evolve::out_dipole; diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h index 86277853c6..6571e72033 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h +++ b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h @@ -30,9 +30,6 @@ class ELEC_evolve // fuxiang add 2021-05-25 static double td_force_dt; - static int td_val_elec_01; - static int td_val_elec_02; - static int td_val_elec_03; static bool td_vext; static std::vector td_vext_dire_case; static bool out_dipole; diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 282c9545a2..5f0bf24f37 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -396,9 +396,6 @@ void Input::Default(void) // tddft //---------------------------------------------------------- td_force_dt = 0.02; - td_val_elec_01 = 1; - td_val_elec_02 = 1; - td_val_elec_03 = 1; td_vext = false; td_vext_dire = "1"; @@ -1483,18 +1480,6 @@ bool Input::Read(const std::string &fn) { read_value(ifs, td_force_dt); } - else if (strcmp("td_val_elec_01", word) == 0) - { - read_value(ifs, td_val_elec_01); - } - else if (strcmp("td_val_elec_02", word) == 0) - { - read_value(ifs, td_val_elec_02); - } - else if (strcmp("td_val_elec_03", word) == 0) - { - read_value(ifs, td_val_elec_03); - } else if (strcmp("td_vext", word) == 0) { read_value(ifs, td_vext); @@ -2909,9 +2894,6 @@ void Input::Bcast() Parallel_Common::bcast_int(vdw_cutoff_period.y); Parallel_Common::bcast_int(vdw_cutoff_period.z); // Fuxiang He add 2016-10-26 - Parallel_Common::bcast_int(td_val_elec_01); - Parallel_Common::bcast_int(td_val_elec_02); - Parallel_Common::bcast_int(td_val_elec_03); Parallel_Common::bcast_double(td_force_dt); Parallel_Common::bcast_bool(td_vext); Parallel_Common::bcast_string(td_vext_dire); diff --git a/source/module_io/input.h b/source/module_io/input.h index 250e4e7ba1..0473bc8a3d 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -381,9 +381,6 @@ class Input // Fuxiang He add 2016-10-26 //========================================================== double td_force_dt; //"fs" - int td_val_elec_01; // valence electron 01 - int td_val_elec_02; // valence electron 02 - int td_val_elec_03; // valence electron 03 bool td_vext; // add extern potential or not std::string td_vext_dire; // vext direction bool out_dipole; // output the dipole or not diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 88a14837c4..996a77c0a6 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -355,9 +355,6 @@ void Input_Conv::Convert(void) //---------------------------------------------------------- #ifdef __LCAO ELEC_evolve::td_force_dt = INPUT.td_force_dt; - ELEC_evolve::td_val_elec_01 = INPUT.td_val_elec_01; - ELEC_evolve::td_val_elec_02 = INPUT.td_val_elec_02; - ELEC_evolve::td_val_elec_03 = INPUT.td_val_elec_03; ELEC_evolve::td_vext = INPUT.td_vext; if (ELEC_evolve::td_vext) { diff --git a/source/module_io/test/INPUT b/source/module_io/test/INPUT index 61a300ed53..8a3f3f0e6b 100644 --- a/source/module_io/test/INPUT +++ b/source/module_io/test/INPUT @@ -273,9 +273,6 @@ exx_opt_orb_tolerence 0 # #Parameters (16.tddft) td_force_dt 0.02 #time of force change -td_val_elec_01 1 #td_val_elec_01 -td_val_elec_02 1 #td_val_elec_02 -td_val_elec_03 1 #td_val_elec_03 td_vext 0 #add extern potential or not td_vext_dire 1 #extern potential direction td_stype 0 #space domain type diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index 469be31773..b3c7ca5a58 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -244,9 +244,6 @@ TEST_F(InputTest, Default) EXPECT_DOUBLE_EQ(INPUT.soc_lambda,1.0); EXPECT_EQ(INPUT.input_error,0); EXPECT_DOUBLE_EQ(INPUT.td_force_dt,0.02); - EXPECT_EQ(INPUT.td_val_elec_01,1); - EXPECT_EQ(INPUT.td_val_elec_02,1); - EXPECT_EQ(INPUT.td_val_elec_03,1); EXPECT_FALSE(INPUT.td_vext); EXPECT_EQ(INPUT.td_vext_dire,"1"); EXPECT_EQ(INPUT.td_stype,0); @@ -571,9 +568,6 @@ TEST_F(InputTest, Read) EXPECT_DOUBLE_EQ(INPUT.soc_lambda,1.0); EXPECT_EQ(INPUT.input_error,0); EXPECT_DOUBLE_EQ(INPUT.td_force_dt,0.02); - EXPECT_EQ(INPUT.td_val_elec_01,1); - EXPECT_EQ(INPUT.td_val_elec_02,1); - EXPECT_EQ(INPUT.td_val_elec_03,1); EXPECT_EQ(INPUT.td_vext,0); // EXPECT_EQ(INPUT.td_vext_dire,"1"); EXPECT_EQ(INPUT.td_stype,0); diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index c1ab4ca37c..ccdae2f987 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -353,9 +353,6 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou ofs << "\n#Parameters (16.tddft)" << std::endl; ModuleBase::GlobalFunc::OUTP(ofs, "td_force_dt", td_force_dt, "time of force change"); - ModuleBase::GlobalFunc::OUTP(ofs, "td_val_elec_01", td_val_elec_01, "td_val_elec_01"); - ModuleBase::GlobalFunc::OUTP(ofs, "td_val_elec_02", td_val_elec_02, "td_val_elec_02"); - ModuleBase::GlobalFunc::OUTP(ofs, "td_val_elec_03", td_val_elec_03, "td_val_elec_03"); ModuleBase::GlobalFunc::OUTP(ofs, "td_vext", td_vext, "add extern potential or not"); ModuleBase::GlobalFunc::OUTP(ofs, "td_vext_dire", td_vext_dire, "extern potential direction"); ModuleBase::GlobalFunc::OUTP(ofs, "out_dipole", out_dipole, "output dipole or not"); From f834b83e79cb553ccf1e182bb4985bc60487855e Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Wed, 8 Mar 2023 13:49:08 +0800 Subject: [PATCH 6/6] delete useless parameters for tddft --- docs/advanced/input_files/input-main.md | 20 +------------------- 1 file changed, 1 insertion(+), 19 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index d3a7efa913..bb6e81ab88 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -41,7 +41,7 @@ [berry_phase](#berry_phase) | [gdir](#gdir) | [towannier90](#towannier90) | [nnkpfile](#nnkpfile) | [wannier_spin](#wannier_spin) - [TDDFT: time dependent density functional theory](#tddft-time-dependent-density-functional-theory) (Under tests) - [td_force_dt](#td_force_dt) | [td_vext](#td_vext) | [td_vext_dire](#td_vext_dire) | [td_stype](#td_stype) | [td_ttype](#td_ttype) | [td_tstart](#td_tstart) | [td_tend](#td_tend) | [td_lcut1](#td_lcut1) | [td_lcut2](#td_lcut2) | [td_gauss_freq](#td_gauss_freq) | [td_guass_phase](#td_gauss_phase) | [td_gauss_sigma](#td_gauss_sigma) | [td_gauss_t0](#td_gauss_t0)| [td_gauss_amp](#td_gauss_amp) | [td_trape_freq](#td_trape_freq) | [td_trape_phase](#td_trape_phase) | [td_trape_t1](#td_trape_t1) | [td_trape_t2](#td_trape_t2) | [td_trape_t3](#td_trape_t3) | [td_trape_amp](#td_trape_amp) | [td_trigo_freq1](#td_trigo_freq1) | [td_trigo_freq2](#td_trigo_freq2) | [td_trigo_phase1](#td_trigo_phase1) | [td_trigo_phase2](#td_trigo_phase2) | [td_trigo_amp](#td_trigo_amp) | [td_heavi_t0](#td_heavi_t0) | [td_heavi_amp](#td_heavi_amp) | [out_dipole](#out_dipole) |[out_efield](#out_efield)| [ocp](#ocp) | [ocp_set](#ocp_set) | [td_val_elec_01](#td_val_elec_01) | [td_val_elec_02](#td_val_elec_02) |[td_val_elec_03](#td_val_elec_03) + [td_force_dt](#td_force_dt) | [td_vext](#td_vext) | [td_vext_dire](#td_vext_dire) | [td_stype](#td_stype) | [td_ttype](#td_ttype) | [td_tstart](#td_tstart) | [td_tend](#td_tend) | [td_lcut1](#td_lcut1) | [td_lcut2](#td_lcut2) | [td_gauss_freq](#td_gauss_freq) | [td_guass_phase](#td_gauss_phase) | [td_gauss_sigma](#td_gauss_sigma) | [td_gauss_t0](#td_gauss_t0)| [td_gauss_amp](#td_gauss_amp) | [td_trape_freq](#td_trape_freq) | [td_trape_phase](#td_trape_phase) | [td_trape_t1](#td_trape_t1) | [td_trape_t2](#td_trape_t2) | [td_trape_t3](#td_trape_t3) | [td_trape_amp](#td_trape_amp) | [td_trigo_freq1](#td_trigo_freq1) | [td_trigo_freq2](#td_trigo_freq2) | [td_trigo_phase1](#td_trigo_phase1) | [td_trigo_phase2](#td_trigo_phase2) | [td_trigo_amp](#td_trigo_amp) | [td_heavi_t0](#td_heavi_t0) | [td_heavi_amp](#td_heavi_amp) | [out_dipole](#out_dipole) |[out_efield](#out_efield)| [ocp](#ocp) | [ocp_set](#ocp_set) | - [DFT+*U* correction](#dftu-correction) (Under development) [dft_plus_u](#dft_plus_u) | [orbital_corr](#orbital_corr) | [hubbard_u](#hubbard_u) | [yukawa_potential](#yukawa_potential) | [yukawa_lambda](#yukawa_lambda) | [omc](#omc) @@ -2183,24 +2183,6 @@ These variables are used to control berry phase and wannier90 interface paramete - **Description**: If ocp is true, the ocp_set is a string to set the number of occupancy, like 1 10 * 1 0 1 representing the 13 band occupancy, 12th band occupancy 0 and the rest 1, the code is parsing this string into an array through a regular expression. - **Default**: none -### td_val_elec_01 - -- **Type**: Integer -- **Description**: Only useful when calculating the dipole. Specifies the number of valence electron associated with the first element. -- **Default**: 1.0 - -### td_val_elec_02 - -- **Type**: Integer -- **Description**: Only useful when calculating the dipole. Specifies the number of valence electron associated with the second element. -- **Default**: 1.0 - -### td_val_elec_03 - -- **Type**: Integer -- **Description**: Only useful when calculating the dipole. Specifies the number of valence electron associated with the third element. -- **Default**: 1.0 - [back to top](#full-list-of-input-keywords) ## Variables useful for debugging