|
| 1 | +#ifndef PREPARE_UNITCELL_H |
| 2 | +#define PREPARE_UNITCELL_H |
| 3 | + |
| 4 | +class UcellTestPrepare |
| 5 | +{ |
| 6 | +public: |
| 7 | + UcellTestPrepare(std::string system_name_in, |
| 8 | + std::string latname_in, |
| 9 | + int lmaxmax_in, |
| 10 | + bool init_vel_in, |
| 11 | + bool selective_dynamics_in, |
| 12 | + bool relax_new_in, |
| 13 | + std::string fixed_axes_in, |
| 14 | + double lat0_in, |
| 15 | + std::valarray<double> latvec_in, |
| 16 | + std::vector<std::string> elements_in, |
| 17 | + std::vector<std::string> pp_files_in, |
| 18 | + std::vector<std::string> pp_types_in, |
| 19 | + std::vector<std::string> orb_files_in, |
| 20 | + std::valarray<int> natom_in, |
| 21 | + std::vector<double> atomic_mass_in, |
| 22 | + std::string coor_type_in, |
| 23 | + std::valarray<double> coordinates_in, |
| 24 | + std::valarray<double> mbl_in = {0}, |
| 25 | + std::valarray<double> velocity_in = {0} |
| 26 | + ): |
| 27 | + system_name(system_name_in), |
| 28 | + latname(latname_in), |
| 29 | + lmaxmax(lmaxmax_in), |
| 30 | + init_vel(init_vel_in), |
| 31 | + selective_dynamics(selective_dynamics_in), |
| 32 | + relax_new(relax_new_in), |
| 33 | + fixed_axes(fixed_axes_in), |
| 34 | + lat0(lat0_in), |
| 35 | + latvec(latvec_in), |
| 36 | + elements(elements_in), |
| 37 | + pp_files(pp_files_in), |
| 38 | + pp_types(pp_types_in), |
| 39 | + orb_files(orb_files_in), |
| 40 | + natom(natom_in), |
| 41 | + atomic_mass(atomic_mass_in), |
| 42 | + coor_type(coor_type_in), |
| 43 | + coordinates(coordinates_in), |
| 44 | + mbl(mbl_in), |
| 45 | + velocity(velocity_in) // velocity assume the existence of mbl in print_stru_file() |
| 46 | + {} |
| 47 | + std::string system_name; |
| 48 | + std::string latname; |
| 49 | + int lmaxmax; |
| 50 | + bool init_vel; |
| 51 | + bool selective_dynamics; |
| 52 | + bool relax_new; |
| 53 | + std::string fixed_axes; |
| 54 | + double lat0; |
| 55 | + std::valarray<double> latvec; |
| 56 | + std::vector<std::string> elements; |
| 57 | + std::vector<std::string> pp_files; |
| 58 | + std::vector<std::string> pp_types; |
| 59 | + std::vector<std::string> orb_files; |
| 60 | + std::valarray<int> natom; |
| 61 | + std::vector<double> atomic_mass; |
| 62 | + std::string coor_type; |
| 63 | + std::valarray<double> coordinates; |
| 64 | + std::valarray<double> mbl; |
| 65 | + std::valarray<double> velocity; |
| 66 | + // ntype |
| 67 | + int ntype; |
| 68 | + int atomic_index; |
| 69 | + |
| 70 | + std::unique_ptr<UnitCell> SetUcellInfo() |
| 71 | + { |
| 72 | + //basic info |
| 73 | + this->ntype = this->elements.size(); |
| 74 | + std::unique_ptr<UnitCell> ucell(new UnitCell); |
| 75 | + ucell->setup(this->latname, |
| 76 | + this->ntype, |
| 77 | + this->lmaxmax, |
| 78 | + this->init_vel, |
| 79 | + this->fixed_axes); |
| 80 | + delete[] ucell->atom_label; |
| 81 | + delete[] ucell->atom_mass; |
| 82 | + delete[] ucell->pseudo_fn; |
| 83 | + delete[] ucell->pseudo_type; |
| 84 | + delete[] ucell->orbital_fn; |
| 85 | + delete[] ucell->magnet.start_magnetization; //mag set here |
| 86 | + ucell->atom_label = new std::string[ucell->ntype]; |
| 87 | + ucell->atom_mass = new double[ucell->ntype]; |
| 88 | + ucell->pseudo_fn = new std::string[ucell->ntype]; |
| 89 | + ucell->pseudo_type = new std::string[ucell->ntype]; |
| 90 | + ucell->orbital_fn = new std::string[ucell->ntype]; |
| 91 | + ucell->magnet.start_magnetization = new double[ucell->ntype]; //mag set here |
| 92 | + ucell->magnet.ux_[0] = 0.0; // ux_ set here |
| 93 | + ucell->magnet.ux_[1] = 0.0; |
| 94 | + ucell->magnet.ux_[2] = 0.0; |
| 95 | + for(int it=0;it<ucell->ntype;++it) |
| 96 | + { |
| 97 | + ucell->atom_label[it] = this->elements[it]; |
| 98 | + ucell->atom_mass[it] = this->atomic_mass[it]; |
| 99 | + ucell->pseudo_fn[it] = this->pp_files[it]; |
| 100 | + ucell->pseudo_type[it] = this->pp_types[it]; |
| 101 | + ucell->orbital_fn[it] = this->orb_files[it]; |
| 102 | + ucell->magnet.start_magnetization[it] = 0.0; //mag set here |
| 103 | + } |
| 104 | + //lattice info |
| 105 | + ucell->lat0 = this->lat0; |
| 106 | + ucell->lat0_angstrom = ucell->lat0 * 0.529177; |
| 107 | + ucell->tpiba = ModuleBase::TWO_PI/ucell->lat0; |
| 108 | + ucell->tpiba2 = ucell->tpiba * ucell->tpiba; |
| 109 | + ucell->latvec.e11 = this->latvec[0]; |
| 110 | + ucell->latvec.e12 = this->latvec[1]; |
| 111 | + ucell->latvec.e13 = this->latvec[2]; |
| 112 | + ucell->latvec.e21 = this->latvec[3]; |
| 113 | + ucell->latvec.e22 = this->latvec[4]; |
| 114 | + ucell->latvec.e23 = this->latvec[5]; |
| 115 | + ucell->latvec.e31 = this->latvec[6]; |
| 116 | + ucell->latvec.e32 = this->latvec[7]; |
| 117 | + ucell->latvec.e33 = this->latvec[8]; |
| 118 | + ucell->a1.x = ucell->latvec.e11; |
| 119 | + ucell->a1.y = ucell->latvec.e12; |
| 120 | + ucell->a1.z = ucell->latvec.e13; |
| 121 | + ucell->a2.x = ucell->latvec.e21; |
| 122 | + ucell->a2.y = ucell->latvec.e22; |
| 123 | + ucell->a2.z = ucell->latvec.e23; |
| 124 | + ucell->a3.x = ucell->latvec.e31; |
| 125 | + ucell->a3.y = ucell->latvec.e32; |
| 126 | + ucell->a3.z = ucell->latvec.e33; |
| 127 | + ucell->GT = ucell->latvec.Inverse(); |
| 128 | + ucell->G = ucell->GT.Transpose(); |
| 129 | + ucell->GGT = ucell->G*ucell->GT; |
| 130 | + ucell->invGGT = ucell->GGT.Inverse(); |
| 131 | + ucell->omega = abs(ucell->latvec.Det())*(ucell->lat0)*(ucell->lat0)*(ucell->lat0); |
| 132 | + //atomic info |
| 133 | + ucell->Coordinate = this->coor_type; |
| 134 | + ucell->atoms = new Atom[ucell->ntype]; |
| 135 | + ucell->set_atom_flag = true; |
| 136 | + this->atomic_index = 0; |
| 137 | + for(int it=0;it<ucell->ntype;++it) |
| 138 | + { |
| 139 | + ucell->atoms[it].label = this->elements[it]; |
| 140 | + ucell->atoms[it].nw = 0; |
| 141 | + ucell->atoms[it].nwl = 2; |
| 142 | + delete[] ucell->atoms[it].l_nchi; |
| 143 | + ucell->atoms[it].l_nchi = new int[ ucell->atoms[it].nwl+1]; |
| 144 | + for(int L=0; L<ucell->atoms[it].nwl+1; L++) |
| 145 | + { |
| 146 | + ucell->atoms[it].l_nchi[L] = 1; |
| 147 | + ucell->atoms[it].nw += (2*L + 1) * ucell->atoms[it].l_nchi[L]; |
| 148 | + } |
| 149 | + ucell->atoms[it].na = this->natom[it]; |
| 150 | + //coordinates and related physical quantities |
| 151 | + delete[] ucell->atoms[it].tau; |
| 152 | + delete[] ucell->atoms[it].tau_original; |
| 153 | + delete[] ucell->atoms[it].taud; |
| 154 | + delete[] ucell->atoms[it].vel; |
| 155 | + delete[] ucell->atoms[it].mag; |
| 156 | + delete[] ucell->atoms[it].angle1; |
| 157 | + delete[] ucell->atoms[it].angle2; |
| 158 | + delete[] ucell->atoms[it].m_loc_; |
| 159 | + delete[] ucell->atoms[it].mbl; |
| 160 | + ucell->atoms[it].tau = new ModuleBase::Vector3<double>[ucell->atoms[it].na]; |
| 161 | + ucell->atoms[it].tau_original = new ModuleBase::Vector3<double>[ucell->atoms[it].na]; |
| 162 | + ucell->atoms[it].taud = new ModuleBase::Vector3<double>[ucell->atoms[it].na]; |
| 163 | + ucell->atoms[it].vel = new ModuleBase::Vector3<double>[ucell->atoms[it].na]; |
| 164 | + ucell->atoms[it].mag = new double[ucell->atoms[it].na]; |
| 165 | + ucell->atoms[it].angle1 = new double[ucell->atoms[it].na]; |
| 166 | + ucell->atoms[it].angle2 = new double[ucell->atoms[it].na]; |
| 167 | + ucell->atoms[it].m_loc_ = new ModuleBase::Vector3<double>[ucell->atoms[it].na]; |
| 168 | + ucell->atoms[it].mbl = new ModuleBase::Vector3<int>[ucell->atoms[it].na]; |
| 169 | + ucell->atoms[it].mass = ucell->atom_mass[it]; // mass set here |
| 170 | + for(int ia=0; ia<ucell->atoms[it].na; ++ia) |
| 171 | + { |
| 172 | + if (ucell->Coordinate == "Direct") |
| 173 | + { |
| 174 | + ucell->atoms[it].taud[ia].x = this->coordinates[this->atomic_index*3+0]; |
| 175 | + ucell->atoms[it].taud[ia].y = this->coordinates[this->atomic_index*3+1]; |
| 176 | + ucell->atoms[it].taud[ia].z = this->coordinates[this->atomic_index*3+2]; |
| 177 | + ucell->atoms[it].tau[ia] = ucell->atoms[it].taud[ia]*ucell->latvec; |
| 178 | + } |
| 179 | + else if (ucell->Coordinate == "Cartesian") |
| 180 | + { |
| 181 | + ucell->atoms[it].tau[ia].x = this->coordinates[this->atomic_index*3+0]; |
| 182 | + ucell->atoms[it].tau[ia].y = this->coordinates[this->atomic_index*3+1]; |
| 183 | + ucell->atoms[it].tau[ia].z = this->coordinates[this->atomic_index*3+2]; |
| 184 | + ModuleBase::Mathzone::Cartesian_to_Direct( |
| 185 | + ucell->atoms[it].tau[ia].x, ucell->atoms[it].tau[ia].y, ucell->atoms[it].tau[ia].z, |
| 186 | + ucell->latvec.e11, ucell->latvec.e12, ucell->latvec.e13, |
| 187 | + ucell->latvec.e21, ucell->latvec.e22, ucell->latvec.e23, |
| 188 | + ucell->latvec.e31, ucell->latvec.e32, ucell->latvec.e33, |
| 189 | + ucell->atoms[it].taud[ia].x, ucell->atoms[it].taud[ia].y, ucell->atoms[it].taud[ia].z); |
| 190 | + } |
| 191 | + ucell->atoms[it].tau_original[ia] = ucell->atoms[it].tau[ia]; |
| 192 | + if(this->init_vel) |
| 193 | + { |
| 194 | + ucell->atoms[it].vel[ia].x = this->velocity[this->atomic_index*3+0]; |
| 195 | + ucell->atoms[it].vel[ia].y = this->velocity[this->atomic_index*3+1]; |
| 196 | + ucell->atoms[it].vel[ia].z = this->velocity[this->atomic_index*3+2]; |
| 197 | + } |
| 198 | + else |
| 199 | + { |
| 200 | + ucell->atoms[it].vel[ia].set(0,0,0); |
| 201 | + } |
| 202 | + ucell->atoms[it].m_loc_[ia].set(0,0,0); |
| 203 | + ucell->atoms[it].angle1[ia] = 0; |
| 204 | + ucell->atoms[it].angle2[ia] = 0; |
| 205 | + if(this->selective_dynamics) |
| 206 | + { |
| 207 | + ucell->atoms[it].mbl[ia].x = this->mbl[this->atomic_index*3+0]; |
| 208 | + ucell->atoms[it].mbl[ia].y = this->mbl[this->atomic_index*3+1]; |
| 209 | + ucell->atoms[it].mbl[ia].z = this->mbl[this->atomic_index*3+2]; |
| 210 | + } |
| 211 | + else |
| 212 | + { |
| 213 | + ucell->atoms[it].mbl[ia] = {1,1,1}; |
| 214 | + } |
| 215 | + ++(this->atomic_index); |
| 216 | + } |
| 217 | + } |
| 218 | + ucell->nat = this->natom.sum(); |
| 219 | + return ucell; |
| 220 | + } |
| 221 | +}; |
| 222 | + |
| 223 | +#endif |
0 commit comments