Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Refactor: rescale vel to temperature when read in vel #4953

Merged
merged 2 commits into from
Aug 12, 2024
Merged

Refactor: rescale vel to temperature when read in vel #4953

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9c512d4
Refactor: rescale vel to temperature when read in vel
YuLiu98 Aug 12, 2024
9ce93b5
Tests: update unittests
YuLiu98 Aug 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 36 additions & 25 deletions source/module_md/md_func.cpp
View file
Open in desktop
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
Open in desktop
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
2 changes: 1 addition & 1 deletion source/module_md/test/fire_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ class FIREtest : public testing::Test

TEST_F(FIREtest, Setup)
{
EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999665, doublethreshold);
EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999994, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 0), 6.0100555286436806e-06, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 1), -1.4746713013791574e-06, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 2), 1.5039983732220751e-06, doublethreshold);
Expand Down
2 changes: 1 addition & 1 deletion source/module_md/test/langevin_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ class Langevin_test : public testing::Test

TEST_F(Langevin_test, setup)
{
EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999665, doublethreshold);
EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999994, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 0), 6.0100555286436806e-06, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 1), -1.4746713013791574e-06, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 2), 1.5039983732220751e-06, doublethreshold);
Expand Down
52 changes: 40 additions & 12 deletions source/module_md/test/md_func_test.cpp
View file
Open in desktop
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
2 changes: 1 addition & 1 deletion source/module_md/test/nhchain_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ class NHC_test : public testing::Test

TEST_F(NHC_test, setup)
{
EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999665, doublethreshold);
EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999994, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 0), 6.0100555286436806e-06, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 1), -1.4746713013791574e-06, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 2), 1.5039983732220751e-06, doublethreshold);
Expand Down
2 changes: 1 addition & 1 deletion source/module_md/test/verlet_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ class Verlet_test : public testing::Test

TEST_F(Verlet_test, setup)
{
EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999665, doublethreshold);
EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999994, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 0), 6.0100555286436806e-06, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 1), -1.4746713013791574e-06, doublethreshold);
EXPECT_NEAR(mdrun->stress(0, 2), 1.5039983732220751e-06, doublethreshold);
Expand Down
Loading