Skip to content

Commit

Permalink
Refactor: rescale vel to temperature when read in vel
Browse files Browse the repository at this point in the history
  • Loading branch information
@YuLiu98
YuLiu98 committed Aug 12, 2024
1 parent 1ddf8e2 commit 9c512d4
Show file tree
Hide file tree
Showing 3 changed files with 91 additions and 37 deletions.
61 changes: 36 additions & 25 deletions source/module_md/md_func.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,15 +83,18 @@ void read_vel(const UnitCell& unit_in, ModuleBase::Vector3<double>* vel)
for (int ia = 0; ia < unit_in.atoms[it].na; ++ia)
{
vel[iat] = unit_in.atoms[it].vel[ia];
if (unit_in.atoms[it].mbl[ia].x == 0) {
if (unit_in.atoms[it].mbl[ia].x == 0)
{
vel[iat].x = 0;
}
if (unit_in.atoms[it].mbl[ia].y == 0) {
}
if (unit_in.atoms[it].mbl[ia].y == 0)
{
vel[iat].y = 0;
}
if (unit_in.atoms[it].mbl[ia].z == 0) {
}
if (unit_in.atoms[it].mbl[ia].z == 0)
{
vel[iat].z = 0;
}
}
++iat;
}
}
Expand All @@ -100,6 +103,28 @@ void read_vel(const UnitCell& unit_in, ModuleBase::Vector3<double>* vel)
return;
}

void rescale_vel(const int& natom,
const double& temperature,
const double* allmass,
const int& frozen_freedom,
ModuleBase::Vector3<double>* vel)
{
double factor = 0.0;
if (3 * natom == frozen_freedom || temperature == 0)
{
factor = 0;
}
else
{
factor = 0.5 * (3 * natom - frozen_freedom) * temperature / kinetic_energy(natom, vel, allmass);
}

for (int i = 0; i < natom; i++)
{
vel[i] = vel[i] * sqrt(factor);
}
}

void rand_vel(const int& natom,
const double& temperature,
const double* allmass,
Expand Down Expand Up @@ -146,20 +171,8 @@ void rand_vel(const int& natom,
}
}

double factor=0.0;
if (3 * natom == frozen_freedom || temperature == 0)
{
factor = 0;
}
else
{
factor = 0.5 * (3 * natom - frozen_freedom) * temperature / kinetic_energy(natom, vel, allmass);
}

for (int i = 0; i < natom; i++)
{
vel[i] = vel[i] * sqrt(factor);
}
// rescale the velocity to the target temperature
rescale_vel(natom, temperature, allmass, frozen_freedom, vel);
}

#ifdef __MPI
Expand Down Expand Up @@ -212,16 +225,14 @@ void init_vel(const UnitCell& unit_in,
<< std::endl;
temperature = t_current;
}
else if (std::fabs(temperature - t_current) > 1e-8)
else
{
std::cout << " INITIAL TEMPERATURE IN INPUT = " << temperature * ModuleBase::Hartree_to_K << " K"
<< std::endl;
std::cout << " READING TEMPERATURE FROM STRU = " << t_current * ModuleBase::Hartree_to_K << " K"
<< std::endl;
std::cout << " INCONSISTENCE, PLEASE CHECK" << std::endl;
std::cout << " ------------------------------------- DONE -----------------------------------------"
<< std::endl;
exit(0);
std::cout << " RESCALE VEL TO INITIAL TEMPERATURE" << std::endl;
rescale_vel(unit_in.nat, temperature, allmass, frozen_freedom, vel);
}
}
else
Expand Down
15 changes: 15 additions & 0 deletions source/module_md/md_func.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,21 @@ void rand_vel(const int& natom,
const int& my_rank,
ModuleBase::Vector3<double>* vel);

/**
* @brief rescale the velocity to the target temperature
*
* @param natom the number of atoms
* @param temperature ion temperature
* @param allmass atomic mass
* @param frozen_freedom the fixed freedom
* @param vel the genarated atomic velocities
*/
void rescale_vel(const int& natom,
const double& temperature,
const double* allmass,
const int& frozen_freedom,
ModuleBase::Vector3<double>* vel);

/**
* @brief calculate energy, forces and virial tensor
*
Expand Down
52 changes: 40 additions & 12 deletions source/module_md/test/md_func_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,9 @@
* - MD_func::init_vel
* - initialize the atomic velocities
*
* - MD_func::rescale_vel
* - rescale the velocity to the target temperature
*
* - MD_func::compute_stress
* - calculate the contribution of classical kinetic energy of atoms to stress
*
Expand Down Expand Up @@ -147,6 +150,33 @@ TEST_F(MD_func_test, readvel)
EXPECT_DOUBLE_EQ(vel[3].z, -2.83313122596e-05);
}

TEST_F(MD_func_test, RescaleVel)
{
int frozen_freedom = 3;
for (int i = 0; i < natom; ++i)
{
allmass[i] = 39.948 / ModuleBase::AU_to_MASS;
vel[i].x = 0.1;
vel[i].y = 0.2;
vel[i].z = 0.3;
}
temperature = 300.0 / ModuleBase::Hartree_to_K;
MD_func::rescale_vel(natom, temperature, allmass, frozen_freedom, vel);

EXPECT_DOUBLE_EQ(vel[0].x, 4.579010775069125e-05);
EXPECT_DOUBLE_EQ(vel[0].y, 9.15802155013825e-05);
EXPECT_DOUBLE_EQ(vel[0].z, 0.00013737032325207373);
EXPECT_DOUBLE_EQ(vel[1].x, 4.579010775069125e-05);
EXPECT_DOUBLE_EQ(vel[1].y, 9.15802155013825e-05);
EXPECT_DOUBLE_EQ(vel[1].z, 0.00013737032325207373);
EXPECT_DOUBLE_EQ(vel[2].x, 4.579010775069125e-05);
EXPECT_DOUBLE_EQ(vel[2].y, 9.15802155013825e-05);
EXPECT_DOUBLE_EQ(vel[2].z, 0.00013737032325207373);
EXPECT_DOUBLE_EQ(vel[3].x, 4.579010775069125e-05);
EXPECT_DOUBLE_EQ(vel[3].y, 9.15802155013825e-05);
EXPECT_DOUBLE_EQ(vel[3].z, 0.00013737032325207373);
}

TEST_F(MD_func_test, InitVelCase1)
{
ucell.init_vel = 1;
Expand All @@ -170,9 +200,7 @@ TEST_F(MD_func_test, InitVelCase3)
ucell.init_vel = 1;
temperature = 310.0 / ModuleBase::Hartree_to_K;

EXPECT_EXIT(MD_func::init_vel(ucell, GlobalV::MY_RANK, false, temperature, allmass, frozen_freedom, ionmbl, vel),
::testing::ExitedWithCode(0),
"");
EXPECT_DOUBLE_EQ(temperature, 310.0 / ModuleBase::Hartree_to_K);
}

TEST_F(MD_func_test, InitVelCase4)
Expand All @@ -198,15 +226,15 @@ TEST_F(MD_func_test, compute_stress)
temperature = 300.0 / ModuleBase::Hartree_to_K;
MD_func::init_vel(ucell, GlobalV::MY_RANK, false, temperature, allmass, frozen_freedom, ionmbl, vel);
MD_func::compute_stress(ucell, vel, allmass, true, virial, stress);
EXPECT_DOUBLE_EQ(stress(0, 0), 5.2064533063673623e-06);
EXPECT_DOUBLE_EQ(stress(0, 1), -1.6467487572445481e-06);
EXPECT_DOUBLE_EQ(stress(0, 2), 1.5039983732220751e-06);
EXPECT_DOUBLE_EQ(stress(1, 0), -1.6467487572445481e-06);
EXPECT_DOUBLE_EQ(stress(1, 1), 2.3806464376131247e-06);
EXPECT_DOUBLE_EQ(stress(1, 2), -1.251414906590483e-06);
EXPECT_DOUBLE_EQ(stress(2, 0), 1.5039983732220751e-06);
EXPECT_DOUBLE_EQ(stress(2, 1), -1.251414906590483e-06);
EXPECT_DOUBLE_EQ(stress(2, 2), 9.6330189688582584e-07);
EXPECT_DOUBLE_EQ(stress(0, 0), 5.2064533063674207e-06);
EXPECT_DOUBLE_EQ(stress(0, 1), -1.6467487572445666e-06);
EXPECT_DOUBLE_EQ(stress(0, 2), 1.5039983732220917e-06);
EXPECT_DOUBLE_EQ(stress(1, 0), -1.6467487572445661e-06);
EXPECT_DOUBLE_EQ(stress(1, 1), 2.380646437613151e-06);
EXPECT_DOUBLE_EQ(stress(1, 2), -1.2514149065904968e-06);
EXPECT_DOUBLE_EQ(stress(2, 0), 1.5039983732220917e-06);
EXPECT_DOUBLE_EQ(stress(2, 1), -1.2514149065904968e-06);
EXPECT_DOUBLE_EQ(stress(2, 2), 9.6330189688583664e-07);
}

TEST_F(MD_func_test, dump_info)
Expand Down

0 comments on commit 9c512d4

Please sign in to comment.