Remove legacy UnitCell velocity initialization

Author: Fei Yang

source/source_md/md_func.cpp

@@ -227,67 +227,6 @@ void rand_vel(const int& natom,
     return;
 }
 
-void init_vel(const UnitCell& unit_in,
-              const int& my_rank,
-              const bool& restart,
-              double& temperature,
-              double* allmass,
-              int& frozen_freedom,
-              ModuleBase::Vector3<int>* ionmbl,
-              ModuleBase::Vector3<double>* vel)
-{
-    ModuleBase::Vector3<int> frozen;
-    get_mass_mbl(unit_in, allmass, frozen, ionmbl);
-    frozen_freedom = frozen.x + frozen.y + frozen.z;
-    if (frozen.x == 0)
-    {
-        ++frozen_freedom;
-    }
-    if (frozen.y == 0)
-    {
-        ++frozen_freedom;
-    }
-    if (frozen.z == 0)
-    {
-        ++frozen_freedom;
-    }
-
-    if (unit_in.init_vel)
-    {
-        std::cout << " Reading velocities from STRU file" << std::endl;
-        read_vel(unit_in, vel);
-        double kinetic = 0.0;
-        double t_current = MD_func::current_temp(kinetic, unit_in.nat, frozen_freedom, allmass, vel);
-        if (restart)
-        {
-            std::cout << " Restart MD, current temperature is " << t_current * ModuleBase::Hartree_to_K << " K"
-                      << std::endl;
-        }
-        else if (temperature < 0)
-        {
-            std::cout << " Autoset the initial tempearture to " << t_current * ModuleBase::Hartree_to_K << " K"
-                      << std::endl;
-            temperature = t_current;
-        }
-        else
-        {
-            std::cout << " Initial temeprature from INPUT is " << temperature * ModuleBase::Hartree_to_K << " K"
-                      << std::endl;
-            std::cout << " Reading temperature from STRU is " << t_current * ModuleBase::Hartree_to_K << " K"
-                      << std::endl;
-            std::cout << " Rescale velocties to initial temperature" << std::endl;
-            rescale_vel(unit_in.nat, temperature, allmass, frozen_freedom, vel);
-        }
-    }
-    else
-    {
-        std::cout << " Random velocities according to initial temperature " 
-                  << temperature * ModuleBase::Hartree_to_K << " K"
-                  << std::endl;
-        rand_vel(unit_in.nat, temperature, allmass, frozen_freedom, frozen, ionmbl, my_rank, vel);
-    }
-}
-
 void init_vel(MDCell& mdcell,
               const bool& init_vel,
               const bool& restart,

source/source_md/md_func.h

@@ -27,26 +27,6 @@ namespace MD_func
  */
 double gaussrand();
 
-/**
- * @brief initialize the atomic velocities
- *
- * @param unit_in unitcell information
- * @param my_rank MPI rank of the processor
- * @param restart whether restart the md
- * @param temperature ion temperature
- * @param allmass atomic mass
- * @param frozen_freedom the fixed freedom
- * @param ionmbl determine whether the atomic freedom is fixed
- * @param vel the genarated atomic velocities
- */
-void init_vel(const UnitCell& unit_in,
-              const int& my_rank,
-              const bool& restart,
-              double& temperature,
-              double* allmass,
-              int& frozen_freedom,
-              ModuleBase::Vector3<int>* ionmbl,
-              ModuleBase::Vector3<double>* vel);
 void init_vel(MDCell& mdcell,
               const bool& init_vel,
               const bool& restart,

source/source_md/test/md_func_test.cpp

@@ -27,9 +27,6 @@
  *   - MD_func::read_vel
  *     - read atomic velocity from STRU
  *
- *   - MD_func::init_vel
- *     - initialize the atomic velocities
- *
  *   - MD_func::rescale_vel
  *     - rescale the velocity to the target temperature
  *
@@ -60,43 +57,6 @@ TEST_F(MD_func_test, gaussrand)
     EXPECT_DOUBLE_EQ(MD_func::gaussrand(), 0.60805637857480721);
 }
 
-TEST_F(MD_func_test, randomvel)
-{
-    ucell.init_vel = 0;
-    temperature = 300 / ModuleBase::Hartree_to_K;
-    MD_func::init_vel(ucell, GlobalV::MY_RANK, false, temperature, allmass, frozen_freedom, ionmbl, vel);
-
-    EXPECT_NEAR(vel[0].x, 9.9105892783200826e-06, doublethreshold);
-    EXPECT_NEAR(vel[0].y, -3.343699576563167e-05, doublethreshold);
-    EXPECT_NEAR(vel[0].z, 9.385130426808701e-05, doublethreshold);
-    EXPECT_NEAR(vel[1].x, -0.00017919300771203808, doublethreshold);
-    EXPECT_NEAR(vel[1].y, 5.7074002254799079e-05, doublethreshold);
-    EXPECT_NEAR(vel[1].z, -3.1088136026582953e-05, doublethreshold);
-    EXPECT_NEAR(vel[2].x, 0.000141316492668737, doublethreshold);
-    EXPECT_NEAR(vel[2].y, -0.00015841124290501442, doublethreshold);
-    EXPECT_NEAR(vel[2].z, 1.900921882689748e-05, doublethreshold);
-    EXPECT_NEAR(vel[3].x, 2.7965925764981002e-05, doublethreshold);
-    EXPECT_NEAR(vel[3].y, 0.00013477423641584702, doublethreshold);
-    EXPECT_NEAR(vel[3].z, -8.177238706840153e-05, doublethreshold);
-}
-
-TEST_F(MD_func_test, getmassmbl)
-{
-    ucell.init_vel = 0;
-    temperature = 300 / ModuleBase::Hartree_to_K;
-    MD_func::init_vel(ucell, GlobalV::MY_RANK, false, temperature, allmass, frozen_freedom, ionmbl, vel);
-
-    for (int i = 0; i < natom; ++i)
-    {
-        EXPECT_DOUBLE_EQ(allmass[i], 39.948 / ModuleBase::AU_to_MASS);
-        EXPECT_TRUE(ionmbl[i].x == 1);
-        EXPECT_TRUE(ionmbl[i].y == 1);
-        EXPECT_TRUE(ionmbl[i].z == 1);
-    }
-
-    EXPECT_TRUE(frozen_freedom == 3);
-}
-
 TEST_F(MD_func_test, readvel)
 {
     MD_func::read_vel(ucell, vel);
@@ -142,64 +102,24 @@ TEST_F(MD_func_test, RescaleVel)
     EXPECT_DOUBLE_EQ(vel[3].z, 0.00013737032325207373);
 }
 
-TEST_F(MD_func_test, InitVelCase1)
-{
-    ucell.init_vel = 1;
-    temperature = -1.0;
-    MD_func::init_vel(ucell, GlobalV::MY_RANK, false, temperature, allmass, frozen_freedom, ionmbl, vel);
-
-    EXPECT_NEAR(temperature, 300.0 / ModuleBase::Hartree_to_K, doublethreshold);
-}
-
-TEST_F(MD_func_test, InitVelCase2)
-{
-    ucell.init_vel = 1;
-    temperature = 300.0 / ModuleBase::Hartree_to_K;
-    MD_func::init_vel(ucell, GlobalV::MY_RANK, false, temperature, allmass, frozen_freedom, ionmbl, vel);
-
-    EXPECT_DOUBLE_EQ(temperature, 300.0 / ModuleBase::Hartree_to_K);
-}
-
-TEST_F(MD_func_test, InitVelCase3)
-{
-    ucell.init_vel = 1;
-    temperature = 310.0 / ModuleBase::Hartree_to_K;
-
-    EXPECT_DOUBLE_EQ(temperature, 310.0 / ModuleBase::Hartree_to_K);
-}
-
-TEST_F(MD_func_test, InitVelCase4)
-{
-    ucell.init_vel = 0;
-    temperature = 300.0 / ModuleBase::Hartree_to_K;
-    MD_func::init_vel(ucell, GlobalV::MY_RANK, false, temperature, allmass, frozen_freedom, ionmbl, vel);
-
-    EXPECT_DOUBLE_EQ(temperature, 300.0 / ModuleBase::Hartree_to_K);
-}
-
-TEST_F(MD_func_test, InitVelCase5)
-{
-    ucell.init_vel = 1;
-    temperature = 330.0 / ModuleBase::Hartree_to_K;
-    MD_func::init_vel(ucell, GlobalV::MY_RANK, true, temperature, allmass, frozen_freedom, ionmbl, vel);
-
-    EXPECT_DOUBLE_EQ(temperature, 330.0 / ModuleBase::Hartree_to_K);
-}
-
 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);
+    const ModuleBase::Vector3<double> test_velocity(0.1, 0.2, 0.3);
+    for (int i = 0; i < natom; ++i)
+    {
+        allmass[i] = 39.948 / ModuleBase::AU_to_MASS;
+        vel[i] = test_velocity;
+    }
     MD_func::compute_stress(ucell, vel, allmass, true, virial, stress);
-    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);
+    EXPECT_DOUBLE_EQ(stress(0, 0), 2.9128300667662788);
+    EXPECT_DOUBLE_EQ(stress(0, 1), 5.8256601335325575);
+    EXPECT_DOUBLE_EQ(stress(0, 2), 8.7384902002988341);
+    EXPECT_DOUBLE_EQ(stress(1, 0), 5.8256601335325575);
+    EXPECT_DOUBLE_EQ(stress(1, 1), 11.651320267065115);
+    EXPECT_DOUBLE_EQ(stress(1, 2), 17.476980400597668);
+    EXPECT_DOUBLE_EQ(stress(2, 0), 8.7384902002988358);
+    EXPECT_DOUBLE_EQ(stress(2, 1), 17.476980400597672);
+    EXPECT_DOUBLE_EQ(stress(2, 2), 26.215470600896506);
 }
 
 TEST_F(MD_func_test, dump_info)