From 2d3d30e2e17d2bf04254f5fd4df052ff8bbf4eb4 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Sun, 25 Jan 2026 17:27:12 +0800 Subject: [PATCH] Refactor&Doc: simplify read_atoms.cpp and update stru.md doc --- docs/advanced/input_files/stru.md | 39 +- docs/quick_start/input.md | 2 +- source/Makefile.Objects | 1 + source/source_cell/CMakeLists.txt | 1 + source/source_cell/read_atoms.cpp | 530 +--------------- source/source_cell/read_atoms_helper.cpp | 596 ++++++++++++++++++ source/source_cell/read_atoms_helper.h | 136 ++++ source/source_cell/test/CMakeLists.txt | 26 + .../test/read_atoms_helper_test.cpp | 483 ++++++++++++++ source/source_cell/test_pw/CMakeLists.txt | 2 +- .../module_deepks/test/CMakeLists.txt | 1 + .../module_deepks/test/Makefile.Objects | 1 + source/source_md/test/CMakeLists.txt | 1 + .../module_pwdft/test/CMakeLists.txt | 1 + 14 files changed, 1319 insertions(+), 501 deletions(-) create mode 100644 source/source_cell/read_atoms_helper.cpp create mode 100644 source/source_cell/read_atoms_helper.h create mode 100644 source/source_cell/test/read_atoms_helper_test.cpp diff --git a/docs/advanced/input_files/stru.md b/docs/advanced/input_files/stru.md index dff5ae3020..a30064351a 100644 --- a/docs/advanced/input_files/stru.md +++ b/docs/advanced/input_files/stru.md @@ -99,6 +99,7 @@ information that comes below. 2. [SG15-ONCV](http://quantum-simulation.org/potentials/sg15_oncv/upf/). 3. [DOJO](http://www.pseudo-dojo.org/). 4. [BLPS](https://github.com/PrincetonUniversity/BLPSLibrary). + 5. For additional pseudopotential options and to view the basic benchmark test results of these pseudopotentials in ABACUS, please refer to the [Benchmarks website](https://kirk0830.github.io/ABACUS-Pseudopot-Nao-Square/pseudopotential/pseudopotential.html) ### NUMERICAL_ORBITAL @@ -110,6 +111,9 @@ information that comes below. ‘Si_gga_8au_60Ry_2s2p1d.orb’ is name of the numerical orbital file. Again here the path is not specified, which means that this file is located in the work directory. Numerical atomic orbitals may be downloaded from the [official website](http://abacus.ustc.edu.cn/pseudo/list.htm). + Recommendation for Pseudopotential and Orbital Sets +For general usage requirements, the APNSv1.0 pseudopotential and orbital set is recommended. You can access it via [AIS square website](https://www.aissquare.com/datasets/detail?pageType=datasets&name=ABACUS-APNS-PPORBs-v1%253Apre-release&id=326) + ### LATTICE_CONSTANT The lattice constant of the system in unit of Bohr. @@ -248,7 +252,7 @@ information that comes below. Several other parameters could be defined after the atom position using key words : - `m` or NO key word: three numbers, which take value in 0 or 1, control how the atom move in geometry relaxation calculations. In example below, the numbers `0 0 0` following the coordinates of the first atom means this atom are *not allowed* to move in all three directions, and the numbers `1 1 1` following the coordinates of the second atom means this atom *can* move in all three directions. - - `v` or `vel` or `velocity`: set the three components of initial velocity of atoms in geometry relaxation calculations(e. g. `v 1.0 1.0 1.0`). + - `v` or `vel` or `velocity`: set the three components of initial velocity of atoms, used only for restarting MD calculations (e.g., `v 1.0 1.0 1.0`). - `mag` or `magmom` : set the start magnetization for each atom. In colinear case only one number should be given. In non-colinear case one have two choice:either set one number for the norm of magnetization here and specify two polar angle later(e. g. see below), or set three number for the xyz commponent of magnetization here (e. g. `mag 0.0 0.0 1.0`). Note that if this parameter is set, the initial magnetic moment setting in the second line will be overrided. - `angle1`: in non-colinear case, specify the angle between z-axis and real spin, in angle measure instead of radian measure - `angle2`: in non-colinear case, specify angle between x-axis and real spin in projection in xy-plane , in angle measure instead of radian measure @@ -305,4 +309,35 @@ information that comes below. 0.0 0.0 0.0 m 0 0 0 mag 1 1 1 0.5 0.5 0.5 m 1 1 1 mag 1 1 1 ``` - However, this autoset will not be vaild once `STRU` specalize a finite magnetic for any single atom. \ No newline at end of file + However, this autoset will not be vaild once `STRU` specalize a finite magnetic for any single atom. + + - `lambda`: Lagrange multiplier vector for spin constraint method. Can specify one value (z-component) or three values for x, y, z components (e.g., `lambda 0.5` or `lambda 0.1 0.2 0.3`). Values are in eV and will be converted to Rydberg internally. Used with spin-constrained DFT (enable with `sc_mag_switch` in INPUT file). + + - `sc`: set the spin constraint target magnetization for each atom. Can specify one value (z-component) or three values for x, y, z components (e.g., `sc 1.0` or `sc 0.5 0.5 1.0`). Used with spin-constrained DFT (enable with `sc_mag_switch` in INPUT file). + +### Important Notes for ATOMIC_POSITIONS + +1. **Coordinate System Selection**: Choose the appropriate coordinate system based on your needs: + - Use `Direct` for fractional coordinates (most common for periodic systems) + - Use `Cartesian_angstrom` when working with molecular structures or experimental data + - Use centered coordinate systems (`Cartesian_angstrom_center_xy/xz/yz/xyz`) for surface or slab calculations where you want to center the structure + +2. **Magnetization Settings**: + - For collinear calculations (`nspin=2`), only specify one magnetization value per atom + - For non-collinear calculations (`nspin=4`), you can specify: + - Three components directly: `mag 1.0 0.0 0.0` (mx, my, mz) + - Magnitude with angles: `mag 1.0 angle1 90 angle2 0` (magnitude, polar angle, azimuthal angle) + - If no magnetization is specified for any atom, ABACUS will automatically set default values (1.0 for nspin=2, or (1,1,1) for nspin=4) + +3. **Movement Constraints**: + - Use `m 1 1 1` to allow the atom to move freely in all directions during relaxation + - Use `m 0 0 0` to fix the atom completely + - Use `m 1 0 1` to allow movement only in x and z directions (useful for constraining surface atoms) + +4. **Keyword Order**: The optional keywords (m, v, mag, angle1, angle2, lambda, sc) can appear in any order after the atomic coordinates, but each keyword should only appear once per atom. + +5. **Common Mistakes to Avoid**: + - Don't mix Direct and Cartesian coordinates in the same STRU file + - Ensure the number of atoms specified matches the actual number of coordinate lines provided + - When using vector magnetization (`mag x y z`), don't also specify angles for the same atom + - Remember that angles are in degrees, not radians \ No newline at end of file diff --git a/docs/quick_start/input.md b/docs/quick_start/input.md index c13d04e342..d7a94449e0 100644 --- a/docs/quick_start/input.md +++ b/docs/quick_start/input.md @@ -95,7 +95,7 @@ O #Name of element > **Note:** users may choose a different name for their structure file using the keyword `stru_file`. The order of the pseudopotential file list and the numerical orbital list (if LCAO is applied) MUST be consistent with that of the atomic type given in `ATOMIC_POSITIONS`. -For a more detailed description of STRU file, please consult [here](../advanced/input_files/stru.md). +> **Important:** When specifying atomic positions, you can use various coordinate systems (Direct, Cartesian, Cartesian_angstrom, etc.) and add optional properties like magnetization, velocity, and movement constraints. See the [detailed STRU documentation](../advanced/input_files/stru.md) for all available options and best practices. ## *KPT* diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 4b8288b5f5..4269b03c99 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -191,6 +191,7 @@ OBJS_CELL=atom_pseudo.o\ read_pp_vwr.o\ unitcell.o\ read_atoms.o\ + read_atoms_helper.o\ print_cell.o\ setup_nonlocal.o\ klist.o\ diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index f9d4f0c535..39cfa88c22 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -15,6 +15,7 @@ add_library( read_pp_vwr.cpp unitcell.cpp read_atoms.cpp + read_atoms_helper.cpp setup_nonlocal.cpp klist.cpp parallel_kpoints.cpp diff --git a/source/source_cell/read_atoms.cpp b/source/source_cell/read_atoms.cpp index cdce715ce7..82be874db5 100644 --- a/source/source_cell/read_atoms.cpp +++ b/source/source_cell/read_atoms.cpp @@ -4,9 +4,10 @@ #include #include "unitcell.h" +#include "read_atoms_helper.h" #include "source_io/module_parameter/parameter.h" -#include "source_cell/print_cell.h" -#include "source_cell/read_stru.h" +#include "print_cell.h" +#include "read_stru.h" #include "source_estate/read_orb.h" #include "source_base/timer.h" #include "source_base/constants.h" @@ -19,46 +20,26 @@ #endif bool unitcell::read_atom_positions(UnitCell& ucell, - std::ifstream &ifpos, - std::ofstream &ofs_running, + std::ifstream &ifpos, + std::ofstream &ofs_running, std::ofstream &ofs_warning) { ModuleBase::TITLE("UnitCell","read_atom_positions"); std::string& Coordinate = ucell.Coordinate; const int ntype = ucell.ntype; - const int nspin = PARAM.inp.nspin; + const int nspin = PARAM.inp.nspin; assert (nspin==1 || nspin==2 || nspin==4); if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifpos, "ATOMIC_POSITIONS")) { ModuleBase::GlobalFunc::READ_VALUE(ifpos, Coordinate); - if(Coordinate != "Cartesian" - && Coordinate != "Direct" - && Coordinate != "Cartesian_angstrom" - && Coordinate != "Cartesian_au" - && Coordinate != "Cartesian_angstrom_center_xy" - && Coordinate != "Cartesian_angstrom_center_xz" - && Coordinate != "Cartesian_angstrom_center_yz" - && Coordinate != "Cartesian_angstrom_center_xyz" - ) + if (!unitcell::validate_coordinate_system(Coordinate, ofs_warning)) { - ModuleBase::WARNING("read_atom_position","Cartesian or Direct?"); - ofs_warning << " There are several options for you:" << std::endl; - ofs_warning << " Direct" << std::endl; - ofs_warning << " Cartesian_angstrom" << std::endl; - ofs_warning << " Cartesian_au" << std::endl; - ofs_warning << " Cartesian_angstrom_center_xy" << std::endl; - ofs_warning << " Cartesian_angstrom_center_xz" << std::endl; - ofs_warning << " Cartesian_angstrom_center_yz" << std::endl; - ofs_warning << " Cartesian_angstrom_center_xyz" << std::endl; - return false; // means something wrong + return false; } - ModuleBase::Vector3 v; - ModuleBase::Vector3 mv; - int na = 0; ucell.nat = 0; //====================================== @@ -68,124 +49,25 @@ bool unitcell::read_atom_positions(UnitCell& ucell, for (int it = 0;it < ntype; it++) { ofs_running << "\n READING ATOM TYPE " << it+1 << std::endl; - - //======================================= - // (1) read in atom label - // start magnetization - //======================================= - ModuleBase::GlobalFunc::READ_VALUE(ifpos, ucell.atoms[it].label); - if(ucell.atoms[it].label != ucell.atom_label[it]) + bool set_element_mag_zero = false; + if (!unitcell::read_atom_type_header(it, ucell, ifpos, ofs_running, + ofs_warning, set_element_mag_zero)) { - ofs_warning << " Label orders in ATOMIC_POSITIONS and ATOMIC_SPECIES sections do not match!" << std::endl; - ofs_warning << " Label read from ATOMIC_POSITIONS is " << ucell.atoms[it].label << std::endl; - ofs_warning << " Label from ATOMIC_SPECIES is " << ucell.atom_label[it] << std::endl; return false; } - ModuleBase::GlobalFunc::OUT(ofs_running, "Atom label", ucell.atoms[it].label); - - bool set_element_mag_zero = false; - ModuleBase::GlobalFunc::READ_VALUE(ifpos, ucell.magnet.start_mag[it]); - -#ifndef __SYMMETRY - //=========================================== - // (2) read in numerical orbital information - // int ucell.atoms[it].nwl - // int* ucell.atoms[it].l_nchi; - //=========================================== - - if ((PARAM.inp.basis_type == "lcao")||(PARAM.inp.basis_type == "lcao_in_pw")) - { - std::string orbital_file = PARAM.inp.orbital_dir + ucell.orbital_fn[it]; - bool normal = elecstate::read_orb_file(it, orbital_file, ofs_running, &(ucell.atoms[it])); - if(!normal) - { - return false; - } - } - else if(PARAM.inp.basis_type == "pw") - { - if ((PARAM.inp.init_wfc.substr(0, 3) == "nao") || PARAM.inp.onsite_radius > 0.0) - { - std::string orbital_file = PARAM.inp.orbital_dir + ucell.orbital_fn[it]; - bool normal = elecstate::read_orb_file(it, orbital_file, ofs_running, &(ucell.atoms[it])); - if(!normal) - { - return false; - } - } - else - { - ucell.atoms[it].nw = 0; - ucell.atoms[it].nwl = 2; - if ( ucell.lmaxmax != 2 ) - { - ucell.atoms[it].nwl = ucell.lmaxmax; - } - ucell.atoms[it].l_nchi.resize(ucell.atoms[it].nwl+1, 0); - for(int L=0; L 0) { - ModuleBase::WARNING("read_atom_positions", " atom number < 0."); - return false; - } - else if (na == 0) - { - std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - std::cout << " Warning: atom number is 0 for atom type: " << ucell.atoms[it].label << std::endl; - std::cout << " If you are confident that this is not a mistake, please ignore this warning." << std::endl; - std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - ofs_running << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - ofs_running << " Warning: atom number is 0 for atom type: " << ucell.atoms[it].label << std::endl; - ofs_running << " If you are confident that this is not a mistake, please ignore this warning." << std::endl; - ofs_running << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - } - else if (na > 0) - { - ucell.atoms[it].tau.resize(na, ModuleBase::Vector3(0,0,0)); - ucell.atoms[it].dis.resize(na, ModuleBase::Vector3(0,0,0)); - ucell.atoms[it].taud.resize(na, ModuleBase::Vector3(0,0,0)); - ucell.atoms[it].boundary_shift.resize(na, ModuleBase::Vector3(0,0,0)); - ucell.atoms[it].vel.resize(na, ModuleBase::Vector3(0,0,0)); - ucell.atoms[it].mbl.resize(na, ModuleBase::Vector3(0,0,0)); - ucell.atoms[it].mag.resize(na, 0); - ucell.atoms[it].angle1.resize(na, 0); - ucell.atoms[it].angle2.resize(na, 0); - ucell.atoms[it].m_loc_.resize(na, ModuleBase::Vector3(0,0,0)); - ucell.atoms[it].lambda.resize(na, ModuleBase::Vector3(0,0,0)); - ucell.atoms[it].constrain.resize(na, ModuleBase::Vector3(0,0,0)); - ucell.atoms[it].mass = ucell.atom_mass[it]; //mohan add 2011-11-07 + unitcell::allocate_atom_properties(ucell.atoms[it], na, ucell.atom_mass[it]); for (int ia = 0;ia < na; ia++) { // modify the reading of frozen ions and velocities -- Yuanbo Li 2021/8/20 + ModuleBase::Vector3 v; + ModuleBase::Vector3 mv; ifpos >> v.x >> v.y >> v.z; mv.x = true ; mv.y = true ; @@ -199,305 +81,27 @@ bool unitcell::read_atom_positions(UnitCell& ucell, ucell.atoms[it].lambda[ia].set(0,0,0); ucell.atoms[it].constrain[ia].set(0,0,0); - std::string tmpid; - tmpid = ifpos.get(); - - if( (int)tmpid[0] < 0 ) - { - std::cout << "read_atom_positions, mismatch in atom number for atom type: " - << ucell.atoms[it].label << std::endl; - exit(1); - } - bool input_vec_mag=false; bool input_angle_mag=false; - // read if catch goodbit before "\n" and "#" - while ( (tmpid != "\n") && (ifpos.good()) && (tmpid !="#") ) - { - tmpid = ifpos.get() ; - // old method of reading frozen ions - char tmp = (char)tmpid[0]; - if ( tmp >= 48 && tmp <= 57 ) - { - mv.x = std::stoi(tmpid); - ifpos >> mv.y >> mv.z ; - } - // new method of reading frozen ions and velocities - if ( tmp >= 'a' && tmp <='z') - { - ifpos.putback(tmp); - ifpos >> tmpid; - } - if ( tmpid == "m" ) - { - ifpos >> mv.x >> mv.y >> mv.z ; - } - else if ( tmpid == "v" ||tmpid == "vel" || tmpid == "velocity" ) - { - ifpos >> ucell.atoms[it].vel[ia].x >> ucell.atoms[it].vel[ia].y >> ucell.atoms[it].vel[ia].z; - } - else if ( tmpid == "mag" || tmpid == "magmom") - { - set_element_mag_zero = true; - double tmpamg=0; - ifpos >> tmpamg; - tmp=ifpos.get(); - while (tmp==' ') - { - tmp=ifpos.get(); - } - - if((tmp >= 48 && tmp <= 57) or tmp=='-') - { - ifpos.putback(tmp); - ifpos >> ucell.atoms[it].m_loc_[ia].y>>ucell.atoms[it].m_loc_[ia].z; - ucell.atoms[it].m_loc_[ia].x=tmpamg; - ucell.atoms[it].mag[ia]=sqrt(pow(ucell.atoms[it].m_loc_[ia].x,2) - +pow(ucell.atoms[it].m_loc_[ia].y,2) - +pow(ucell.atoms[it].m_loc_[ia].z,2)); - input_vec_mag=true; - - } - else - { - ifpos.putback(tmp); - ucell.atoms[it].mag[ia]=tmpamg; - } - } - else if ( tmpid == "angle1") - { - ifpos >> ucell.atoms[it].angle1[ia]; - ucell.atoms[it].angle1[ia]=ucell.atoms[it].angle1[ia]/180 *ModuleBase::PI; - input_angle_mag=true; - set_element_mag_zero = true; - } - else if ( tmpid == "angle2") - { - ifpos >> ucell.atoms[it].angle2[ia]; - ucell.atoms[it].angle2[ia]=ucell.atoms[it].angle2[ia]/180 *ModuleBase::PI; - input_angle_mag=true; - set_element_mag_zero = true; - } - else if ( tmpid == "lambda") - { - double tmplam=0; - ifpos >> tmplam; - tmp=ifpos.get(); - while (tmp==' ') - { - tmp=ifpos.get(); - } - if((tmp >= 48 && tmp <= 57) or tmp=='-') - { - ifpos.putback(tmp); - ifpos >> ucell.atoms[it].lambda[ia].y>>ucell.atoms[it].lambda[ia].z; - ucell.atoms[it].lambda[ia].x=tmplam; - } - else - { - ifpos.putback(tmp); - ucell.atoms[it].lambda[ia].z=tmplam; - } - ucell.atoms[it].lambda[ia].x /= ModuleBase::Ry_to_eV; - ucell.atoms[it].lambda[ia].y /= ModuleBase::Ry_to_eV; - ucell.atoms[it].lambda[ia].z /= ModuleBase::Ry_to_eV; - } - else if ( tmpid == "sc") - { - double tmplam=0; - ifpos >> tmplam; - tmp=ifpos.get(); - while (tmp==' ') - { - tmp=ifpos.get(); - } - if((tmp >= 48 && tmp <= 57) or tmp=='-') - { - ifpos.putback(tmp); - ifpos >> ucell.atoms[it].constrain[ia].y>>ucell.atoms[it].constrain[ia].z; - ucell.atoms[it].constrain[ia].x=tmplam; - } - else - { - ifpos.putback(tmp); - ucell.atoms[it].constrain[ia].z=tmplam; - } - } - } - // move to next line - while ( (tmpid != "\n") && (ifpos.good()) ) - { - tmpid = ifpos.get(); - } - std::string mags; - - // ---------------------------------------------------------------------------- - // recalcualte mag and m_loc_ from read in angle1, angle2 and mag or mx, my, mz - if(input_angle_mag) - {// angle1 or angle2 are given, calculate mx, my, mz from angle1 and angle2 and mag - ucell.atoms[it].m_loc_[ia].z = ucell.atoms[it].mag[ia] * - cos(ucell.atoms[it].angle1[ia]); - if(std::abs(sin(ucell.atoms[it].angle1[ia])) > 1e-10 ) - { - ucell.atoms[it].m_loc_[ia].x = ucell.atoms[it].mag[ia] * - sin(ucell.atoms[it].angle1[ia]) * cos(ucell.atoms[it].angle2[ia]); - ucell.atoms[it].m_loc_[ia].y = ucell.atoms[it].mag[ia] * - sin(ucell.atoms[it].angle1[ia]) * sin(ucell.atoms[it].angle2[ia]); - } - } - else if (input_vec_mag) - {// mx, my, mz are given, calculate angle1 and angle2 from mx, my, mz - double mxy=sqrt(pow(ucell.atoms[it].m_loc_[ia].x,2)+pow(ucell.atoms[it].m_loc_[ia].y,2)); - ucell.atoms[it].angle1[ia]=atan2(mxy,ucell.atoms[it].m_loc_[ia].z); - if(mxy>1e-8) - { - ucell.atoms[it].angle2[ia]=atan2(ucell.atoms[it].m_loc_[ia].y,ucell.atoms[it].m_loc_[ia].x); - } - } - else// only one mag is given, assume it is z + // Parse optional properties + if (!unitcell::parse_atom_properties(ifpos, ucell.atoms[it], ia, mv, + input_vec_mag, input_angle_mag, + set_element_mag_zero)) { - ucell.atoms[it].m_loc_[ia].x = 0; - ucell.atoms[it].m_loc_[ia].y = 0; - ucell.atoms[it].m_loc_[ia].z = ucell.atoms[it].mag[ia]; + return false; } - if(nspin==4) - { - if(!PARAM.inp.noncolin) - { - //collinear case with nspin = 4, only z component is used - ucell.atoms[it].m_loc_[ia].x = 0; - ucell.atoms[it].m_loc_[ia].y = 0; - } - //print only ia==0 && mag>0 to avoid too much output - //print when ia!=0 && mag[ia] != mag[0] to avoid too much output - // 'A || (!A && B)' is equivalent to 'A || B',so the following - // code is equivalent to 'ia==0 || (...)' - if(ia==0 || (ucell.atoms[it].m_loc_[ia].x != ucell.atoms[it].m_loc_[0].x - || ucell.atoms[it].m_loc_[ia].y != ucell.atoms[it].m_loc_[0].y - || ucell.atoms[it].m_loc_[ia].z != ucell.atoms[it].m_loc_[0].z)) - { - //use a stringstream to generate string: "concollinear magnetization of element it is:" - std::stringstream ss; - ss << "Magnetization for this type"; - if(ia!=0) - { - ss<<" (atom"<0 to avoid too much output - //print when ia!=0 && mag[ia] != mag[0] to avoid too much output - if(ia==0 || (ucell.atoms[it].mag[ia] != ucell.atoms[it].mag[0])) - { - //use a stringstream to generate string: "cocollinear magnetization of element it is:" - std::stringstream ss; - ss << "magnetization of element " << it+1; - if(ia!=0) - { - ss<<" (atom"< 1e-5) - { - autoset_mag = 0; - break; - } - } - } - if (autoset_mag) - { - if(nspin==4) - { - for (int it = 0;it < ntype; it++) - { - for (int ia = 0;ia < ucell.atoms[it].na; ia++) - { - ucell.atoms[it].m_loc_[ia].x = 1.0; - ucell.atoms[it].m_loc_[ia].y = 1.0; - ucell.atoms[it].m_loc_[ia].z = 1.0; - ucell.atoms[it].mag[ia] = sqrt(pow(ucell.atoms[it].m_loc_[ia].x,2) - +pow(ucell.atoms[it].m_loc_[ia].y,2) - +pow(ucell.atoms[it].m_loc_[ia].z,2)); - ModuleBase::GlobalFunc::OUT(ofs_running,"Autoset magnetism for this atom", 1.0, 1.0, 1.0); - } - } - } - else if(nspin==2) - { - for (int it = 0;it < ntype; it++) - { - for (int ia = 0;ia < ucell.atoms[it].na; ia++) - { - ucell.atoms[it].mag[ia] = 1.0; - ucell.atoms[it].m_loc_[ia].x = ucell.atoms[it].mag[ia]; - ModuleBase::GlobalFunc::OUT(ofs_running,"Autoset magnetism for this atom", 1.0); - } - } - } - } - // End Autoset magnetization + // Auto-set magnetization if needed + unitcell::autoset_magnetization(ucell, nspin, ofs_running); } // end scan_begin - //check if any atom can move in MD - if(!ucell.if_atoms_can_move() && PARAM.inp.calculation=="md" && PARAM.inp.esolver_type!="tddft") - { - ModuleBase::WARNING("read_atoms", "no atoms can move in MD simulations!"); - return false; - } - - ofs_running << std::endl; - ModuleBase::GlobalFunc::OUT(ofs_running,"TOTAL ATOM NUMBER",ucell.nat); - ofs_running << std::endl; - - if (ucell.nat == 0) - { - ModuleBase::WARNING("read_atom_positions","no atoms found in the system!"); - return false; - } - - // mohan add 2010-06-30 - unitcell::check_dtau(ucell.atoms,ucell.ntype, ucell.lat0, ucell.latvec); - - if (unitcell::check_tau(ucell.atoms, ucell.ntype, ucell.lat0)) - { - unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,ofs_running); - return true; - } - return false; + // Final validation and output + return unitcell::finalize_atom_positions(ucell, ofs_running, ofs_warning); }//end read_atom_positions diff --git a/source/source_cell/read_atoms_helper.cpp b/source/source_cell/read_atoms_helper.cpp new file mode 100644 index 0000000000..4fc3dfe6cb --- /dev/null +++ b/source/source_cell/read_atoms_helper.cpp @@ -0,0 +1,596 @@ +#include "read_atoms_helper.h" +#include "source_io/module_parameter/parameter.h" +#include "source_base/global_function.h" +#include "source_base/constants.h" +#include "source_base/mathzone.h" +#include "read_stru.h" +#include "print_cell.h" +#include "source_estate/read_orb.h" +#include +#include +#include + +namespace { + // Magic number constants for character code checks + constexpr char DIGIT_START = '0'; // ASCII 48 + constexpr char DIGIT_END = '9'; // ASCII 57 + constexpr char LOWER_A = 'a'; + constexpr char LOWER_Z = 'z'; + constexpr char MINUS_SIGN = '-'; +} + +namespace unitcell { + +bool validate_coordinate_system(const std::string& Coordinate, + std::ofstream& ofs_warning) +{ + if(Coordinate != "Cartesian" + && Coordinate != "Direct" + && Coordinate != "Cartesian_angstrom" + && Coordinate != "Cartesian_au" + && Coordinate != "Cartesian_angstrom_center_xy" + && Coordinate != "Cartesian_angstrom_center_xz" + && Coordinate != "Cartesian_angstrom_center_yz" + && Coordinate != "Cartesian_angstrom_center_xyz" + ) + { + ModuleBase::WARNING("read_atom_position","Cartesian or Direct?"); + ofs_warning << " There are several options for you:" << std::endl; + ofs_warning << " Direct" << std::endl; + ofs_warning << " Cartesian_angstrom" << std::endl; + ofs_warning << " Cartesian_au" << std::endl; + ofs_warning << " Cartesian_angstrom_center_xy" << std::endl; + ofs_warning << " Cartesian_angstrom_center_xz" << std::endl; + ofs_warning << " Cartesian_angstrom_center_yz" << std::endl; + ofs_warning << " Cartesian_angstrom_center_xyz" << std::endl; + return false; + } + return true; +} + +void allocate_atom_properties(Atom& atom, int na, double mass) +{ + atom.tau.resize(na, ModuleBase::Vector3(0,0,0)); + atom.dis.resize(na, ModuleBase::Vector3(0,0,0)); + atom.taud.resize(na, ModuleBase::Vector3(0,0,0)); + atom.boundary_shift.resize(na, ModuleBase::Vector3(0,0,0)); + atom.vel.resize(na, ModuleBase::Vector3(0,0,0)); + atom.mbl.resize(na, ModuleBase::Vector3(0,0,0)); + atom.mag.resize(na, 0); + atom.angle1.resize(na, 0); + atom.angle2.resize(na, 0); + atom.m_loc_.resize(na, ModuleBase::Vector3(0,0,0)); + atom.lambda.resize(na, ModuleBase::Vector3(0,0,0)); + atom.constrain.resize(na, ModuleBase::Vector3(0,0,0)); + atom.mass = mass; +} + +void set_atom_movement_flags(Atom& atom, int ia, + const ModuleBase::Vector3& mv) +{ + if(!PARAM.inp.fixed_atoms) + { + atom.mbl[ia] = mv; + } + else + { + atom.mbl[ia] = 0.0; + atom.mbl[ia].print(); + } +} + +void autoset_magnetization(UnitCell& ucell, int nspin, + std::ofstream& ofs_running) +{ + const int ntype = ucell.ntype; + + // Check if any atom has non-zero magnetization + int autoset_mag = 1; + for (int it = 0; it < ntype; it++) + { + for (int ia = 0; ia < ucell.atoms[it].na; ia++) + { + if(std::abs(ucell.atoms[it].mag[ia]) > 1e-5) + { + autoset_mag = 0; + break; + } + } + } + + if (autoset_mag) + { + if(nspin==4) + { + for (int it = 0; it < ntype; it++) + { + for (int ia = 0; ia < ucell.atoms[it].na; ia++) + { + ucell.atoms[it].m_loc_[ia].x = 1.0; + ucell.atoms[it].m_loc_[ia].y = 1.0; + ucell.atoms[it].m_loc_[ia].z = 1.0; + ucell.atoms[it].mag[ia] = sqrt(pow(ucell.atoms[it].m_loc_[ia].x,2) + +pow(ucell.atoms[it].m_loc_[ia].y,2) + +pow(ucell.atoms[it].m_loc_[ia].z,2)); + ModuleBase::GlobalFunc::OUT(ofs_running,"Autoset magnetism for this atom", 1.0, 1.0, 1.0); + } + } + } + else if(nspin==2) + { + for (int it = 0; it < ntype; it++) + { + for (int ia = 0; ia < ucell.atoms[it].na; ia++) + { + ucell.atoms[it].mag[ia] = 1.0; + ucell.atoms[it].m_loc_[ia].x = ucell.atoms[it].mag[ia]; + ModuleBase::GlobalFunc::OUT(ofs_running,"Autoset magnetism for this atom", 1.0); + } + } + } + } +} + +bool finalize_atom_positions(UnitCell& ucell, + std::ofstream& ofs_running, + std::ofstream& ofs_warning) +{ + // Check if any atom can move in MD + if(!ucell.if_atoms_can_move() && PARAM.inp.calculation=="md" && PARAM.inp.esolver_type!="tddft") + { + ModuleBase::WARNING("read_atoms", "no atoms can move in MD simulations!"); + return false; + } + + ofs_running << std::endl; + ModuleBase::GlobalFunc::OUT(ofs_running,"TOTAL ATOM NUMBER",ucell.nat); + ofs_running << std::endl; + + if (ucell.nat == 0) + { + ModuleBase::WARNING("read_atom_positions","no atoms found in the system!"); + return false; + } + + // Check atom positions + unitcell::check_dtau(ucell.atoms, ucell.ntype, ucell.lat0, ucell.latvec); + + if (unitcell::check_tau(ucell.atoms, ucell.ntype, ucell.lat0)) + { + print_tau(ucell.atoms, ucell.Coordinate, ucell.ntype, ucell.lat0, ofs_running); + return true; + } + return false; +} + +ModuleBase::Vector3 calculate_lattice_center( + const ModuleBase::Matrix3& latvec, + const std::string& center_mode) +{ + ModuleBase::Vector3 latcenter(0.0, 0.0, 0.0); + + if (center_mode == "xy" || center_mode == "xyz") + { + latcenter.x = (latvec.e11 + latvec.e21 + latvec.e31) / 2.0; + latcenter.y = (latvec.e12 + latvec.e22 + latvec.e32) / 2.0; + } + + if (center_mode == "xz" || center_mode == "xyz") + { + latcenter.x = (latvec.e11 + latvec.e21 + latvec.e31) / 2.0; + latcenter.z = (latvec.e13 + latvec.e23 + latvec.e33) / 2.0; + } + + if (center_mode == "yz" || center_mode == "xyz") + { + latcenter.y = (latvec.e12 + latvec.e22 + latvec.e32) / 2.0; + latcenter.z = (latvec.e13 + latvec.e23 + latvec.e33) / 2.0; + } + + return latcenter; +} + +void transform_atom_coordinates(Atom& atom, int ia, + const std::string& Coordinate, + const ModuleBase::Vector3& v, + const ModuleBase::Matrix3& latvec, + double lat0, + ModuleBase::Vector3& latcenter) +{ + if(Coordinate=="Direct") + { + // change v from direct to cartesian, + // the unit is GlobalC::sf.ucell.lat0 + atom.taud[ia] = v; + atom.tau[ia] = v * latvec; + } + else if(Coordinate=="Cartesian") + { + atom.tau[ia] = v; // in unit ucell.lat0 + } + else if(Coordinate=="Cartesian_angstrom") + { + atom.tau[ia] = v / ModuleBase::BOHR_TO_A / lat0; + } + else if(Coordinate=="Cartesian_angstrom_center_xy") + { + latcenter = calculate_lattice_center(latvec, "xy"); + atom.tau[ia] = v / ModuleBase::BOHR_TO_A / lat0 + latcenter; + } + else if(Coordinate=="Cartesian_angstrom_center_xz") + { + latcenter = calculate_lattice_center(latvec, "xz"); + atom.tau[ia] = v / ModuleBase::BOHR_TO_A / lat0 + latcenter; + } + else if(Coordinate=="Cartesian_angstrom_center_yz") + { + latcenter = calculate_lattice_center(latvec, "yz"); + atom.tau[ia] = v / ModuleBase::BOHR_TO_A / lat0 + latcenter; + } + else if(Coordinate=="Cartesian_angstrom_center_xyz") + { + latcenter = calculate_lattice_center(latvec, "xyz"); + atom.tau[ia] = v / ModuleBase::BOHR_TO_A / lat0 + latcenter; + } + else if(Coordinate=="Cartesian_au") + { + atom.tau[ia] = v / lat0; + } + + // Convert to direct coordinates if using Cartesian + if(Coordinate=="Cartesian" || + Coordinate=="Cartesian_angstrom" || + Coordinate=="Cartesian_angstrom_center_xy" || + Coordinate=="Cartesian_angstrom_center_xz" || + Coordinate=="Cartesian_angstrom_center_yz" || + Coordinate=="Cartesian_angstrom_center_xyz" || + Coordinate=="Cartesian_au") + { + double dx=0.0; + double dy=0.0; + double dz=0.0; + ModuleBase::Mathzone::Cartesian_to_Direct(atom.tau[ia].x, + atom.tau[ia].y, + atom.tau[ia].z, + latvec.e11, latvec.e12, latvec.e13, + latvec.e21, latvec.e22, latvec.e23, + latvec.e31, latvec.e32, latvec.e33, + dx,dy,dz); + + atom.taud[ia].x = dx; + atom.taud[ia].y = dy; + atom.taud[ia].z = dz; + } +} + +void process_magnetization(Atom& atom, int it, int ia, + int nspin, bool input_vec_mag, + bool input_angle_mag, + std::ofstream& ofs_running) +{ + // Recalculate mag and m_loc_ from read in angle1, angle2 and mag or mx, my, mz + if(input_angle_mag) + { + // angle1 or angle2 are given, calculate mx, my, mz from angle1 and angle2 and mag + atom.m_loc_[ia].z = atom.mag[ia] * cos(atom.angle1[ia]); + if(std::abs(sin(atom.angle1[ia])) > 1e-10) + { + atom.m_loc_[ia].x = atom.mag[ia] * + sin(atom.angle1[ia]) * cos(atom.angle2[ia]); + atom.m_loc_[ia].y = atom.mag[ia] * + sin(atom.angle1[ia]) * sin(atom.angle2[ia]); + } + } + else if (input_vec_mag) + { + // mx, my, mz are given, calculate angle1 and angle2 from mx, my, mz + double mxy=sqrt(pow(atom.m_loc_[ia].x,2)+pow(atom.m_loc_[ia].y,2)); + atom.angle1[ia]=atan2(mxy,atom.m_loc_[ia].z); + if(mxy>1e-8) + { + atom.angle2[ia]=atan2(atom.m_loc_[ia].y,atom.m_loc_[ia].x); + } + } + else + { + // only one mag is given, assume it is z + atom.m_loc_[ia].x = 0; + atom.m_loc_[ia].y = 0; + atom.m_loc_[ia].z = atom.mag[ia]; + } + + if(nspin==4) + { + if(!PARAM.inp.noncolin) + { + // collinear case with nspin = 4, only z component is used + atom.m_loc_[ia].x = 0; + atom.m_loc_[ia].y = 0; + } + // print only ia==0 && mag>0 to avoid too much output + // print when ia!=0 && mag[ia] != mag[0] to avoid too much output + if(ia==0 || (atom.m_loc_[ia].x != atom.m_loc_[0].x + || atom.m_loc_[ia].y != atom.m_loc_[0].y + || atom.m_loc_[ia].z != atom.m_loc_[0].z)) + { + std::stringstream ss; + ss << "Magnetization for this type"; + if(ia!=0) + { + ss<<" (atom"<0 to avoid too much output + // print when ia!=0 && mag[ia] != mag[0] to avoid too much output + if(ia==0 || (atom.mag[ia] != atom.mag[0])) + { + std::stringstream ss; + ss << "magnetization of element " << it+1; + if(ia!=0) + { + ss<<" (atom"<& mv, + bool& input_vec_mag, + bool& input_angle_mag, + bool& set_element_mag_zero) +{ + std::string tmpid; + tmpid = ifpos.get(); + + if( (int)tmpid[0] < 0 ) + { + std::cout << "read_atom_positions, mismatch in atom number for atom type: " + << atom.label << std::endl; + exit(1); + } + + // read if catch goodbit before "\n" and "#" + while ( (tmpid != "\n") && (ifpos.good()) && (tmpid !="#") ) + { + tmpid = ifpos.get(); + // old method of reading frozen ions + char tmp = (char)tmpid[0]; + if ( tmp >= DIGIT_START && tmp <= DIGIT_END ) + { + mv.x = std::stoi(tmpid); + ifpos >> mv.y >> mv.z; + } + // new method of reading frozen ions and velocities + if ( tmp >= LOWER_A && tmp <= LOWER_Z) + { + ifpos.putback(tmp); + ifpos >> tmpid; + } + if ( tmpid == "m" ) + { + ifpos >> mv.x >> mv.y >> mv.z; + } + else if ( tmpid == "v" ||tmpid == "vel" || tmpid == "velocity" ) + { + ifpos >> atom.vel[ia].x >> atom.vel[ia].y >> atom.vel[ia].z; + } + else if ( tmpid == "mag" || tmpid == "magmom") + { + set_element_mag_zero = true; + double tmpamg=0; + ifpos >> tmpamg; + tmp=ifpos.get(); + while (tmp==' ') + { + tmp=ifpos.get(); + } + + if((tmp >= DIGIT_START && tmp <= DIGIT_END) or tmp==MINUS_SIGN) + { + ifpos.putback(tmp); + ifpos >> atom.m_loc_[ia].y>>atom.m_loc_[ia].z; + atom.m_loc_[ia].x=tmpamg; + atom.mag[ia]=sqrt(pow(atom.m_loc_[ia].x,2) + +pow(atom.m_loc_[ia].y,2) + +pow(atom.m_loc_[ia].z,2)); + input_vec_mag=true; + + } + else + { + ifpos.putback(tmp); + atom.mag[ia]=tmpamg; + } + } + else if ( tmpid == "angle1") + { + ifpos >> atom.angle1[ia]; + atom.angle1[ia]=atom.angle1[ia]/180 *ModuleBase::PI; + input_angle_mag=true; + set_element_mag_zero = true; + } + else if ( tmpid == "angle2") + { + ifpos >> atom.angle2[ia]; + atom.angle2[ia]=atom.angle2[ia]/180 *ModuleBase::PI; + input_angle_mag=true; + set_element_mag_zero = true; + } + else if ( tmpid == "lambda") + { + double tmplam=0; + ifpos >> tmplam; + tmp=ifpos.get(); + while (tmp==' ') + { + tmp=ifpos.get(); + } + if((tmp >= DIGIT_START && tmp <= DIGIT_END) or tmp==MINUS_SIGN) + { + ifpos.putback(tmp); + ifpos >> atom.lambda[ia].y>>atom.lambda[ia].z; + atom.lambda[ia].x=tmplam; + } + else + { + ifpos.putback(tmp); + atom.lambda[ia].z=tmplam; + } + atom.lambda[ia].x /= ModuleBase::Ry_to_eV; + atom.lambda[ia].y /= ModuleBase::Ry_to_eV; + atom.lambda[ia].z /= ModuleBase::Ry_to_eV; + } + else if ( tmpid == "sc") + { + double tmplam=0; + ifpos >> tmplam; + tmp=ifpos.get(); + while (tmp==' ') + { + tmp=ifpos.get(); + } + if((tmp >= DIGIT_START && tmp <= DIGIT_END) or tmp==MINUS_SIGN) + { + ifpos.putback(tmp); + ifpos >> atom.constrain[ia].y>>atom.constrain[ia].z; + atom.constrain[ia].x=tmplam; + } + else + { + ifpos.putback(tmp); + atom.constrain[ia].z=tmplam; + } + } + } + // move to next line + while ( (tmpid != "\n") && (ifpos.good()) ) + { + tmpid = ifpos.get(); + } + + return true; +} + +bool read_atom_type_header(int it, UnitCell& ucell, + std::ifstream& ifpos, + std::ofstream& ofs_running, + std::ofstream& ofs_warning, + bool& set_element_mag_zero) +{ + //======================================= + // (1) read in atom label + // start magnetization + //======================================= + ModuleBase::GlobalFunc::READ_VALUE(ifpos, ucell.atoms[it].label); + + if(ucell.atoms[it].label != ucell.atom_label[it]) + { + ofs_warning << " Label orders in ATOMIC_POSITIONS and ATOMIC_SPECIES sections do not match!" << std::endl; + ofs_warning << " Label read from ATOMIC_POSITIONS is " << ucell.atoms[it].label << std::endl; + ofs_warning << " Label from ATOMIC_SPECIES is " << ucell.atom_label[it] << std::endl; + return false; + } + ModuleBase::GlobalFunc::OUT(ofs_running, "Atom label", ucell.atoms[it].label); + + set_element_mag_zero = false; + ModuleBase::GlobalFunc::READ_VALUE(ifpos, ucell.magnet.start_mag[it]); + +#ifndef __SYMMETRY + //=========================================== + // (2) read in numerical orbital information + // int ucell.atoms[it].nwl + // int* ucell.atoms[it].l_nchi; + //=========================================== + + if ((PARAM.inp.basis_type == "lcao")||(PARAM.inp.basis_type == "lcao_in_pw")) + { + std::string orbital_file = PARAM.inp.orbital_dir + ucell.orbital_fn[it]; + bool normal = elecstate::read_orb_file(it, orbital_file, ofs_running, &(ucell.atoms[it])); + if(!normal) + { + return false; + } + } + else if(PARAM.inp.basis_type == "pw") + { + if ((PARAM.inp.init_wfc.substr(0, 3) == "nao") || PARAM.inp.onsite_radius > 0.0) + { + std::string orbital_file = PARAM.inp.orbital_dir + ucell.orbital_fn[it]; + bool normal = elecstate::read_orb_file(it, orbital_file, ofs_running, &(ucell.atoms[it])); + if(!normal) + { + return false; + } + } + else + { + ucell.atoms[it].nw = 0; + ucell.atoms[it].nwl = 2; + if ( ucell.lmaxmax != 2 ) + { + ucell.atoms[it].nwl = ucell.lmaxmax; + } + ucell.atoms[it].l_nchi.resize(ucell.atoms[it].nwl+1, 0); + for(int L=0; L +#include +#include "unitcell.h" +#include "source_base/vector3.h" +#include "source_base/matrix3.h" + +namespace unitcell { + +/** + * @brief Validate coordinate system type + * @param Coordinate The coordinate system string to validate + * @param ofs_warning Output stream for warnings + * @return true if valid, false otherwise + */ +bool validate_coordinate_system(const std::string& Coordinate, + std::ofstream& ofs_warning); + +/** + * @brief Allocate and initialize atom property vectors + * @param atom The atom object to allocate properties for + * @param na Number of atoms + * @param mass Atomic mass + */ +void allocate_atom_properties(Atom& atom, int na, double mass); + +/** + * @brief Set atom movement constraints based on fixed_atoms parameter + * @param atom The atom object + * @param ia Atom index + * @param mv Movement vector (1=movable, 0=fixed) + */ +void set_atom_movement_flags(Atom& atom, int ia, + const ModuleBase::Vector3& mv); + +/** + * @brief Set default magnetization if not explicitly specified + * @param ucell Unit cell object + * @param nspin Number of spin components + * @param ofs_running Output stream for running information + */ +void autoset_magnetization(UnitCell& ucell, int nspin, + std::ofstream& ofs_running); + +/** + * @brief Perform final validation and output + * @param ucell Unit cell object + * @param ofs_running Output stream for running information + * @param ofs_warning Output stream for warnings + * @return true if validation passes, false otherwise + */ +bool finalize_atom_positions(UnitCell& ucell, + std::ofstream& ofs_running, + std::ofstream& ofs_warning); + +/** + * @brief Calculate lattice center for different centering modes + * @param latvec Lattice vectors + * @param center_mode Centering mode: "xy", "xz", "yz", or "xyz" + * @return Lattice center coordinates + */ +ModuleBase::Vector3 calculate_lattice_center( + const ModuleBase::Matrix3& latvec, + const std::string& center_mode); + +/** + * @brief Convert between different coordinate systems + * @param atom The atom object + * @param ia Atom index + * @param Coordinate Coordinate system type + * @param v Input position vector + * @param latvec Lattice vectors + * @param lat0 Lattice constant + * @param latcenter Lattice center (output parameter) + */ +void transform_atom_coordinates(Atom& atom, int ia, + const std::string& Coordinate, + const ModuleBase::Vector3& v, + const ModuleBase::Matrix3& latvec, + double lat0, + ModuleBase::Vector3& latcenter); + +/** + * @brief Convert between magnetization representations and output + * @param atom The atom object + * @param it Atom type index + * @param ia Atom index + * @param nspin Number of spin components + * @param input_vec_mag Whether vector magnetization was input + * @param input_angle_mag Whether angle magnetization was input + * @param ofs_running Output stream for running information + */ +void process_magnetization(Atom& atom, int it, int ia, + int nspin, bool input_vec_mag, + bool input_angle_mag, + std::ofstream& ofs_running); + +/** + * @brief Parse optional atom properties (mag, angle1, angle2, lambda, sc, m, v) + * @param ifpos Input file stream + * @param atom The atom object + * @param ia Atom index + * @param mv Movement vector (output parameter) + * @param input_vec_mag Whether vector magnetization was input (output parameter) + * @param input_angle_mag Whether angle magnetization was input (output parameter) + * @param set_element_mag_zero Whether to reset element magnetization (output parameter) + * @return true if parsing succeeds, false otherwise + */ +bool parse_atom_properties(std::ifstream& ifpos, + Atom& atom, int ia, + ModuleBase::Vector3& mv, + bool& input_vec_mag, + bool& input_angle_mag, + bool& set_element_mag_zero); + +/** + * @brief Read atom type metadata (label, magnetization, orbital info, atom count) + * @param it Atom type index + * @param ucell Unit cell object + * @param ifpos Input file stream + * @param ofs_running Output stream for running information + * @param ofs_warning Output stream for warnings + * @param set_element_mag_zero Whether to reset element magnetization (output parameter) + * @return true if reading succeeds, false otherwise + */ +bool read_atom_type_header(int it, UnitCell& ucell, + std::ifstream& ifpos, + std::ofstream& ofs_running, + std::ofstream& ofs_warning, + bool& set_element_mag_zero); + +} // namespace unitcell + +#endif // READ_ATOMS_HELPER_H diff --git a/source/source_cell/test/CMakeLists.txt b/source/source_cell/test/CMakeLists.txt index 9e8ba3ce6f..708679c32c 100644 --- a/source/source_cell/test/CMakeLists.txt +++ b/source/source_cell/test/CMakeLists.txt @@ -21,6 +21,7 @@ list(APPEND cell_simple_srcs ../read_stru.cpp ../read_atom_species.cpp ../read_atoms.cpp + ../read_atoms_helper.cpp ../atom_spec.cpp ../atom_pseudo.cpp ../pseudo.cpp @@ -93,6 +94,31 @@ AddTest( ../../source_base/parallel_common.cpp ../../source_base/parallel_comm.cpp ../parallel_kpoints.cpp ../../source_base/tool_quit.cpp ../../source_base/global_variable.cpp ../../source_base/global_file.cpp ../../source_base/global_function.cpp ../../source_base/memory.cpp ../../source_base/timer.cpp ../../source_base/parallel_reduce.cpp ) +# Add unit test for read_atoms_helper +AddTest( + TARGET MODULE_CELL_read_atoms_helper_test + LIBS parameter ${math_libs} base device + SOURCES read_atoms_helper_test.cpp + ../read_atoms_helper.cpp + ../read_stru.cpp + ../print_cell.cpp + ../atom_spec.cpp + ../unitcell.cpp + ../update_cell.cpp + ../bcast_cell.cpp + ../atom_pseudo.cpp + ../pseudo.cpp + ../read_pp.cpp + ../read_pp_complete.cpp + ../read_pp_upf201.cpp + ../read_pp_upf100.cpp + ../read_pp_vwr.cpp + ../read_pp_blps.cpp + ../read_atom_species.cpp + ../sep.cpp + ../sep_cell.cpp +) + find_program(BASH bash) add_test(NAME MODULE_CELL_bcast_atom_pseudo_test diff --git a/source/source_cell/test/read_atoms_helper_test.cpp b/source/source_cell/test/read_atoms_helper_test.cpp new file mode 100644 index 0000000000..1cd8084784 --- /dev/null +++ b/source/source_cell/test/read_atoms_helper_test.cpp @@ -0,0 +1,483 @@ +#include "gtest/gtest.h" +#include "gmock/gmock.h" +#include "../read_atoms_helper.h" +#include "source_base/vector3.h" +#include "source_base/matrix3.h" +#include +#include + +// Mock implementations for missing functions that are not in the linked sources +namespace elecstate { + bool read_orb_file(int it, std::string& orbital_file, std::ofstream& ofs_running, Atom* atom) { + // Mock implementation - just return true + return true; + } +} + +// Mock output class methods +void output::printM3(std::ofstream& ofs, const std::string& description, const ModuleBase::Matrix3& m) { + // Mock implementation +} + +void output::printrm(std::ofstream& ofs, const std::string& description, const ModuleBase::matrix& m, const double& limit) { + // Mock implementation +} + +// Mock InfoNonlocal class +InfoNonlocal::InfoNonlocal() {} +InfoNonlocal::~InfoNonlocal() {} + +// Mock Magnetism class +Magnetism::Magnetism() {} +Magnetism::~Magnetism() {} + +// Mock read_atom_positions function (we're testing the helpers, not the main function) +namespace unitcell { + bool read_atom_positions(UnitCell& ucell, std::ifstream& ifpos, + std::ofstream& ofs_running, std::ofstream& ofs_warning) { + // Mock implementation + return true; + } +} + +// Test fixture for read_atoms_helper tests +class ReadAtomsHelperTest : public ::testing::Test +{ +protected: + void SetUp() override + { + // Create temporary output streams + ofs_warning.open("test_warning.log"); + ofs_running.open("test_running.log"); + } + + void TearDown() override + { + ofs_warning.close(); + ofs_running.close(); + // Clean up temporary files + std::remove("test_warning.log"); + std::remove("test_running.log"); + } + + std::ofstream ofs_warning; + std::ofstream ofs_running; +}; + +// Test validate_coordinate_system function +TEST_F(ReadAtomsHelperTest, ValidateCoordinateSystem_ValidInputs) +{ + EXPECT_TRUE(unitcell::validate_coordinate_system("Direct", ofs_warning)); + EXPECT_TRUE(unitcell::validate_coordinate_system("Cartesian", ofs_warning)); + EXPECT_TRUE(unitcell::validate_coordinate_system("Cartesian_angstrom", ofs_warning)); + EXPECT_TRUE(unitcell::validate_coordinate_system("Cartesian_au", ofs_warning)); + EXPECT_TRUE(unitcell::validate_coordinate_system("Cartesian_angstrom_center_xy", ofs_warning)); + EXPECT_TRUE(unitcell::validate_coordinate_system("Cartesian_angstrom_center_xz", ofs_warning)); + EXPECT_TRUE(unitcell::validate_coordinate_system("Cartesian_angstrom_center_yz", ofs_warning)); + EXPECT_TRUE(unitcell::validate_coordinate_system("Cartesian_angstrom_center_xyz", ofs_warning)); +} + +TEST_F(ReadAtomsHelperTest, ValidateCoordinateSystem_InvalidInputs) +{ + EXPECT_FALSE(unitcell::validate_coordinate_system("Invalid", ofs_warning)); + EXPECT_FALSE(unitcell::validate_coordinate_system("direct", ofs_warning)); // case sensitive + EXPECT_FALSE(unitcell::validate_coordinate_system("", ofs_warning)); + EXPECT_FALSE(unitcell::validate_coordinate_system("Cartesian_angstrom_center", ofs_warning)); +} + +// Test calculate_lattice_center function +TEST_F(ReadAtomsHelperTest, CalculateLatticeCenterXY) +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = 10.0; latvec.e12 = 0.0; latvec.e13 = 0.0; + latvec.e21 = 0.0; latvec.e22 = 10.0; latvec.e23 = 0.0; + latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 10.0; + + auto center = unitcell::calculate_lattice_center(latvec, "xy"); + + EXPECT_DOUBLE_EQ(center.x, 5.0); + EXPECT_DOUBLE_EQ(center.y, 5.0); + EXPECT_DOUBLE_EQ(center.z, 0.0); +} + +TEST_F(ReadAtomsHelperTest, CalculateLatticeCenterXZ) +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = 10.0; latvec.e12 = 0.0; latvec.e13 = 0.0; + latvec.e21 = 0.0; latvec.e22 = 10.0; latvec.e23 = 0.0; + latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 10.0; + + auto center = unitcell::calculate_lattice_center(latvec, "xz"); + + EXPECT_DOUBLE_EQ(center.x, 5.0); + EXPECT_DOUBLE_EQ(center.y, 0.0); + EXPECT_DOUBLE_EQ(center.z, 5.0); +} + +TEST_F(ReadAtomsHelperTest, CalculateLatticeCenterYZ) +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = 10.0; latvec.e12 = 0.0; latvec.e13 = 0.0; + latvec.e21 = 0.0; latvec.e22 = 10.0; latvec.e23 = 0.0; + latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 10.0; + + auto center = unitcell::calculate_lattice_center(latvec, "yz"); + + EXPECT_DOUBLE_EQ(center.x, 0.0); + EXPECT_DOUBLE_EQ(center.y, 5.0); + EXPECT_DOUBLE_EQ(center.z, 5.0); +} + +TEST_F(ReadAtomsHelperTest, CalculateLatticeCenterXYZ) +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = 10.0; latvec.e12 = 0.0; latvec.e13 = 0.0; + latvec.e21 = 0.0; latvec.e22 = 10.0; latvec.e23 = 0.0; + latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 10.0; + + auto center = unitcell::calculate_lattice_center(latvec, "xyz"); + + EXPECT_DOUBLE_EQ(center.x, 5.0); + EXPECT_DOUBLE_EQ(center.y, 5.0); + EXPECT_DOUBLE_EQ(center.z, 5.0); +} + +TEST_F(ReadAtomsHelperTest, CalculateLatticeCenterNonCubic) +{ + ModuleBase::Matrix3 latvec; + latvec.e11 = 8.0; latvec.e12 = 0.0; latvec.e13 = 0.0; + latvec.e21 = 2.0; latvec.e22 = 6.0; latvec.e23 = 0.0; + latvec.e31 = 1.0; latvec.e32 = 1.0; latvec.e33 = 10.0; + + auto center = unitcell::calculate_lattice_center(latvec, "xyz"); + + EXPECT_DOUBLE_EQ(center.x, (8.0 + 2.0 + 1.0) / 2.0); + EXPECT_DOUBLE_EQ(center.y, (0.0 + 6.0 + 1.0) / 2.0); + EXPECT_DOUBLE_EQ(center.z, (0.0 + 0.0 + 10.0) / 2.0); +} + +// Test allocate_atom_properties function +TEST_F(ReadAtomsHelperTest, AllocateAtomProperties) +{ + Atom atom; + int na = 5; + double mass = 12.0; + + unitcell::allocate_atom_properties(atom, na, mass); + + EXPECT_EQ(atom.tau.size(), na); + EXPECT_EQ(atom.dis.size(), na); + EXPECT_EQ(atom.taud.size(), na); + EXPECT_EQ(atom.boundary_shift.size(), na); + EXPECT_EQ(atom.vel.size(), na); + EXPECT_EQ(atom.mbl.size(), na); + EXPECT_EQ(atom.mag.size(), na); + EXPECT_EQ(atom.angle1.size(), na); + EXPECT_EQ(atom.angle2.size(), na); + EXPECT_EQ(atom.m_loc_.size(), na); + EXPECT_EQ(atom.lambda.size(), na); + EXPECT_EQ(atom.constrain.size(), na); + EXPECT_DOUBLE_EQ(atom.mass, mass); +} + +// Test transform_atom_coordinates for Direct coordinates +TEST_F(ReadAtomsHelperTest, TransformAtomCoordinatesDirect) +{ + Atom atom; + atom.tau.resize(1); + atom.taud.resize(1); + + ModuleBase::Vector3 v(0.5, 0.5, 0.5); + ModuleBase::Matrix3 latvec; + latvec.e11 = 10.0; latvec.e12 = 0.0; latvec.e13 = 0.0; + latvec.e21 = 0.0; latvec.e22 = 10.0; latvec.e23 = 0.0; + latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 10.0; + + double lat0 = 1.0; + ModuleBase::Vector3 latcenter; + + unitcell::transform_atom_coordinates(atom, 0, "Direct", v, latvec, lat0, latcenter); + + EXPECT_DOUBLE_EQ(atom.taud[0].x, 0.5); + EXPECT_DOUBLE_EQ(atom.taud[0].y, 0.5); + EXPECT_DOUBLE_EQ(atom.taud[0].z, 0.5); + EXPECT_DOUBLE_EQ(atom.tau[0].x, 5.0); + EXPECT_DOUBLE_EQ(atom.tau[0].y, 5.0); + EXPECT_DOUBLE_EQ(atom.tau[0].z, 5.0); +} + +// Test transform_atom_coordinates for Cartesian coordinates +TEST_F(ReadAtomsHelperTest, TransformAtomCoordinatesCartesian) +{ + Atom atom; + atom.tau.resize(1); + atom.taud.resize(1); + + ModuleBase::Vector3 v(5.0, 5.0, 5.0); + ModuleBase::Matrix3 latvec; + latvec.e11 = 10.0; latvec.e12 = 0.0; latvec.e13 = 0.0; + latvec.e21 = 0.0; latvec.e22 = 10.0; latvec.e23 = 0.0; + latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 10.0; + + double lat0 = 1.0; + ModuleBase::Vector3 latcenter; + + unitcell::transform_atom_coordinates(atom, 0, "Cartesian", v, latvec, lat0, latcenter); + + EXPECT_DOUBLE_EQ(atom.tau[0].x, 5.0); + EXPECT_DOUBLE_EQ(atom.tau[0].y, 5.0); + EXPECT_DOUBLE_EQ(atom.tau[0].z, 5.0); + EXPECT_DOUBLE_EQ(atom.taud[0].x, 0.5); + EXPECT_DOUBLE_EQ(atom.taud[0].y, 0.5); + EXPECT_DOUBLE_EQ(atom.taud[0].z, 0.5); +} + +// Test process_magnetization for nspin=2 +TEST_F(ReadAtomsHelperTest, ProcessMagnetizationNspin2) +{ + Atom atom; + atom.mag.resize(1); + atom.m_loc_.resize(1); + atom.angle1.resize(1); + atom.angle2.resize(1); + + atom.mag[0] = 2.0; + atom.m_loc_[0].set(0, 0, 0); + + unitcell::process_magnetization(atom, 0, 0, 2, false, false, ofs_running); + + // For nspin=2, only z component should be set + EXPECT_DOUBLE_EQ(atom.m_loc_[0].x, 0.0); + EXPECT_DOUBLE_EQ(atom.m_loc_[0].y, 0.0); + EXPECT_DOUBLE_EQ(atom.m_loc_[0].z, 2.0); + EXPECT_DOUBLE_EQ(atom.mag[0], 2.0); +} + +// Test process_magnetization for nspin=4 with vector input +TEST_F(ReadAtomsHelperTest, ProcessMagnetizationNspin4VectorInput) +{ + Atom atom; + atom.mag.resize(1); + atom.m_loc_.resize(1); + atom.angle1.resize(1); + atom.angle2.resize(1); + + atom.m_loc_[0].set(1.0, 1.0, 1.0); + atom.mag[0] = sqrt(3.0); + + // Set noncolin to true to allow non-collinear magnetization + // Note: This requires PARAM to be properly initialized + + unitcell::process_magnetization(atom, 0, 0, 4, true, false, ofs_running); + + // Angles should be calculated from vector components + EXPECT_GT(atom.angle1[0], 0.0); + EXPECT_GT(atom.angle2[0], 0.0); +} + +// Test process_magnetization with angle input +TEST_F(ReadAtomsHelperTest, ProcessMagnetizationAngleInput) +{ + Atom atom; + atom.mag.resize(1); + atom.m_loc_.resize(1); + atom.angle1.resize(1); + atom.angle2.resize(1); + + atom.mag[0] = 2.0; + atom.angle1[0] = M_PI / 2.0; // 90 degrees + atom.angle2[0] = 0.0; + atom.m_loc_[0].set(0, 0, 0); + + // Note: For nspin=4, if noncolin is false (default), x and y components are zeroed + // So we test with nspin=2 instead to verify the angle calculation works + unitcell::process_magnetization(atom, 0, 0, 2, false, true, ofs_running); + + // For nspin=2, only z component is used, which should be mag[0] * cos(angle1) + // With angle1 = PI/2, cos(PI/2) = 0 + EXPECT_NEAR(atom.m_loc_[0].z, 0.0, 1e-10); + EXPECT_DOUBLE_EQ(atom.mag[0], atom.m_loc_[0].z); +} + +// Test parse_atom_properties with movement flags +TEST_F(ReadAtomsHelperTest, ParseAtomPropertiesMovementFlags) +{ + std::string input_str = "1.0 2.0 3.0 m 1 0 1\n"; + std::istringstream iss(input_str); + + // Create a temporary file for testing + std::ofstream temp_file("test_input.tmp"); + temp_file << input_str; + temp_file.close(); + + std::ifstream ifpos("test_input.tmp"); + + Atom atom; + atom.label = "C"; + atom.vel.resize(1); + atom.mag.resize(1); + atom.m_loc_.resize(1); + atom.angle1.resize(1); + atom.angle2.resize(1); + atom.lambda.resize(1); + atom.constrain.resize(1); + + ModuleBase::Vector3 mv(1, 1, 1); + bool input_vec_mag = false; + bool input_angle_mag = false; + bool set_element_mag_zero = false; + + // Skip the position coordinates + double x, y, z; + ifpos >> x >> y >> z; + + bool result = unitcell::parse_atom_properties(ifpos, atom, 0, mv, + input_vec_mag, input_angle_mag, + set_element_mag_zero); + + EXPECT_TRUE(result); + EXPECT_EQ(mv.x, 1); + EXPECT_EQ(mv.y, 0); + EXPECT_EQ(mv.z, 1); + + ifpos.close(); + std::remove("test_input.tmp"); +} + +// Test parse_atom_properties with velocity +TEST_F(ReadAtomsHelperTest, ParseAtomPropertiesVelocity) +{ + std::string input_str = "1.0 2.0 3.0 v 0.1 0.2 0.3\n"; + + std::ofstream temp_file("test_input.tmp"); + temp_file << input_str; + temp_file.close(); + + std::ifstream ifpos("test_input.tmp"); + + Atom atom; + atom.label = "C"; + atom.vel.resize(1); + atom.mag.resize(1); + atom.m_loc_.resize(1); + atom.angle1.resize(1); + atom.angle2.resize(1); + atom.lambda.resize(1); + atom.constrain.resize(1); + + ModuleBase::Vector3 mv(1, 1, 1); + bool input_vec_mag = false; + bool input_angle_mag = false; + bool set_element_mag_zero = false; + + // Skip the position coordinates + double x, y, z; + ifpos >> x >> y >> z; + + bool result = unitcell::parse_atom_properties(ifpos, atom, 0, mv, + input_vec_mag, input_angle_mag, + set_element_mag_zero); + + EXPECT_TRUE(result); + EXPECT_DOUBLE_EQ(atom.vel[0].x, 0.1); + EXPECT_DOUBLE_EQ(atom.vel[0].y, 0.2); + EXPECT_DOUBLE_EQ(atom.vel[0].z, 0.3); + + ifpos.close(); + std::remove("test_input.tmp"); +} + +// Test parse_atom_properties with scalar magnetization +TEST_F(ReadAtomsHelperTest, ParseAtomPropertiesScalarMag) +{ + std::string input_str = "1.0 2.0 3.0 mag 2.5\n"; + + std::ofstream temp_file("test_input.tmp"); + temp_file << input_str; + temp_file.close(); + + std::ifstream ifpos("test_input.tmp"); + + Atom atom; + atom.label = "C"; + atom.vel.resize(1); + atom.mag.resize(1); + atom.m_loc_.resize(1); + atom.angle1.resize(1); + atom.angle2.resize(1); + atom.lambda.resize(1); + atom.constrain.resize(1); + + ModuleBase::Vector3 mv(1, 1, 1); + bool input_vec_mag = false; + bool input_angle_mag = false; + bool set_element_mag_zero = false; + + // Skip the position coordinates + double x, y, z; + ifpos >> x >> y >> z; + + bool result = unitcell::parse_atom_properties(ifpos, atom, 0, mv, + input_vec_mag, input_angle_mag, + set_element_mag_zero); + + EXPECT_TRUE(result); + EXPECT_DOUBLE_EQ(atom.mag[0], 2.5); + EXPECT_TRUE(set_element_mag_zero); + EXPECT_FALSE(input_vec_mag); + + ifpos.close(); + std::remove("test_input.tmp"); +} + +// Test parse_atom_properties with vector magnetization +TEST_F(ReadAtomsHelperTest, ParseAtomPropertiesVectorMag) +{ + std::string input_str = "1.0 2.0 3.0 mag 1.0 2.0 3.0\n"; + + std::ofstream temp_file("test_input.tmp"); + temp_file << input_str; + temp_file.close(); + + std::ifstream ifpos("test_input.tmp"); + + Atom atom; + atom.label = "C"; + atom.vel.resize(1); + atom.mag.resize(1); + atom.m_loc_.resize(1); + atom.angle1.resize(1); + atom.angle2.resize(1); + atom.lambda.resize(1); + atom.constrain.resize(1); + + ModuleBase::Vector3 mv(1, 1, 1); + bool input_vec_mag = false; + bool input_angle_mag = false; + bool set_element_mag_zero = false; + + // Skip the position coordinates + double x, y, z; + ifpos >> x >> y >> z; + + bool result = unitcell::parse_atom_properties(ifpos, atom, 0, mv, + input_vec_mag, input_angle_mag, + set_element_mag_zero); + + EXPECT_TRUE(result); + EXPECT_DOUBLE_EQ(atom.m_loc_[0].x, 1.0); + EXPECT_DOUBLE_EQ(atom.m_loc_[0].y, 2.0); + EXPECT_DOUBLE_EQ(atom.m_loc_[0].z, 3.0); + EXPECT_NEAR(atom.mag[0], sqrt(1.0 + 4.0 + 9.0), 1e-10); + EXPECT_TRUE(input_vec_mag); + EXPECT_TRUE(set_element_mag_zero); + + ifpos.close(); + std::remove("test_input.tmp"); +} + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/source/source_cell/test_pw/CMakeLists.txt b/source/source_cell/test_pw/CMakeLists.txt index 48893a14ae..e05ecbc4cc 100644 --- a/source/source_cell/test_pw/CMakeLists.txt +++ b/source/source_cell/test_pw/CMakeLists.txt @@ -10,7 +10,7 @@ install(FILES unitcell_test_pw_para.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) AddTest( TARGET MODULE_CELL_unitcell_test_pw LIBS parameter ${math_libs} base device - SOURCES unitcell_test_pw.cpp ../unitcell.cpp ../read_atoms.cpp ../atom_spec.cpp ../update_cell.cpp ../bcast_cell.cpp + SOURCES unitcell_test_pw.cpp ../unitcell.cpp ../read_atoms.cpp ../read_atoms_helper.cpp ../atom_spec.cpp ../update_cell.cpp ../bcast_cell.cpp ../atom_pseudo.cpp ../pseudo.cpp ../read_pp.cpp ../read_pp_complete.cpp ../read_pp_upf201.cpp ../read_pp_upf100.cpp ../read_stru.cpp ../read_atom_species.cpp ../read_pp_vwr.cpp ../read_pp_blps.cpp ../../source_io/output.cpp diff --git a/source/source_lcao/module_deepks/test/CMakeLists.txt b/source/source_lcao/module_deepks/test/CMakeLists.txt index 36327d3f1d..fae7adc5b8 100644 --- a/source/source_lcao/module_deepks/test/CMakeLists.txt +++ b/source/source_lcao/module_deepks/test/CMakeLists.txt @@ -7,6 +7,7 @@ add_executable( ../../../source_cell/atom_spec.cpp ../../../source_cell/atom_pseudo.cpp ../../../source_cell/read_atoms.cpp + ../../../source_cell/read_atoms_helper.cpp ../../../source_cell/read_stru.cpp ../../../source_cell/print_cell.cpp ../../../source_cell/read_atom_species.cpp diff --git a/source/source_lcao/module_deepks/test/Makefile.Objects b/source/source_lcao/module_deepks/test/Makefile.Objects index 42b0143270..31b32cef5a 100644 --- a/source/source_lcao/module_deepks/test/Makefile.Objects +++ b/source/source_lcao/module_deepks/test/Makefile.Objects @@ -57,6 +57,7 @@ read_pp_blps.o\ unitcell.o\ check_atomic_stru.o\ read_atoms.o\ +read_atoms_helper.o\ read_cell_pseudopots.o\ setup_nonlocal.o diff --git a/source/source_md/test/CMakeLists.txt b/source/source_md/test/CMakeLists.txt index 516fc52b0c..aadccfd449 100644 --- a/source/source_md/test/CMakeLists.txt +++ b/source/source_md/test/CMakeLists.txt @@ -12,6 +12,7 @@ list(APPEND depend_files ../../source_cell/read_atom_species.cpp ../../source_cell/atom_pseudo.cpp ../../source_cell/read_atoms.cpp + ../../source_cell/read_atoms_helper.cpp ../../source_cell/pseudo.cpp ../../source_cell/read_pp.cpp ../../source_cell/read_pp_complete.cpp diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index 5922b02567..8a05dfdf2a 100644 --- a/source/source_pw/module_pwdft/test/CMakeLists.txt +++ b/source/source_pw/module_pwdft/test/CMakeLists.txt @@ -43,6 +43,7 @@ AddTest( ../../../source_cell/read_stru.cpp ../../../source_cell/read_atom_species.cpp ../../../source_cell/read_atoms.cpp + ../../../source_cell/read_atoms_helper.cpp ../../../source_cell/read_pp.cpp ../../../source_cell/read_pp_complete.cpp ../../../source_cell/read_pp_upf100.cpp